Novel Cellulose Based Foam-Formed Products: Applications andNumerical StudiesbyPouyan JahangiriB.Sc. Mechanical Engineering, Amirkabir University of Technology, 2010A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of Applied ScienceinTHE FACULTY OF GRADUATE AND POSDOCTORAL STUDIES(Mechanical Engineering)The University Of British Columbia(Vancouver)November 2013c? Pouyan Jahangiri, 2013AbstractA novel methodology using cellulose fibres in foam laid media is proposed in order to producebiodegradable, low-density porous materials called foam-paper. Foam-forming is a process inwhich paper-making fibres are located among the bubbles created by aqueous solution of surfac-tant. Finally, the suspension is de-watered and a 3D structure of fibre-network is made. Due tothe 3D porous structure, high specific volume and capability of cellulose to participate in chemicalreactions, these new products can be applied in heat insulation, sound absorption and aerosol fil-tration. Simplicity of production, biodegradability and being economically affordable are the mostnoticeable factors which prefer foam-papers to the other products with similar applications. In thecurrent work, the first set of experiments is done on morphological characteristics of foam-papers.The effect of manufacturing condition (foam air-content) and fibre morphology such as fibre lengthand coarseness on the characteristics of the final product is studied. The influence of fibre specificsurface (by adding valley beaten fibres), fibrillation (by using PFI refined fibres) and different pulptypes is also investigated on tensile-index and specific volume (bulk) of the final products. Thesecond set of experiments is run to specify the novel applications of foam-papers in sub-micronaerosol filtration, heat insulation and sound absorption. In order to perceive the physics of theflow though the grossly disordered geometry of foam-papers and also for providing insight into thestructure of these novel products, numerical simulations based on Lattice-Boltzmann technique arecarried out. For this purpose, 2D micro X-ray tomographic images are taken from the cross-sectionof some foam-paper samples to reconstruct their 3D structure. Finally, a model based on randomcylinders and the frequency distribution of fibres in the thickness of the samples is proposed toreduce the huge amount of memory and large amount of CPU time.iiPrefaceThe current thesis is about proposing a novel methodology using cellulose fibres in foam laid mediain order to produce biodegradable, low-density porous materials called foam-paper which has takenplace in the period of 2011 to 2013.The first chapter of the study is about answering to three fundamental questions to informwhy this research was carried out. The questions are mainly about ?Why we start making foam-papers??, ?How we came up with the idea of making foam-papers?? and ?What we did in ourresearch to introduce these new products??. The need for writing this chapter is identified by Dr.D. Mark Martinez and Pouyan Jahangiri gathered the information to answer these question.In the second chapter, the major characteristics of foam-papers are studied. The research isconducted on measuring some physical and mechanical properties of the foam-papers. Dr. ArioMadani proposed the need to study the effect of fibre morphology on the properties of the finalproducts. The experiments was continued by Pouyan Jahangiri to understand the effect of fibremorphology foam-papers properties. This chapter is supervised by Dr. James A. Olson and Dr. D.Mark Martinez.In the third chapter, the applications of foam-papers on sub-micron aerosol filtration, thermalinsulation and sound absorption were investigated. Dr. James A. Olson proposed to design exper-iments to measure the foam-paper properties for these three applications. The experimental workwere carried out by Pouyan Jahangiri. The contribution of UBC aerosol lab (Courtesy of Dr. StevenN. Rogak) and UBC research centre for Acoustic Research (Courtesy of Dr. Murray R. Hudgson)were significant. This chapter was authored by Pouyan Jahangiri and supervised by Dr. James A.Olson and Dr. D. Mark Martinez.In the fourth chapter, the purpose of the chapter was to understand the physics of the flow insidethe porous media and finally to simulate the flow inside the geometry of foam-papers (Courtesy ofDr. A. B. Phillion UBC Okanagan campus). The Lattice-Boltzmann method is proposed by Dr. D.Mark Martinez for the flow simulation. Pouyan Jahangiri carried out the numerical analysis andiiideveloped a 3D unsteady solver based on the Lattice-Boltzmann method for parallel machines. Amodel geometry based on randomly oriented cylinders and fibre concentration in foam-papers wereproposed by Pouyan Jahangiri. This chapter was authored by Pouyan Jahangiri and supervised byDr. D. Mark Martinez and Dr. James A. Olson.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Foam-paper research: a brief introduction . . . . . . . . . . . . . . . . . . . . . . 32 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1 Foam structure, characteristics and rheology . . . . . . . . . . . . . . . . . . . . . 62.1.1 Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1.2 Characteristics and stability . . . . . . . . . . . . . . . . . . . . . . . . . 92.1.3 Rheology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.2 Foam-forming and its application in papermaking . . . . . . . . . . . . . . . . . . 143 Foam-paper: experiments and applications . . . . . . . . . . . . . . . . . . . . . . . 173.1 Production and characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17v3.1.1 Production details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183.1.2 Morphology and characteristics . . . . . . . . . . . . . . . . . . . . . . . 203.2 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2.1 Aerosol filtration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2.1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.2.1.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2.1.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . 323.2.1.3.1 Effect of air-content and pulp type . . . . . . . . . . . 343.2.1.3.2 Effect of VB weight ratio (fibre specific surface) . . . . 353.2.1.3.3 Effect of average fibre length . . . . . . . . . . . . . . 363.2.1.3.4 Effect of nano-fibrillated fibres . . . . . . . . . . . . . 373.2.2 Heat insulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2.2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 393.2.2.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.2.2.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . 433.2.3 Acoustics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.2.3.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.2.3.2 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 503.2.3.3 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . 514 Foam-paper: Lattice-Boltzmann simulation . . . . . . . . . . . . . . . . . . . . . . . 554.1 Kinetic theory and lattice gas cellular automata . . . . . . . . . . . . . . . . . . . 594.1.1 Kinetic theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 594.1.2 A brief review of lattice gas cellular automata (LGCA) . . . . . . . . . . . 614.2 Lattice-Boltzmann method: hydrodynamics . . . . . . . . . . . . . . . . . . . . . 624.2.1 LBM versus Navier-Stokes solvers . . . . . . . . . . . . . . . . . . . . . . 654.2.2 Applications of LBM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.2.3 Source of errors in LBM . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.2.4 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 684.2.4.1 No-slip boundary condition . . . . . . . . . . . . . . . . . . . . 694.2.4.1.1 Bounce-back algorithm . . . . . . . . . . . . . . . . . 694.2.4.1.2 Filippova-Hanel-Mei (FHM) algorithm . . . . . . . . . 704.2.4.2 Open boundary conditions . . . . . . . . . . . . . . . . . . . . . 724.2.4.2.1 Inlet/Outlet boundary condition . . . . . . . . . . . . . 72vi4.2.4.2.2 Symmetry boundary condition . . . . . . . . . . . . . 744.2.5 Numerical stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.3 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 754.3.1 Benchmark studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.3.1.1 Flow inside a 2D channel . . . . . . . . . . . . . . . . . . . . . 764.3.1.2 Flow around a 2D circular cylinder . . . . . . . . . . . . . . . . 784.3.1.3 Pressure driven flow around bundle of 2D cylinders . . . . . . . 814.3.1.4 Flow inside a uniform 3D random fibrous structure . . . . . . . . 824.3.2 Simulation of flow inside model geometry of a foam-paper . . . . . . . . . 855 Conclusions and final remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92viiList of TablesTable 3.1 Foam-Paper samples: pulp types and additives. . . . . . . . . . . . . . . . . . . 19Table 3.2 Morphological characteristics of the pulp used to make foam-paper. . . . . . . . 19Table 3.3 Pressure drop for the air-dried NBSK and CTMP samples used in 6.66 l pm . . . 35Table 3.4 Specifications of the apparatus(TPS-2500-S) which is used to measure thermalconductivity of foam-papers. . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Table 3.5 Characteristics of NBSK foam-paper samples for thermal conductivity test. . . . 44Table 3.6 Air flow resistivity of air-dried NBSK and freeze-dried NLF foam-papers com-pared to commercial sound absorbers. . . . . . . . . . . . . . . . . . . . . . . . 54Table 4.1 3D Benchmark test case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83Table 4.2 Foam-paper simulation cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . 86viiiList of FiguresFigure 1.1 Total paper consumption by Paper Grade, North America (2009) [1] . . . . . . 2Figure 1.2 Total paper consumption by Printing Grade, North America (2000-2009) [1] . . 2Figure 2.1 Image of a wet foam. (UBC Pulp and Paper Centre) . . . . . . . . . . . . . . . 7Figure 3.1 NBSK foam-paper at 40% air-content. . . . . . . . . . . . . . . . . . . . . . . 20Figure 3.2 (a) Bulk vs. foam Air-Content (b) Tensile-Index vs. foam Air-Content forNBSK and CTMP foam-papers. . . . . . . . . . . . . . . . . . . . . . . . . . 21Figure 3.3 (a,c) Bulk vs. foam Air-Content (b,d) Tensile-Index vs. foam Air-Content forvarious weight ratio of hardwood and PFI refined NBSK samples. . . . . . . . 22Figure 3.4 (a) Bulk vs. foam Air-Content (b) Tensile-Index vs. foam Air-Content forvarious weight ratio of NBSK valley beaten foam-papers. . . . . . . . . . . . . 23Figure 3.5 FQA results for (a) hardwood foam-papers (b) VB foam-papers . . . . . . . . 24Figure 3.6 Tensile-Index vs. Bulk. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 3.7 (a) Normalized vs. fibre aspect-ratio (b) FQA results for samples made bythree different average fibre length (courtesy of Dr. Ario Madani) . . . . . . . 25Figure 3.8 Porosity vs. foam Air-Content (a) different pulps (b) hardwood additive (c)valley-beaten additive. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25Figure 3.9 Gravitational settling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27Figure 3.10 Diffusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Figure 3.11 Impaction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Figure 3.12 Interception. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 3.13 Filtration efficiency versus particle diameter. (the picture is made by author.) . 30Figure 3.14 The aerosol circuit for measuring filtration efficiency and filter pressure-drop. . 32ixFigure 3.15 NBSK and CTMP (a) Normalized permeability vs. foam Air-Content (b)Darcy number vs. Solidity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 3.16 Micro X-ray tomographic image of NBSK foam-paper produced in a 50%foam air-content (Courtesy of Dr. Andre Phillion, The University of BritishColumbia Okanagan campus). . . . . . . . . . . . . . . . . . . . . . . . . . . 33Figure 3.17 (a) Permeability vs. Hardwood weight ratio (b) Permeability vs. pulp Freeness(c) Permeability vs. VB weight ratio. . . . . . . . . . . . . . . . . . . . . . . 33Figure 3.18 Filtration efficiency (a) NBSK in different foam air-contents (b) NBSK vs.CTMP pulp in 45% foam air-content. . . . . . . . . . . . . . . . . . . . . . . 34Figure 3.19 Filtration efficiency for different weight ratio of VB pulp in NBSK foam-paperstructure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35Figure 3.20 Micro X-ray tomographic images for the foam-paper samples at air-content of40% (a) 0% VB (b) 10% VB (c) 50% VB. (Courtesy of Dr. Andre Phillion,The University of British Columbia Okanagan campus). . . . . . . . . . . . . 36Figure 3.21 Effect of average fibre length on filtration efficiency. . . . . . . . . . . . . . . 36Figure 3.22 Filtration efficiency for NLF in different foam air-content. . . . . . . . . . . . 37Figure 3.23 Diameter distribution of NLF (The figure is made by author and the data arefrom: http://www.eftfibers.com) . . . . . . . . . . . . . . . . . . . . . . . . . 38Figure 3.24 Filter quality factor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Figure 3.25 Thermal conductivity of NBSK foam-paper vs. foam air-content in three work-ing temperatures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 3.26 Thermal conductivity of NBSK foam-paper vs. working temperature for dif-ferent porosities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 3.27 Thermal conductivity of NBSK foam-papers vs. temperature comparing tocommercial heat insulators. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46Figure 3.28 Three microphone impedance tube for acoustics tests. . . . . . . . . . . . . . . 51Figure 3.29 Acoustic properties: (a) NBSK (50% AC) in different thickness (b) Effect ofmulti-layer and compressed samples (c) Effect of air-content. . . . . . . . . . . 52Figure 3.30 Acoustic properties: results of freeze-dried fibrillated fibres (NLF). . . . . . . 53Figure 3.31 Acoustic properties of foam-paper products compared to commercial acousticmaterials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54Figure 4.1 Curved boundary condition and lattice orientation. . . . . . . . . . . . . . . . 71Figure 4.2 Unknown and known probability distribution functions in D2Q9 lattice. . . . . 73xFigure 4.3 (a) Lattice-Boltzmann simulation of a 2D channel flow at Re= 4.0. (b) Residuals. 76Figure 4.4 L2 norm errors of x-velocity for five different grid resolutions. . . . . . . . . . 77Figure 4.5 Errors made by applying Zou-He boundary conditions . . . . . . . . . . . . . 77Figure 4.6 Compressibility errors (a) Solution (b) RMS error vs. lattice Ma number. . . . 78Figure 4.7 Compressibility instability plots at Ma = 1. . . . . . . . . . . . . . . . . . . . 78Figure 4.8 (a) Solid-nodes which have a connection at least to one fluid-node (b) Value ofgeometric parameter ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79Figure 4.9 Drag coefficient vs. Red number. . . . . . . . . . . . . . . . . . . . . . . . . . 80Figure 4.10 Drag coefficient vs. cylinder resolution at Red = 10. . . . . . . . . . . . . . . 80Figure 4.11 Velocity magnitude around circular cylinder at low and moderate Red . . . . . . 81Figure 4.12 Normalized permeability vs. porosity . . . . . . . . . . . . . . . . . . . . . . 82Figure 4.13 Sample uniform 3D fibrous geometry. . . . . . . . . . . . . . . . . . . . . . . 83Figure 4.14 Normalized permeability vs. porosity (3D case). . . . . . . . . . . . . . . . . 84Figure 4.15 X-velocity: 3D uniform random fibrous media. . . . . . . . . . . . . . . . . . 84Figure 4.16 Tomographic images of side and top view of a foam-paper. . . . . . . . . . . . 85Figure 4.17 Fibre concentration frequency in the thickness of foam-paper. . . . . . . . . . 85Figure 4.18 Random numbers generated from the distribution of fibre frequency in thethickness of a foam-paper. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 4.19 Convergence criteria (a) Residuals (b) Reinlet (c) Mamax. . . . . . . . . . . . . 87Figure 4.20 Comparing the x-velocity simulation of a foam-paper with a uniform 3D ran-dom fibrous media at same porosity. . . . . . . . . . . . . . . . . . . . . . . . 87Figure 4.21 Convergence criteria (a) Residuals (b) Reinlet (c) Mamax. . . . . . . . . . . . . 88xiNomenclatureA Area P PressureAC Air-Content Q FlowrateC Concentration QH Heat flowrateck Discrete velocity in kth direction Q12 Collision operatorcs Lattice speed of sound ~q FluxDp Pore diameter strength ~qH Heat fluxdp Particle diameter Weight Re Reynolds numberF External force ReP Pore Reynolds numberf Probability distribution function (PDF) Red Reynolds number based on diameterf1 PDF of the first particle ~r Particle position in phase-spacef12 Two particles PDF S Entropyf123 Three particles PDF Ss Specific surfacef1? Post collision PDF T Temperaturef e Equilibrium PDF TI Tensile-Indexf ne Non-equilibrium PDF TH Hot plate temperaturef ? Superficial PDF in FHM method TC Cold plate temperaturefk PDF in kth direction Ss Specific surfacef?k Post collision PDF in kth direction T TemperatureH H function Vs Superficial velocityK Permeability VT Thermal velocityk Thermal conductivity Vrel Relative velocityL Length ~v Particle velocity in phase-spaceLt Thickness ? Dynamic viscosityMa Ma number ? Filtration efficiencym mass ? Porosityxii~n Normal unitary vectorT Tortuosity? Density? Scattering function?s Solid angelKB Boltzmann constant? Relaxation timewk Weighting function in kth direction? Geometric parameter in FHM algorithm? Extrapolation parameter in FHM algorithmxiiiAcknowledgmentsI would like to express my deepest appreciation to my supervisors, Professor James A. Olson andProfessor D. Mark Martinez who continually and convincingly conveyed a sprit of adventure in re-gard to research. Without their guidance and persistent help this dissertation would have not beenpossible.I would like to thank Dr. Steven N. Rogak and Dr. Murray R. Hodgson for providing experimentalfacilities in UBC Aerosol Lab and UBC Centre of Acoustic Research.I sincerely thank Anupam Biswas for his indispensable support, suggestions and scientific discus-sions during this project.I would like to thank Sadaf S. Zeinoddini, Banda Logawa, Yash Sharma, Amin Engarnevis, JamesMontgomery, Mina Samimi, Mostafa Mohsenvand, Dr. Reza Korei and Dr. Ario Madani for theirgreat help and support during this research.I would like to thank George Soong and Nici Darychuk for their continuous assistance in UBCPulp and Paper Centre and other members of the Centre for all the social activities that helped tomaintain a positive atmosphere.I would like to thank the Natural Science and Engineering Research Council of Canada (NSERC)and all members of Green Fiber Networks society for their help and support.Finally, I would like to thank my family and all my friends for their unconditional support.Pouyan JahangirixivDedicationTo my mother, the best person I have ever met.xvChapter 1IntroductionClimate change and pollution caused by use of non-renewable fossil fuel products has deleteriousimpact on human health, in addition, continuously increasing price of oil raises the need of substi-tuting the oil based products by green renewable ones. As the world?s environmental consciousnessand knowledge about the importance of the green products is progressing, the demand for these ma-terials is increasing. The first and the oldest eco-label certificate, is awarded by the Blue Angel,launched in Germany[1] in 1978 for services and products which are environmentally compatible.Blue Angel covers about 10000 products in 80 different categories[1]. After this revolutionary idea,many countries started issuing their own national and international certificates. Today, numerouscountries like Japan, India, Malaysia and so many other Asian countries have their own certificatesand standards for green products (up to 50 different types of different green certificates and labels).The main message of these certificates is to inform the customers about the importance of the finalproducts being environmentally compatible as it helps in local and global contamination reduction.Biodegradable substances are materials which are able to be decomposed in natural environment.One of the most important sub-branches of biodegradable materials is called green-renewablewhich can be derived from existing plants and trees. These substances can be recycled and re-produced again. Cellulose is one of the most abundant, cost effective renewable resources that canbe used as an initial substance for making green products. In 2009 printing grade paper comprisedthe largest (up to 66% of the whole pulp and paper industry) share of all paper grades consumedin the North-America followed by container-board (used for packaging applications), followed bytissues as illustrated in Figure 1.1.166%26% 8%010203040506070Printing Grade Containerboard TissueFigure 1.1: Total paper consumption by Paper Grade, North America (2009) [1]2000 2001 2002 2003 2004 2005 2006 2007 2008 2009468101214161820YearMillionsofTonnesNewsprintUncoated Freesheet (Book, Copy paper)Coated Freesheet (Direct mail, Report and ...)Figure 1.2: Total paper consumption by Printing Grade, North America (2000-2009) [1]From 2006 to 2009, total North America paper consumption in printing grades (Newsprint,Uncoated Freesheet and Coated Freesheet) decreased by 24% as depicted in Figure 1.2 . Basedon this information, it can be concluded that the profitability of traditional markets of paper grades(mainly printing grades) as a major part of North America?s forest and pulp and paper industries,2are diminishing and regaining seems improbable. The renovation in forest and pulp and paperindustries in North America (especially in Canada), will take place if the production line policiesare revolutionized. One way for getting closer to this purpose is to switch from the current paperproduction lines to more profitable and innovative ones.1.1 Foam-paper research: a brief introductionIn follow up to one of the major missions of the Canada?s strategic network for the developmentof innovative green wood fibre products, this research is conducted on production of novel eco-friendly lignocellulosic fibre-based materials in foam medium. In this process, the honeycombthree dimensional geometry of foam is used as a template in order to obtain an open structure.Since the final product has a 3D porous structure, has extremely low density and is uniquely soft,it can be applied as aerosol filters, heat insulators, impact resistors and sound absorbers.The Procedure of making foam-paper looks almost similar to regular paper-making process. Inthe process of making handsheet, first a suspension of pulp and water in a specified consistency(between 0.5% to 1.0%) is made. After that, a specified volume of pulp slurry for target grammageof 100 g/m2 based on TAPPI standard procedure is de-watered under vacuum. At this stage, fibresare deposited on one layer of forming fabric and eventually, the handsheet is oven dried. In foam-paper making procedure, the pulp slurry is prepared by adding up a percentage of standard surfaceactive agent (surfactant) under strong shear force which results in uniform bubbles in mixture. Thisuniform bubble medium plays the role of a 3D template for the fibres that are already located inbetween the bubbles. The resulting shear thinning suspension is de-watered under vacuum and thesample is air dried without pressing for 24 to 48 hours in CTH 1 room. The final product has open3D structure, high porosity, high permeability and is extremely light.In the first study, an attempt is made to comprehend general characteristics of foam-papers.For this purpose, two different pulp types, NBSK2 and CTMP are used to make the samples indifferent foam air-contents. It is generally figured out that by increasing foam air-content, bulkand porosity of samples increased which can be considered as the main advantage of the method.However, it has to be mentioned that the disadvantage of the current method is noticeable lossof some mechanical properties like tensile-index. To overcome this problem two different fibrous1Constant Temperature and Humidity2Northern Bleached Softwood Kraft3strengthening additives are added to the structure of foam-papers (valley-beaten NBSK fibres andPFI refined NBSK) and effect of fibre specific surface and freeness are studied. Improvement in re-sults of tensile-index is observed in both cases. It is shown that high specific surface valley-beatenfibres show much more improvement (up to 70% tensile-index of the standard handsheet) but bulkreduction is the penalty which is paid for this achievement. In order to be more precise about theeffect of fibre morphology on the final product characteristics, fractionation is applied and the ef-fect of different fibre length is investigated on the bulk and tensile-index of final products.In the second study, three different applications of foam-papers are studied experimentally. Thefirst one is the application of foam-papers as sub-micron aerosol filters. Pressure drop and perme-ability of foamed filters are measured using the apparatus in UBC aerosol laboratory as shownin Figure 3.14. Scanning Mobility Particle Sizer (SMPS) is also connected to the same setup tomeasure the filtration efficiency of foam-papers in sub-micron particle range (between 10 nm and600 nm). Effects of fibre type (using different pulps: NBSK and CTMP), fibre specific surface andfreeness (using NBSK valley beaten fibres and PFI refined NBSK respectively) and morpholog-ical characteristics of fibres such as fibre length and coarseness (using fractionated NBSK pulp)on filtration efficiency of foam-papers are studied. To understand the effect of crowding number,a brief study is done on some new samples of foamed filters made up of nanofibrillated fibres inpretty high foam air-contents (50% to 70%) using freeze-drying technique. Finally the results arecompared to some available data of commercial aerosol filters.The second one is the application of foam-papers in heat insulation. Samples of the same gram-mage (100 g/m2) but different porosities of NBSK foam-paper made up of NBSK pulp is usedin the experiment. Thermal conductivity of samples is measured using TPS3 standard analysismethod at three different temperatures and the results are compared to some available results ofcommercial heat insulators. The third one is the application of foam-papers in acoustics. Samplesof NBSK pulp made in three different foam air-contents are used for sound absorption coefficientmeasurement in low frequencies. The experiment is also done for multi-layer foam-papers andfoam-papers made up of nanofibrillated fibres. The results of absorption coefficients are comparedto commercial acoustic materials.In the third study, the physics of the flow inside grossly disordered geometry of foam-papersis the matter of interest. For this purpose, micro X-ray tomography technique is used to take 2D3Transient Plane Source4images of the cross-section of sample. Initially, some image processing techniques are appliedto de-noise and binarize 2D images. After that, 2D images are reconstructed to a 3D matrix formaking the actual geometry of the foam-paper. Other useful information such as density distri-bution in the thickness and porosity of foam-papers is also computed based on the results of themicro X-ray tomograph. After reconstructing the geometry, a 3D parallel flow solver based onLattice-Boltzmann method (LBM) is developed using C/C++ language on Linux platform. At thebeginning, the code is validated for some well-known 2D and 3D problems and the results arecompared to some CFD 4 results of FLUENT simulations and some available literature for low andmoderate Reynolds numbers. Then the results of permeability of the simulation for some fibrousporous geometries made by 3D random cylinders are compared to some previous LBM simula-tions and mathematical models. Finally, a model for the geometry of foam-papers based on thedata from micro X-ray tomography (density distribution function in thickness) using 3D randomcylinders is proposed and used instead of using the actual geometry for reducing the large amountof CPU time. The flow is solved for the computer generated geometry of a foam-paper and valuesof permeability are computed.4Computational Fluid Dynamics5Chapter 2Background2.1 Foam structure, characteristics and rheologyBefore discussing about foam forming process and its application in paper making, it seems ap-propriate to have a brief review on rheology and characteristics of foam, dynamics of foam, foamdrainage and foam flow in porous media.2.1.1 StructureMost of the time, the word ?Foam? is used for aqueous foam i.e. a substance that encompassesuniform or non-uniform concentrated dispersion of gas bubbles (non-continuous phase) in a volumefraction of water (continuous phase) which contains a percentage of surfactant[2]. Usually, inaqueous foams, volume fraction of bubbles called foam air-content is greater than volume fractionof water. Foams are broadly divided into ?wet? and ?dry? based on volume fraction of water. Theair-content of wet foams is normally less than 80% and the shape of the bubbles is very close to acomplete sphere as it is shown in Figure 2.1 but the air-content of dry foams is greater than 90%(the volume fraction of liquid part is less than 10%) and the shape of the bubbles is very similar topolyhedral.6Figure 2.1: Image of a wet foam. (UBC Pulp and Paper Centre)The 3D random cellular morphology of the foam is long matter of interest for mathematicians.In fact, foam mathematics is a multi-scale system in many cases. The first scale is bubble scale. Inthis scale foam is studied as ?idealized foam?, the geometry of foam is a good candidate for math-ematical problems of 3D tessellations and minimal surfaces problem. In 1887 Lord Kelvin askedhis famous question which was ?What is the most efficient bubble foam??. Indeed, dividing a 3Dspace into equi-volume cells using minimum surface area between them is the matter of questionand is referred to as ?Kelvin?s Problem?. Kelvin proposed his famous model based on bi-truncatedcubic honeycomb, which is called ?Kelvin?s Structure?[3] which is accepted as the solution for thisproblem until it is disproved by Denis Weaire, German physicist and his student Robert Phelan in1993[3] and the model was called Weaire ? Phelan structure. This model is known to be the mostoptimal unit cell of a perfectly ordered foam (ideal foam). In mid 19th century, Belgian physicist,Joseph Plateau describes the structure of soap films based on his experimental observations. Manyother different examples with almost the same pattern in nature also conform Plateau?s law[4]. Inthe second scale which is smaller than size of a bubble, the continuum liquid that splits the flat facesbetween two bubbles is called the lamella. The area that three lamella intersect is called Plateauborder of the foam. The third scale, is the scale for gas-liquid interface at the film separation sur-face which most of the time is occupied by surfactant molecules to stabilize the foam at larger scale.Major properties of foam such as stability and viscosity is specified by bubble diameter, bubblesize distribution and air-content (related to foam density). Numerous experimental and numericalstudies have been carried out on the effect of bubble size and bubble size distribution on charac-teristics of the final foam. It is important to study the effect of these parameters since they have7significant effect on the final product?s quality. For instance Isarin et al.[5] mention in their arti-cle that in textile technology, the quality of final product is strongly dependent on the bubble sizedistribution of foam. Gidoa et al.[6] perform an experimental study on measuring bubble size ofan aqueous foam in two different foam air-contents which are made by applying high shear forceusing rotary mixers. They generally measure that increasing the angular speed of the mixer causesa reduction in average bubble radii. They also observe that the mean bubble size is slightly smallerin higher air-content foam than the lower air-content one, however, no result is presented about thebubble size distribution. The experimental study on the bubble size distribution of foam which isproduced by the same technique as Gidoa et al.[6] is conducted by Kroezen[7, 8]. The author findsthat the average bubble radius decreases as the applied mean shear force increases. In follow upto this study, a non-linear relationship between bubble size distribution and mean bubble radius isfound for different foam air-contents in various flow regimes (from laminar to turbulence). Hirt etal.[9] investigate the effect of viscosity of the liquid phase of foam on the average bubble radiusfor low air-content foams. The author finds that by increasing the viscosity of liquid phase, themean bubble size diminishes especially in the foams with lower air-content. An interesting studyis done by Chang et al.[10] on measuring mean bubble diameter and its distribution by freezingsmall samples of foam. The author shows that freezing the aqueous foam (foam in solid state) willnot affect the bubble size and its distribution. So this can be considered as easier approach forinvestigating on the foam structure. Isarin et al.[5] utilized image analysis techniques in order tocompute bubble radii and bubble size distribution in an aqueous foam. This study is one of the ear-liest studies which used image analysis technique to calculate foam parameters. The author finds agood reproducibility in their results comparing to the other available experimental studies. Manyother researchers have done studies on using different imaging techniques such as Diffusing-WaveSpectroscopy (DWS), X-ray tomography, Neutron Scattering and Magnetic Resonance Imaging(MRI) to reconstruct a clear 3D image of foam structure. The internal bubble architecture hasbeen also studied computationally using sequential attempts of evolution of the minimum surfaceenergy both in deterministic (surface evolver technique[11]) or probabilistic methods (Ising modeland Potts model[12]).There are also different types of foams in which the continuous part can be a solid or can beother liquids such as liquid metal or polymeric liquid. Solid foams are made by dispersion ofgas bubbles inside a solid medium. These materials are actually solid porous structures which areconsidered as cellular light and bulky engineering materials with vast applications in different in-dustries. They are categorized into two major divisions, closed-cell foams and open-cell foams.8Closed-cell foams are those which do not have interconnected pores so they normally have moredensity and less air volume fraction. Open-cell foams are those which have interconnected poresand consequently, this category is more permeable, softer and cheaper (needs less amount of ma-terial) than closed-cell foams. Metal foams and polymeric foams are two examples of the mostimportant industrial solid foams. Metal foams are cellular structures containing a large volumegas-fraction which are made by solid metals mainly aluminum alloys. These materials have vastapplication in making orthopedic tools. Polymeric foams are materials containing a polymericmatrix with interconnected pores or air bubbles (closed pores). They are used in a broad varietyof applications such as insulation materials, disposable packaging and ultra soft materials. In thisthesis, the word ?Foam? points to wet foam which is a result of surfactant agitation in aqueousenvironment.2.1.2 Characteristics and stabilityFoam stability is capability of foam to persist for noticeable length of time without losing its 3Dstructure and collapsing into two detached constituent phases. In order to be more accurate, fromthermodynamic point of view, aqueous foams are inherently unstable and breaking down happensin the direction of diminishing total surface free energy. Surface active agents can contribute todelaying foam decay to a great extent. Without using surfactant, though bubbles can be generatedby applying high shear force on fluid, the final two phase fluid will not last long due to gravitywhich causes liquid drainage, Laplace pressure which causes mass transfer from smaller to largerbubbles and osmotic pressure which causes an internal drainage in the foam structure from lamellato Plateau borders due to local concentration difference[13]. In the first mechanism, the liquid isdischarged due to its own weight. The more fluid is draining, the thinner the lamella is becomingand eventually this will result in rupture of foam bubbles. In the second mechanism, which mostlyhappens in poly-disperse1 foams, spontaneous diffusion of gas through the interface from smallerbubbles to neighbouring larger bubbles happens[13, 14]. This phenomenon is already predicted byYoung - Laplace equation. Based on this equation, capillary pressure is inversely proportional tothe effective radius of interfacial film i.e. smaller bubbles cause higher capillary pressure whicheventually results in the inter-diffusion of gas from smaller bubbles to the adjacent larger ones. Thisprocedure spontaneously increase the average size of the bubbles in foam media (by shrinkage ofsmaller bubbles and expansion of the larger adjacent bubbles) and makes narrower foam film whichresults in bubble bursting and failure in the 3D structure of foam. In the third mechanism, since the1Non-uniform bubble size distribution9lamella is slightly curved, Plateau border has lower pressure than lamella. The pressure gradientresults in liquid drainage from lamella to Plateau border which makes the lamella very narrow andfinally the foam structure will collapse.Both bubble formation and foam stability can significantly be improved by adding surfaceactive agents (surfactants) to the liquid. Surfactants are materials that reduce the surface tension of aliquid or the interfacial tension of two liquids. Surfactants are amphiphilic i.e. their molecules haveboth hydrophilic (polar) group in the water solution and hydrophobic (non-polar) group away fromwater. Surfactant molecules are accumulated in the gas-liquid interface and contributes to stabilizethe detached dispersed gas phase in the liquid phase. The addition of the surfactant moleculesinto the interface provides an expanding force against the surface tension[15] which diminishesthe existing tension at the gas/liquid interface and finally reduces the effective surface tension.Surfactants contribute to bubble formation by decreasing the total surface energy. The energyrequired to make a change in surface area is given by:dG = ?dA (2.1)where dG is surface free energy and ? is surface tension. The greater the amount of adsorption ofsurfactant molecules, the greater the decrease in ? . A decrease in ? results in diminishing the totalsurface free energy, so bubbles need less energy to form. It should be also noted that, increasing thesurfactant concentration will not always result in decreasing surface tension. The critical micelleconcentration (CMC) is defined as the concentration of surfactant above which the surface tensionremains relatively constant due to the formation of micelles. Before reaching CMC, considerablechange in surface tension can be observed[2, 16].Surfactants can also contribute to make more stable foam. As mentioned before, a destabilizingmechanism exists due to the curvature in lamella and pressure gradient between Plateau region andlamella. The continuous liquid phase in lamella is attached to the bubble surfaces by the hydrophilicpart of surfactant molecules. Presence of surfactant in this part is very important which can reducethe surface tension in the foam lamella. Reducing the surface tension incur reducing in Laplacepressure as predicted by Young-Laplace equation.?P = ?(~?.n?) (2.2)By diminishing the Laplace pressure, the amount of fluid which would continuously drain to10Plateau border reduces. This helps the foam to become more stable. DLVO theory2[17] can explainfoam stability in more detail. Based on DLVO theory, the forces which play vital role in stabilizingfoam are Van-der-Waals attraction between the surfaces of bubbles and the electrostatic repulsionif the lamella contains ions. In the absence of surfactant (no ions in lamella) the Van-der-Waalsattractive force causes bubbles to get closer to each other and the resulting pressure drains the fluidof lamella to Plateau border. Consequently, this will lead to thinning the lamella and merging bub-bles, which eventually results in destabilizing the foam. However, increasing the concentration ofsurfactant incur increasing the charge density in lamella which increases the effect of electrostaticrepulsive forces. The final trade off between attractive and repulsive forces increases the stabilityof the final foam if the concentration of surfactant is sufficient.2.1.3 RheologyAlthough the structure and stability of foam are investigated in general literature, rheology of foamas a complex non-Newtonian fluid is also a matter of concern. By increasing the application of foamin different engineering processes, studying foam flow both from fundamental points of view[18?20] and engineering applications view point becomes more important. For instance Krug[18],Wendorff and Ainley[19] and Lincicome[20] study foam flow and rheology in different pipes andco-centric cylinders. Many other works deal with the application of foam in oil industry such asapplication of foam to move and displace oil through a granular structure. Experimental studies onfoam flow and rheology is started by Sibree[21] and Grove et al.[22] on soap-foam and fire-fightingfoam. Different types of viscometers and other non-viscometer apparatus have been setup by dif-ferent researchers to investigate various rheological parameters of foam such as fractional forcewhich delays the foam flow, viscosity and stress-strain relationship by changing the concentrationand type of surfactant and different techniques for foam preparation which has a major effect onthe rheological properties of the final foam.The rheology of foam is one of the most interesting fluid mechanics problems due to the specificliquid/gas dispersion structure. Based on the results of the works cited earlier, a single rheologicalbehaviour is not reported. The reason behind it is related to geometrical scale of bubbles which arecomparable to the geometrical scales of fluid phase and measuring apparatus.Thus, a small change in the structure of the foam (for example due to foam drainage and collapseduring the test) results in a noticeable change in the measurement of the final rheological parame-2Derjaguin-Landau-Verwey-Overbeek theory11ters, flow behaviours and finally leads to significant statistical fluctuations. The parameters whichhave main effect on the foam flow and rheology are the ratio of average bubble volume to theflow channel size, bubble size distribution, foam interaction with channel walls, type of fluid andsurfactant, isotropy of bubble distribution, absolute pressure and the flow regime (depends on theReynolds number). It should be noted that in most of surveys, final results are affected by someof the parameters that their influence is not taken into account. Rotational viscometry[14] andcapillary tube viscometry[14] are two major methods to measure the rheological characteristics ofaqueous foams. Marsden and Khan[23] use a rotary viscometer (modified fan VG meter) to mea-sure apparent viscosity of a continuous flow of foam in the range of 70% to 96% air-content. Theyfind that foam is a shear-thinning fluid and the values of apparent viscosity are between 55 cP and500 cP. They also investigate the effect of air-content and rotational speed on the foam apparentviscosity. Based on Marsden and Khan[23] results, apparent viscosity increases by increasing thefoam air-content at the same rotational speed and decreases by increasing the rotational speed inthe same air-content.Later, Wenzel et al.[24] study the effect of slip of foam flow on a smooth surface in a co-axial cylindrical viscometer and later on the cone viscometer. Wenzel et al.[24] conduct theirexperiments for numerous dried foams and plot the growing trend of shear stress by increasingthe angular velocity of viscometer. Patton et al.[25] utilize the method of capillary tube viscom-etry by studying foam motion in long capillary tubes at upstream and releasing to atmosphere atdownstream. They apply this method due to the utilization and applicability of such rheologicalresults to understand the physics of the foam flow in porous media. Patton et al.[25] find that thevalue of normalized pressure drop (with respect to inlet pressure) is close to 1. This implies thatas the foam continues flowing along capillary tubes, it also considerably expands which results intwo essential conclusions. First, bubble size distribution changes along the capillary tube (largerbubbles are observed at the outlet). Second, it is found that the foam pressure varies non-linearlyalong the capillary tube; However, drainage of the liquid part of foam leads to bubble rupture andeventually changing the rheology of foam inside capillary tubes. The latter effect is also dependson the diameter of the test tube and wall slip condition. This geometry dependency of foam rheo-logical parameters (here apparent viscosity) is investigated by many researchers[24?26]. Patton etal.[25] also figure out that the relationship between dP/L and the foam flow rate in the same tubediameter is non-linear. Hirasaki and Lawson[26] prove both theoretically and experimentally thatthe rheological characteristics of foam are strongly dependent on the diameter of the capillary tubeat the same flow regime.12Many other researchers studied foam dynamics with a different point of view. Foam can alsobe viewed as a continuous bulk non-Newtonian fluid. Different non-Newtonian fluid models forthe stress-strain relations such as power law and Bingham plastic have been tried to predict the be-haviour of foam as a continuous fluid by comparing to the experimental results which are geometrydependent. This dependency of models to geometry seems to arise due to unidentified effect of atleast one essential parameter in the experiments as discussed before. This problem is overcomeand demonstrated by Kraynic[27] by combining two different theories concerning shear deforma-tion for foam-flow inside a tube. This research motivated scientists to write more general modelsfor shear deformation of flowing foam. For a 2D ideal mono-disperse foam, Princen[28] predict ayield stress which causes a deformation. To be more precise, before reaching this value only elas-tic deformation could happen; However after reaching the critical value, the cells (represents thebubbles trapped in the fluid phase of the foam) in the row which is closer to the applied shear stressshift their location by one cell in the direction of applied shear force. This situation is dramaticallyidealized and compared to the 3D poly-dispersed structure of the actual foam but such yield stresscan also be predicted for the real foam. Princen[28] also finds a relation for the yield stress ofideal two dimensional mono-disperse foam. He finds that the values of yield stress is proportionalto the values of surface tension and inversely proportional to the values of bubble diameter. Thedimensionless factor of proportionality is a value close to unity and is strongly dependent on thefoam air-content. Hirasaki and Lawson[26] modify Princen?s[28] relation for the yield stress toreach a critical capillary number for the two dimensional ideal foam. For this purpose, they use theidea of hydraulic diameter of the wetted perimeter of the foam lamella which osculates the tubewall. They also extend the analysis of single bubble motion to motion of a bubble chain which aresplit by lamella.Another interesting topic in foam dynamics is foam flow in grossly disordered and porousstructures. This topic is of great concern in many industries especially in oil industry. The majorequation that relates the flow rate to the pressure drop of porous media irrespective of the internalgeometry of the media is Darcy?s equation. Darcy?s equation is experimentally determined byDarcy and also can be derived from Navier-Stokes equation for homogeneous media. This equationis valid for the flow in continuum regime and in the range of low to moderate Reynolds numbers.Equation 2.3 represents the Darcy equation.13~q =?k?~?P (2.3)where k is permeability, ? is dynamic viscosity, ~?P is pressure gradient and~q is the flux (dischargeper unit area). In the case of foam flow in porous media, apparent viscosity is being used inEquation 2.3 and q?P is called mobility. It is observed that, the mobility of foam flow throughporous medium for the foams with smaller average bubble size (cell size) is much lower[14]. Thereason for that is experimentally studied by Heller[14] and Hirasaki and Lawson[26]. This is dueto the fact that the apparent viscosity of foam is extremely geometry dependent[14, 26]. Hirasakiand Lawson[26] point out that by increasing the average bubble size in a capillary tube, apparentviscosity of the foam decreases dramatically which leads to increasing the foam mobility inside thecapillary tube. The same scenario happens in foam flow in porous media.2.2 Foam-forming and its application in papermakingBefore starting discussing about the application and advantages of foam-forming process in paper-making, it is wise to present a brief introduction on paper-making process. Paper is a universalproduct with vast applications mostly in writing, packaging, hygiene products and so on. Omittingthe details of different paper products, the process for manufacturing of paper is generalized insome steps.In the first step, after breaking the wood into wood chips and making fibres, the fibres which areuseful for paper-making process are separated and beaten down to pulp using various refining tech-niques. In the second step extra treatment is done on pulp. Mechanical, chemical, biological andother parameters of the final paper is controlled by adding specific chemical premixes to the pulp.In the third step, a dilute suspension of pulp (at low consistency of about 0.5% to 1.0%) in water isprepared and drained though a screen called forming fabric. Forming fabric is mostly made up ofmono-filament plastic strands and its weave pattern must be designed in a way that could guaranteethe minimum resistance to water drainage and maximum retention of fibre fines. In this step, pulpslurry is drained on a forming fabric under vacuum and a mat of randomly distributed interwovenfibres deposit on the fabric. The size of the final paper sheet depends on the size of the formingfabric. After that, extra moisture of paper is removed first by applying pressure and then drying bypassing hot air flow. The third step in paper making is carried out automatically on paper-machinesin paper-making mills.14One of the major disadvantages of paper-making process is requirement of the huge amountsof water which are used to transfer fibres from the head-box to the paper-making section in paper-machine. The more amount of water is used for fibre transformation, the more water needs to bedrained and consequently more energy is used for the water removal process, which is significantat the industrial scale. One way to reduce this amount of energy is to switch the fibre suspendingmedium from water to a fluid which can allow less energy consumption for its drainage. Aqueousfoam is introduced to pulp and paper industry in 1972 as a good alternative to water[29]. Sincemore than 70% of foam is air, it needs much less energy to drain the suspension in forming sec-tion. The second important advantage of foam is the 3D honeycomb structure which prevents fibreflocculation and also provides the flexibility of making porous, bulky and light papers with openstructures for making new products with specific applications. This happens because the fibreslocate in between the bubbles during the process.The idea of applying foam in paper-making initially is proposed by Radvan and Gatward[29]in 1972 and Smith and Punton[30] at 1974. Radvan and Gatward[29] for the first time proposed theidea of using foam as a medium for suspending fibres in paper-making. They find that a uniformdispersion of very long and fine fibres even at relatively high concentrations can be made with-out fibre flocculation. They also propose a discontinuous foam forming unit attached to a smallpaper machine and they called it Radfoam process[29]. Tringham[31] applies Radfoam processand explores different new properties of papers made by this process. He finds that papers whichare made by Radfoam process are typically very uniform and have more open structures comparedto the normal papers. Tringham[31] also finds that The specific volume (bulk) of the hand-sheetsusing Radfoam process is 20 to 30 times greater than a standard handsheet for softwood kraft pulp.Smith et al.[30, 32] use Radfoam process to make new papers in foam media and studied othercharacteristics of the new products. They confirm that foam forming process is a good method tomake uniform paper (without fibre flocculation) using longer fibres. They conclude that surfacetension and bubble spacing in foam change the properties of the final material however chemicaleffects are negligible. Moreover, It is showed in their research that beaten pulp increase the tensilestrength of foam paper due to increasing the specific surface of fibres.Recently, Lappalainen and Lehmonen[33] at VTT Technical Research Centre of Finland studycharacteristics (like bubble size distribution) and dynamics of foam containing cellulose fibres.They use a high speed imaging technique to collect various images of foam-fibre mixture (foam15bubbles with radius larger than 25?m can be recognized in their imaging technique). Image pro-cessing techniques are combined with Circular Hough Transform for determination and analysisof bubble size distribution of the foam-fibre suspension. Poranen et al.[34] at VTT Technical Re-search Centre of Finland produce new cellulose based 3D porous products in foam environmentusing an automatized foam handsheet maker in a pilot scale foam forming research environmentcalled SUORA. They use various types of pulps and different basis weights to make their products.They study some morphological characteristics of the foam formed materials and also present anautomatized foam paper maker. To our knowledge the last two works done in VTT are the initialworks which have been seriously followed and published in the last few decades.16Chapter 3Foam-paper: experiments andapplications3.1 Production and characteristicsThe essential concern of paper-making is all about uniformity of final paper. The source of non-uniformity is related to fibre flocculation in the pulp slurry. Fibre flocculation is the consequence ofnumerous parameters but the most essential ones are fibre length, consistency of pulp suspensionand the viscosity of suspending medium. One way of circumventing this problem is to make a verydilute aqueous suspension of pulp however, energy consumption and production rate limitationsregarding water drainage, prohibit using high-dilution processes. Another approach to solve thisproblem is replacing the aqueous suspending medium by foam. Since 60 to 70 percent of foam thatis being used for this purpose is air, much lower amount of water needs to be drained comparingto the high-dilution processes. At the same time, since fibres are located in the foam lamella, fibreflocculation will be prevented even in case of relatively long fibres. Another advantage of usingfoam is that since the structure of foam is 3D and porous, it is a good candidate to be utilized asmedium to make 3D porous paper called foam-paper with specific industrial applications. Beforestudying various applications, lets initially start explaining the production detail and characteristicsof foam-paper.173.1.1 Production detailsIn this study, a method very similar to Radfoam process[29] is applied to make the foam suspen-sion of fibres. The suspension of foam and fibre is prepared using a foam maker unit based onthe apparatus similar to the one reported by Smith et al.[30, 32]. The foam maker unit has twomixing cell part and each cell has a Caframo Brinkmann RZR1 type mixer to apply shear forcein order to reach higher foam air-content1. In the first mixing cell, normally 5 ml of surfactant isadded to 1000 ml of the suspension but for reaching higher foam air-content the concentration ofthe surfactant is increased. The foam air-content is varied from 10% to 50% by changing the mixerangular velocity from 280-2200 RPM2. The higher air-contents from 50% to 70% are reached byincreasing the concentration of surfactant. Between the first cell unit and the second cell unit thereis a vent at the bottom to let the foam pass from first cell to the second one. In the second cell, thefoam bubbles get smaller and more uniform since extra shear force is applied and also cellulosefibres are added to the foam simultaneously. The consistency of the pulp solution that are usedfor making foam-papers was 0.5% based on the TAPPI recommendation for standard handsheetmaking. Figure3.1 shows a sample of a final product.Two different methods of drying are applied to make different samples in different foam air-contents for various purposes and also study the effect of different drying techniques on the finalproducts. The first method is air-drying without pressing. In this method, the suspension of foam-fibre is then filtered in a 15 cm diameter Buchner funnel under 9.8 kPa of vacuum. One layer ofregular fabric and two layers of plastic forming fabric is used inside the funnel for making thesamples. De-watered samples are then dried without pressing at 24?C and 56% humidity in CTH3room for 24 to 48 hours. In this method of drying, two different pulps are used. Some additivesare also used for special purposes which will be explained in the next section. The second methodis freeze-drying. In this method, after making the foam-fibre suspension, it will be frozen in liquidNitrogen (?196?C) and de-watered using a freeze-drier apparatus. Freeze-drier works by reducingthe pressure inside the container to be sublimated from solid phase to gas phase immediately. Thedesired grammage of both air-dried and freeze-dried samples are 100g/m2. Table 3.1 shows thedetail of the pulp and additives used the make samples for different applications.1Air-Content is the volume fraction of air in foam. It can be measured using a cylinder in lab by subtracting the totalvolume of foam with the fluid volume per total volume of foam.2Round Per Minute3Constant Temperature and Humidity18Table 3.1: Foam-Paper samples: pulp types and additives.Pulp type Additives Drying method AC(%)NBSK ? Air-Dried 10%-55%NBSK 10%-50% Eucalyptus(Hardwood) Air-Dried 10%-55%NBSK 10%-50% NBSK Valley Beaten Air-Dried 10%-55%NBSK 100% PFI refined Air-Dried 20%-60%CTMP ? Air-Dried 10%-55%NBSK ? Freeze-Dried 20%-70%NLF4 ? Freeze-Dried 50%-70%In order to be more precise about the effect of fibre morphology on the foam-papers properties,fractionation5 method is applied to separate longer fibres from the shorter ones. The effect of fibrelength and coarseness are studied on foam-papers properties and aerosol filtration tests. For thispurpose, fibres are filtered based on their length in different stages of Bauer-Mcnett fractionator.The screen sizes in this fractionator are 16, 30 and 50 based on ASTM standard (labelled as LWM16 ,LWM30 and LWM50). The length distribution, average fibre coarseness and the average fibre lengthare analysed using FQA6. Morphological characteristics of the fibres that are used to make thefoam-papers are briefly reported in Table 3.2. However more detailed information about the lengthdistribution of fibres will be reported in the next section.Table 3.2: Morphological characteristics of the pulp used to make foam-paper.Pulp type LW (mm) d f (?m) Coarseness (mg/m) Fine(%)NBSK 2.45 27.6 0.17 2.50CTMP 1.8 30.3 0.52 6.15NLF 0.1-0.3 0.2-0.3 ? ?Eucalyptus 0.76 ? 0.04 2.72LWM16 - NBSK 3.29 28.7 0.19 25.4LWM30 - NBSK 2.72 27.8 0.14 29.8LWM50 - NBSK 1.81 26.6 0.15 ?3Nanofibrillated Lyocell Fibre5Fractionation is a technique to separate long fibres from short fibres.6Fibre Quality Analyser193.1.2 Morphology and characteristicsThe morphological characteristics of the fibres that are used to make the foam paper, are studied.The research in this section is conducted in finding the relation between the manufacturing condi-tion and the final-product condition. Moreover, the effect of fibre morphology on the final-productproperties is investigated. In this research, foam air-content is considered as manufacturing con-dition and parameters such as specific volume (bulk), porosity and tensile index are considered asfinal-product condition.The first experiment which is actually one of the main purposes of this research, is done onmeasuring the specific volume (bulk) of foam-paper samples. Figure 3.2a shows that by increasingfoam air-content, the values of bulk dramatically increases, up to 50 times for NBSK and up to 20times for CTMP, compared to the standard handsheet. CTMP shows increase in specific surfaceand this can be ascribed to the presence of the fines. Longer fibres create an open structure andfines are placed in this open space and since there has not been a very noticeable amount of water,the hydraulic forces during the drainage do not transfer the fines. This effect contributes to increaseof the specific surface and consequently lower bulk at higher air-contents. Generally, the penaltyfor increasing bulk is losing the mechanical properties (such as tensile-index) to a great extent.Figure 3.2b shows that by increasing the air-content up to 55% the values for the tensile-indexdiminish down to about 10% of standard handsheet for CTMP and about 4% of standard handsheetfor NBSK pulp. CTMP has higher value of tensile-index due to the presence of high specific areafibres compared to NBSK.Figure 3.1: NBSK foam-paper at 40% air-content.20??????????????CTMP10 20 30 40 50 60102030405060Air-ContentH%LHBulkL?HBulkLHandsheet? NBSK(a)10 20 30 40 50 6010?210?1100Air?ContentM(%)TI/TI HandsheetNBSKCTMP(b)Figure 3.2: (a) Bulk vs. foam Air-Content (b) Tensile-Index vs. foam Air-Content for NBSKand CTMP foam-papers.In order to solve the strength problem, different weight percentages of refined Eucalyptus (hard-wood) fibres is added to the structure of NBSK foam-papers. Figure 3.3b shows that by increasingthe weight ratio of hardwood in foam-paper structures, tensile-index increases slightly which alsoleads to a slight reduction in the values of bulk. It should be noted that a laboratory PFI refinerwith 5000 revolution is used to refine hardwood Eucalyptus pulp. The freeness of hardwood pulpafter refining is reduced from 417.7 CSF to 335 CSF. PFI refining increase the percentage of exter-nal fibrils and fines that contributes to increasing the tensile-index of foam-papers. However bothhigher amounts of fibrils and shorter fibres lead to lower values of bulk (Figure 3.5a).To study the effect of fibrillation on both bulk and tensile-index more precisely, PFI refiningapplied to NBSK pulps and samples with three different values of freeness are made from newrefined pulps. Figures 3.3c and 3.3d show the effect of freeness on bulk and tensile-index respec-tively. It can be concluded that by decreasing the values of fibre freeness (increasing the refiningenergy), the values of tensile-index increases but this is quite opposite for the values of bulk.21000????????????@@@?10% Eucalyptus030% Eucalyptus@50% Eucalyptus10 20 30 40 50 60102030405060Air-ContentH%LHBulkL?HBulkLHandsheet? NBSK(a)10 20 30 40 50 6010?210?1100Air?ContentE(p)TI/TI HandsheetNBSK10pEEucalyptus30pEEucalyptus50pEEucalyptus(b)??????????????10 20 30 40 50 60102030405060Air-ContentH%LHBulkL?HBulkLHandsheet? NBSK HCSF = 760L? PFI refined NBSK HCSF = 155L? PFI refined NBSK HCSF = 518L(c)10 20 30 40 50 6010?210?1100Air?ContentF(%)TI/TI HandsheetNBSKF(CSFF=F760)NBSKF(CSFF=F518)NBSKF(CSFF=F155)(d)Figure 3.3: (a,c) Bulk vs. foam Air-Content (b,d) Tensile-Index vs. foam Air-Content forvarious weight ratio of hardwood and PFI refined NBSK samples.To study the effect of specific surface more rigorously, three different weight ratio of NBSKvalley beaten (VB) fibres are added to the structure of standard foam-papers. As shown in Figure3.4b the tensile-index of the samples improved up to 70% of standard handsheet for 50% weightratio of VB samples. It is not surprising that the values of bulk is reduced by factor of 2.5 inthe same sample. The more amount of VB fibres in the structure of a foam-paper lead to thehigher average surface area of the fibres and consequently stronger bondings. At the same timethe average fibre length is reducing (Table 3.2). Since longer fibres are able to support the 3D22structure, increasing the percentage of shorter fibres in the structure of samples leads to diminishingthe values of bulk (Figure 3.5b). Samples containing different weight ratios of VB fibres are provedto be more effective at increasing mechanical properties of the final foam-papers than hardwoodsamples as illustrated in Figure 3.4. Since preserving high bulk and high tensile-index is a tradeoff, samples containing 30% of VB fibres show reasonable bulk and tensile-index simultaneously.The trade off between tensile-index and bulk for various foam-papers is illustrated in Figure 3.6.???????????????????10% Valley -Beaten?30% Valley -Beaten?50% Valley -Beaten10 20 30 40 50 60102030405060Air-ContentH%LHBulkL?HBulkLHandsheet? NBSK(a)10 20 30 40 50 6010?210?1100Air?ContentV(%)TI/TI HandsheetNBSK10%VVB30%VVB50%VVB(b)Figure 3.4: (a) Bulk vs. foam Air-Content (b) Tensile-Index vs. foam Air-Content for variousweight ratio of NBSK valley beaten foam-papers.Figure 3.5 shows the average fibre length distribution for both hardwood and VB foam-papers.For samples with additives two maximum points is observed which shows the presence of twofrequent fibre length in the mixed samples. It can be concluded that the average length for valleybeaten fibres are much lower than the hardwood ones. Other morphological parameters (results ofFQA) are reported in Table 3.2.230 1 2 3 4 5 600.511.522.533.544.55Lw (mm)Frequencyo(%)NBSK10%oHardwood50%oHardwood(a)0 1 2 3 4 5 600.511.522.533.544.555.5Lw (mm)Frequencyt(%)NBSK10%tValleyt?tBeaten50%tValleyt?tBeaten(b)Figure 3.5: FQA results for (a) hardwood foam-papers (b) VB foam-papers?????@@@000??????????????????????????10% NBSK Valley-Beaten?30% NBSK Valley-Beaten?50% NBSK Valley-Beaten?10% Eucalyptus030% Eucalyptus@50% Eucalyptus?CTMP0 20 40 60 80 100 12005101520Bulk Hg?cm3LTensile-IndexHN.m?gL? PFI refined NBSKHCSF = 518L? PFI refined NBSK HCSF = 155L? NBSKFigure 3.6: Tensile-Index vs. Bulk.To understand the effect of fibre morphology on the properties of final product, fractionatedNBSK pulp is used to make samples with different average fibre length. Figure3.7a shows thatby increasing the fibre aspect-ratio, the values of normalized bulk increase at different foam air-24contents at the same grammage. This happens because much more space is required for the samenumber of longer fibres to locate among each other than shorter ones.70 80 90 100 11010203040506070Fibre aspect-rat io ( L Wdf )Bulk/Bulk Handsheet30% Air ? Content40% Air ? Content50% Air ? Content(a)?? ? ?????????? ? ?????????????????????????????????????????????????????????????????????????????? ??? ? ? ? ?? ???????????????????????????? ??????? ?????? ??? ? ???? ???????? ????????????????????????????????????????????????????????????????????? ? ?????? ??? ? ? ???????? ??? ???????????????????????????????????????????????????????????????????????????????????????????????? ?? ????????????? ?????? ????????????????????????????????????????????????????????????????????????? ??????????????????? ??LWM16?LWM30?LWM500 1 2 3 4 5 60.00.51.01.52.02.53.0LWHmmLFrequencyH%L? NBSK(b)Figure 3.7: (a) Normalized vs. fibre aspect-ratio (b) FQA results for samples made by threedifferent average fibre length (courtesy of Dr. Ario Madani)The other important parameter is porosity of final product. porosity is defined as the ratio ofvolume of pores to the total volume of sample. Figure 3.8 shows how porosity of final productschanges by changing the foam air-content. Figure 3.8 generally shows that foam-papers are materi-als with very high porosity. By increasing foam air-content the porosity of final products increasesup to 99.7% for NBSK foam-paper.20 25 30 35 40 45 50 55 600.930.940.950.960.970.980.991Aird?dContentd(%)PorosityNBSKCTMP(a)30 35 40 45 50 550.960.9650.970.9750.980.9850.990.9951AirS?SContentS(%)Porosity NBSK10%SEucalyptus30%SEucalyptus50%SEucalyptus(b)30 35 40 45 50 550.960.9650.970.9750.980.9850.990.9951AirV?VContentV%Porosity NBSK10%VValleyV?VBeaten30%VValleyV?VBeaten50%VValleyV?VBeaten(c)Figure 3.8: Porosity vs. foam Air-Content (a) different pulps (b) hardwood additive (c)valley-beaten additive.253.2 ApplicationsIn this section three potential industrial applications of foam-papers are studied. The 3D porousand filter-like structure of foam-papers let these materials to be applied in filtration of sub-micronaerosol particles. The pressure drop of these products are noticeably lower than the standard hand-sheet and also the filter quality factor of some samples are comparable to the commercial filters.In this section, first the pressure drop, permeability and filtration efficiency of foam-papers aremeasured and then the effect of different additives on filtration efficiency and filter pressure dropis studied. Since the porosities of foam-papers are quite high, they can also be applied as heatinsulators and sound absorbers some of which can compete with their commercial counterparts.The major purpose of this study is to figure out how the essential target parameters vary withdifferent manufacturing conditions. It is also important to study the effect of various additives onthese parameters and make efficient foam-papers for that specific purpose.3.2.1 Aerosol filtrationAn air filter is known as a medium composed of fibrous materials which is able to capture and re-move airborne particles. The importance of using aerosol filters is raised by increasing the demandfor clean air. Controlling the air quality plays a vital role in several industries such as buildingindoor air quality (for HVAC7 purpose) and internal combustion engines. Sometimes increasingthe concentration of particulates in the air causes very serious health problems, so it is important tomake very efficient filters to capture airborne particles in micro and nano scales. Since the internalstructure of foam-papers is made up of random orientation of entangled fibres and their specificvolume is high (which are associated with low pressure drop), we can assume foam-papers arecandidates for aerosol filtration applications. In this section, a research on filtration properties ofvarious foam-papers is carried out.3.2.1.1 BackgroundA non-woven fibrous filter media is commonly used to capture airborne particulates. Fibres in thesize range between 100 nm and 100 ?m criss-cross to form a web of several layers. This structuremostly contains air in order to increase breathability (Reduce pressure drop). Since the particlestend to go through very small holes in this porous structure, they may not become trapped. Rather,while trying to navigate through the layers of filter media, a particle will be attached to a fibre dueto a numerous different mechanisms. The main mechanisms are inertial impaction, gravitational7Heating, Ventilation and Air Conditioning26settling, Brownian diffusion, electrostatic deposition, and interception[35, 36]. In order to explainhow a particle deposits in the filter media, it is necessary to first consider the flow of air stream.The path of the air around a fibre is described in terms of flow streamlines. Deviation of particlescarried by the air from the streamlines depends upon the aerodynamic diameter of particles. Thelarger the size, the greater the tendency to deviate from air streamlines due to higher inertia of aparticle. Several authors investigated that the rate of particle deposition is minimum for the parti-cles of an intermediate size (between 100 nm to 200 nm)[35, 36]. The deposition mechanisms aremore efficient for either very small or very large particles.Gravitational settling or sedimentation is settling of particles that fall down due to the grav-ity. This is the dominant mechanism governing the deposition of very large particles particularlythose in the range greater than few microns, while Brownian diffusion has the greatest effect onvery small particles less than a hundred nano-meter in size[37]. Brownian motion is the processby which aerosol particles move randomly due to collisions with the gas molecules. Smaller parti-cles randomly move across streamlines due to Brownian motion until they touch surface of a fibre.Figures 3.9 and 3.10 show the filtration of very large and very small particles due to gravitationalsettling and diffusion mechanisms respectively.FiberGravityStreamlines Gravitational SettlingFigure 3.9: Gravitational settling.27FiberStreamlines DiffusionFigure 3.10: Diffusion.Most of respirable particles are much smaller to be filtered out by gravitational settling. Res-pirable particles above 400 nm in diameter are not likely to be influenced by diffusion mechanismdue to very small diffusion coefficients. They may be captured by interception and inertial im-paction [38]. Inertial impaction is the deposition mechanism while particles are not able to followthe curved streamlines due to their high values of inertia. Impaction of a certain particle with anobstacle can be predicted by Stokes number. Stokes number is a non-dimensional group whichdepends on particle size, velocity, and drag force exerted by flow, and characteristic size of theobstacle. For the values of Stokes number greater than 1, the particle will collide with fibres[37].Interception occurs when small particles follow the streamlines, but the streamlines will naturallytake the particle close enough to the obstacle to contact the surface of the fibre. If the particle flowstoo close to a fibre, collision may happen. Figures 3.11 and 3.12 show the filtration mechanismsfor impaction and interception respectively.FiberStreamlines Inertial impactionFigure 3.11: Impaction.28FiberStreamlines InterceptionFigure 3.12: Interception.Electrostatic deposition is an entirely different particulate capture mechanism. By startingproducing resin-wool filter in 1930, many of the filters which used in different applications areelectrostatically charged[39]. Although it may be difficult to quantify the charge on either theparticle or the fibre, electrostatic attraction can be an extremely effective capturing mechanism.A charged particle will be attracted to fibre which has an opposite charge. Figure 3.13 showsthe total filtration efficiency of a typical fibrous filter media as a function of particle diameter.As seen in this figure, the efficiency is a strong function of particle size. Air filters are highlyefficient at removing particles with a diameter less than 0.1?m (by diffusion) and greater than0.6 ?m (by Impaction). Particles with diameters between 0.1 to 0.6 ?m have significantly lowerfiltration efficiency. Particles in this range are too large to be effectively pushed around by diffusionand too small to be effectively captured by impaction. It should be noted that filtration efficiency(which is measured from the concentration of particles before and after a filter), filter pressure-drop(the difference between the pressure in upstream and downstream) and filter quality-factor (takesboth the effect of efficiency and pressure drop into account) are three essential parameters that isimportant to be measured for a filter as it is shown in Equations.? =Cupstream?CdownstreamCupstream(3.1)QF =? ln(1??)?P(3.2)29Diffusion Capture Diffusion and Interception CaptureImpaction and Interception CaptureFigure 3.13: Filtration efficiency versus particle diameter. (the picture is made by author.)In pulp and paper industry, some researchers are trying to make efficient filters out of pulp inorder to find a low-cost and biodegradable replacement for polymeric filters. One of the earliestresearch works in this field is conducted by Rodman and Lessmann[40] in 1988. They reviewdifferent types of filter media that are commonly used in automotive filtration systems. They alsopresent the filter-papers which their micro-structure constructed by different types of fibres. Theeffect of different parameters is investigated on the filtration efficiency of the filters and some tra-ditional design approaches for improving the efficiencies are discussed. They did filtration test forpaper-filters in the particle diameter range between 5 and 30 ?m and concluded that adding micro-fibres in the structure of paper-filters can increase filtration efficiency up to 45%. For the sake ofthe filtration of sub-micron particles, Mao et al.[41] develop a method to make very efficient filters(up to 95% efficiency) using hardwood and softwood pulps with a pressure drop in the range of thecommercial filters. They add various percentage of high specific surface valley-beaten fibres intothe structure of paper-filters and use freeze-drying technique to retain the surface fibrillation intothe dry state. The filter structure made up of this method has very high filtration efficiency(effectof microfibrils) and low pressure drop (long NBSK fibre provide high permeability). Macfarlane etal.[42] use similar method to what Mao et al. use and carried out a comprehensive research on thethe effect of different freez-drying techniques and various parameters of paper-filter on the filtra-tion efficiency and pressure-drop of filters. They find a new parameters in their freeze-dried filters30called filter sectioning and layering. Duo to what they found, a non-uniformity in the fibre packingdensity is observed in the z-direction of filter-papers because of the effect of cold wall. This non-uniformity is increased by increasing the percentage of valley-beaten pulp in the structure of thepaper-filters.3.2.1.2 MethodologyThe filtration efficiency and filter pressure drop of samples are measured in a small aerosol circuitas it is shown in Figure 3.14 . Each foam-paper is cut into the 35 mm diameter circular sample tofit into the sample-holder in the circuit. Aerosol particles are generated by the solution of SodiumChloride (NaCl) containing 0.1% NaCl in the water using a collision nebulizer. The aerosol par-ticles is dried by passing them into a tunnel full of humidity absorber materials (desiccant) andfinally dragged through the filter sample using a normal vacuum pump. The flow-rate of the airthrough the samples are set into 6.66 L/min which is equal to the superficial gas velocity of 11.5cm/s. The pressure-drop of the samples are measured by sensing the pressure signal before andafter filter using a digital pressure sensor. The concentration of particles before and after filter (andconsequently filtration efficiency) is measured using scanning mobility particle sizer (SMPS)[36].The aerodynamic diameter of the particles using for all filtration tests are between 14 nm to 670nm8.In this survey, the application of the filtration of different foam-papers (variations are due tovarious pulp types and two methods of drying) is studied. The effect of fibre morphology and foamair-content on the filtration properties and air permeability of the final products are investigated.The experiments are started by measuring the filtration efficiency and pressure-drop of air-driedNBSK foam-papers (standard foam-paper) in different foam air-contents. The effect of fibre lengthand coarseness on the filtration properties of foam papers studied by using air-dried NBSK foam-papers made up of fractionated NBSK pulp. Filters made by air-dried CTMP pulp are used tostudy the effect of pulp type on the filtration properties of the products. Freeze-drying technique isapplied to make new samples out of Nano-fibrillated Lyocell Fibres (NLF) to find out how nano-fibres change the foam-papers structure and consequently filtration properties of the new products.8http://dustmonitors.ru/d/68562/d/spms 3936.pdf31AirCollisionNebulizerSampleHolderSMPSPressureDropPumpSMPSFigure 3.14: The aerosol circuit for measuring filtration efficiency and filter pressure-drop.3.2.1.3 Results and discussionFigure 3.15a shows that by increasing foam air-content, permeability of both air-dried NBSK andCTMP samples increases up to 9000 times greater than the standard handsheet at the highest air-content. This high values of permeability of samples are reasonable because of the very large openpaths among the fibres which are illustrated in the tomographic images. CTMP samples show lowerpermeability because of the higher specific surface and more percentage of fines in its structure.?????????CTMP20 30 40 50 601002005001000200050001 ? 104Air- ContentH%LHKL?HKLHandsheet? NBSK(a)10?210?210?11001011021??K/d2 fNBSKCTMP(b)Figure 3.15: NBSK and CTMP (a) Normalized permeability vs. foam Air-Content (b) Darcynumber vs. Solidity.32Figure 3.16: Micro X-ray tomographic image of NBSK foam-paper produced in a 50% foamair-content (Courtesy of Dr. Andre Phillion, The University of British ColumbiaOkanagan campus).0 10 20 30 40 5010?610?5Hardwood (%)Permeability (cm2 )50% AC30% AC(a)0 200 400 600 80000.511.522.5x 10?5Freeness (CSF)Permeability (cm2 )50% Air ? Content(b)0 10 20 30 40 5010?710?610?510?4VB (%)Permeability (cm2 )50% AC30% AC(c)Figure 3.17: (a) Permeability vs. Hardwood weight ratio (b) Permeability vs. pulp Freeness(c) Permeability vs. VB weight ratio.However, permeability of samples with refined Eucalyptus (hardwood) additives are reducedas is shown in Figure 3.17a. This happens due to two main reasons. First, the average length ofEucalyptus fibres are lower than the average length of NBSK ones (Figure 3.5a). Using longerfibres generally results in samples with higher permeability since longer fibres need more spacethan shorter fibres to make a fibrous network in the same grammage. Second, the percentage offibrils in refined Eucalyptus is much more than the NBSK pulp. Figure 3.17b shows the effect offibrillation on permeability of final samples. The results show that by decreasing the fibre freeness(increasing the amount of fibrils in the structure of foam-paper), permeability of samples decreases33to a great extent. Permeability reduction is much more pronounced for the samples with variousvalley-beaten weight ratio as illustrated in Figure3.17c. This happens due to the presence of veryhigh specific surface fibres in the structure of these series of foam-papers. The higher the specificsurface of fibres, the more they can block the fluid flow and consequently the less permeability willbe obtained.0 100 200 300 400 500 600020406080100ParticleCdiameterC(nm)EfficiencyC(%)NBSK (45%CAC)NBSKC(10%CAC)(a)0 100 200 300 400 500 600020406080100ParticleCdiameterC(nm)EfficiencyC(%)NBSK(45%CAC)CTMPC(45%CAC)(b)Figure 3.18: Filtration efficiency (a) NBSK in different foam air-contents (b) NBSK vs.CTMP pulp in 45% foam air-content.3.2.1.3.1 Effect of air-content and pulp typeFigure 3.18a shows the results of filtration efficiency in the range of sub-micron particles for NBSKfoam-papers at 10% and 45% air-contents. The results show that by increasing the air-content, thefiltration efficiency of standard foam-paper diminishes until it loses its filtration properties at 45%air-content. Though NBSK sample at 10% air-content show reasonable filtration efficiency, itshows much higher pressure drop per unit length which leads to having poor filter quality factor.For this reason, CTMP is used for obtaining better filtration efficiency at higher air-contents. CTMPpulp shows better filtration properties at the same air-content compared to NBSK pulp as shown inFigure 3.18b. Though the values of pressure drop per unit length of CTMP sample at 45% AC arebetter than NBSK sample at 10% AC, it still shows very low quality factor. The values of pressuredrop per unit length are reported in Table 3.3 for these test.34Table 3.3: Pressure drop for the air-dried NBSK and CTMP samples used in 6.66 l pmSample Pressure drop (Pa) Thickness (mm) dPL (Pa/m)NBSK (10% AC) 25.80 1.6 16125.0NBSK (30% AC) 11.94 4.1 2912.2NBSK (45% AC) 10.4 9.3 1120.1NBSK (54% AC) 9.77 11.4 857.20NBSK (60% AC) 8.76 15.5 565.50CTMP (45% AC) 125.4 10.9 11504.00 100 200 300 400 500 600020406080100ParticleCdiameterC(nm)EfficiencyC(%)NBSKC(45%CAC)50%CVBC(45%CAC)30%CVBC(45%CAC)10%CVBC(45%CAC)MERVC14MERVC8Figure 3.19: Filtration efficiency for different weight ratio of VB pulp in NBSK foam-paperstructure.3.2.1.3.2 Effect of VB weight ratio (fibre specific surface)Figure 3.19 shows the effect of VB weight ratio on the filtration efficiency of foam-papers. Byincreasing the weight ratio of NBSK VB pulp in the structure of foam-paper, the filtration efficiencydramatically increases. This happens due to the presence of high specific surface fibres. Micro X-ray tomographic images are also confirm the mentioned claim as illustrated in Figure 3.20. Theresults also compared to the results of the commercial filters, MERV-8 and MERV-14. According35to the pressure drop per unit length of the VB samples, the results of 10% VB and 30% VB arecomparable to the results of MERV-8 and MERV-14 respectively. 50% VB samples show a veryhigh pressure drop and very low filter quality factor.(a))(b) (c)Figure 3.20: Micro X-ray tomographic images for the foam-paper samples at air-content of40% (a) 0% VB (b) 10% VB (c) 50% VB. (Courtesy of Dr. Andre Phillion, TheUniversity of British Columbia Okanagan campus).0 100 200 300 400 500 600020406080100ParticleWdiameterW(nm)EfficiencyW(/)f=W1.81Wmm, ?f =W0.14Wmg/m)=W0.15Wmg/m)W=W2.50Wmm, ?NBSKW(45/WACW?LNBSKW(45/WACW?LWFigure 3.21: Effect of average fibre length on filtration efficiency.3.2.1.3.3 Effect of average fibre lengthFigure 3.21 compares the effect of average fibre length on filtration efficiency of NBSK samplesat almost the same air-content (45%), grammage, fibre coarseness9 and average fibre diameter10.90.14 mg/m for shorter fibres and 0.15 mg/m for longer fibres1026.8 ?m for shorter fibres and 27.6 ?m for longer fibres36The results represents higher filtration efficiency for the samples made up of shorter fibres. Sincethe grammage of the samples are the same, the average number of the fibres per unit volume ofthe samples made by shorter fibres are greater than the longer ones. This effect contributes toincreasing the tortuosity of the internal structure, hence the probability of particle-fibre collisionswill increase. Moreover, using longer fibres in the structure of foam-papers leads to having largeropen paths as elaborated before. The larger the internal void paths in the structure of foam-papers,the less probability of particle-fibre collisions. Though the filtration efficiency improved usingshorter fibres, a slight increase in the pressure drop of the samples is observed. Another problemof foam-papers made by shorter fibres is that the mechanical properties of these samples are verypoor. They easily lose their structure at high flowrates.0 100 200 300 400 500 600020406080100ParticleAdiameterA(nm)EfficiencyA(%)NanofibrillatedALyocellAFibreA(35%AAC)NanofibrillatedALyocellAFibreA(55%AAC)MERVA14Figure 3.22: Filtration efficiency for NLF in different foam air-content.3.2.1.3.4 Effect of nano-fibrillated fibresIn this section, nano-fibrillated cellulose is added to the foam-media to make new foam-paperproducts out of nano-fibrillated fibres. Freeze-drying applied to slurry of NLF and foam in order topreserve the fibre-foam internal structure and avoiding hornification. Figure 3.22 shows fantasticfiltration efficiency results for the samples made by NLF in two different foam air-contents. Thereason for this can be answered by looking into the morphology of nanofibrillated fibres. Themain difference between a actual NBSK fibre and a nanofibrillated cellulose fibre is in the values37of length and diameter of fibres as shown in Table 3.2. The average diameter of NLF is about300 nm (Figure 3.23) and the average length is between 0.1 mm to 0.3 mm. This means thatthe number of nanofibrillated fibres per unit volume is much more than the products made up ofNBSK macrofibres. Moreover, the large aspect ratio of fibres leads to a very high tortuous fibrenetworks as elucidated before. These two effects together leads to capturing between 80% to 95%of sub-micron particles in different air-contents. The pressure drop for the samples are a bit highcompared to the commercial filters (about 150 Pa for 55% air-content and about 320 Pa for 35%air-content). However, the results of filter quality factor illustrated in Figure 3.24 are reasonablecompared to MERV-14. The lower air-content samples of freeze-dried NFL are not recommendedfor the filtration purpose because of high values of pressure-drop (For air-content 20% the pressure-drop is more than 750 Pa). Figure 3.24 shows that the freeze-dried foam-papers made by NLF in35% and 55% foam air-contents and air-dried foam-papers with 30% valley beaten fibre weightratio have reasonable quality factor comparing to the commercial filters. Although the pressuredrop values for both freeze-dried NLF foam-papers are high but since these filters have very highfiltration efficiency, the quality factors of them are in the rang of commercial filters. Concerningthe 30% VB air-dried filters, having low values of pressure drop (about 50 Pa) and the averagefiltration efficiency of about 45% (between 10 nm and 670 nm particle diameter) incur a reasonablequality factor.300 600 900 12005101520NLF diameter (nm)Frequency (%)Figure 3.23: Diameter distribution of NLF (The figure is made by author and the data arefrom: http://www.eftfibers.com)380 100 200 300 400 500 60000.511.522.53ParticleVdiameterVRnm8QualityVFactorVR1/cmVH 2O830CVVBNLFVR35CVAC8NLFVR55CVAC8MERVV14MERVV8Figure 3.24: Filter quality factor.3.2.2 Heat insulationIn this section heat insulation properties of foam-papers are studied. Thermal conductivity is mea-sured in a range of working temperature for each foam-paper. Moreover, The relation betweenmanufacturing condition(foam air-content) and thermal conductivity of the final product is studied.Finally, the results of the thermal conductivity of foam-papers made using various air-contents arecompared to the results of commercial heat insulators.3.2.2.1 BackgroundThermal conductivity testing is a material property measurement associated with the difficulty inwhich heat is conducted through a specific type of material. Thermal conductivity is a pure materialproperty independent of other different parameters like area of conduction or the thickness of the39material. Thermal conductivity is defined as the ratio of heat flux to the associated thermal gradient.~qH =?k(~?T ) (3.3)In one dimension the equation can be simplified as following. This measurement can be studied asthe thermal conduction between two parallel, isothermal surfaces of area ?A? and at temperaturesTH (temperature of hot plate) and TC (temperature of cold plate) separated by a layer of the samplematerial having a thickness L with a steady state heat power Q. Thermal conductivity, k, is thusdefined as:k =Q.L(TH ?TC)A(3.4)Owing to recent developments in geological sciences, soil and building technology and oil indus-try, heat transfer in disordered structures and conductivity of porous materials are becoming veryimportant. As an example, experiment done by Tseng et al.[43] can be pointed out. The thermalconductivity of the polyurethane foam is investigated both theoretically and experimentally for thedevelopment of liquid hydrogen storage tanks in temperature range of 20 K to 300 K in their re-search.The effect of porosity of this sort of materials plays an essential role in their thermal conductivity.The relation between porosity and thermal conductivity is studied by several researchers. The au-thors try to find a relation between thermal conductivity coefficient and porosity of different porousmaterials. However some difficulties arise since thermal conductivity of porous materials is not afunction of porosity alone. Other parameters such as pore shape and type of the material of thesolid state (metallic or polymeric or ... ) play important roles in the heat transfer mechanism. Thementioned mechanism is also different for a particular porous material at high and low porosities.So that the model which is accurate for a particular material with low porosity may not be rigorousfor the same material with high porosity.Several models have been proposed for approximating the thermal conductivity of porous materialsversus porosity. The simplest model which is also quite rigorous for porosities less than 0.1 andgreater than 0.9 is proposed by Loeb[44]. More precisely, for a mono-disperse porous materialwith closed pores having the porosity less than 0.15, MaxwellEucken[45] model is recommended.This model is not valid for higher values of porosity, since the connectivity of pores (open poros-ity) change the mechanism of heat transfer. For the high porosity condition, Russells relation[46]is recommended. Russell?s relation has been developed based on parallel/series resistor approachused to explain the effect of cubic pores in a medium[46]. Reichenauer et al.[47] finds that pore40size is also a very important factor that affects thermal conductivity of porous materials. If thepore size is very small (less than 500 nm), thermal conductivity of such a material is less than thethermal conductivity of the same material (with the same porosity) with larger pore sizes. Thisphenomenon happens due to Knudsen effect.Recently, Smith et al.[48] carried out a complete experimental and numerical research about thethermal conductivity of monolithic porous materials (tin oxide, alumina, and zirconium ceramics)with both high and low porosities. They found an accurate relation between thermal conductivity ofa porous material and porosity of the same material based on the thermal conductivity of the solidphase and the gaseous phase (air). Smith et al.[48] also figured out that the thermal conductivitymodels for porosities less than 0.65 and higher than this value should be different in order to haveaccurate results.3.2.2.2 MethodologyThere are number of possible ways to measure thermal conductivity, each of them being suitable fora limited range of materials, depending on thermal properties, working temperature and other dif-ferent factors. Steady-state and transient methods are two different classes for measuring thermalconductivity of a sample material. In general, steady-state techniques perform the measurementwhile the temperature of the sample material doesn?t experience any serious variation with respectto time. This makes the signal analysis much more straightforward because of implying a quasi-constant signal. The only difficulty is that a well-designed and engineered experimental setup isusually needed. Transient techniques perform measurement in the condition that the temperatureis changing with respect to time (for example in the process of heating up the sample). This mea-surement is relatively faster compared to the steady-state one since reaching steady-state conditionis not required. However, the time dependency of the final signal with respect to time makes themathematical analysis of the final signal more difficult. The following methods are well known intransient measurement:? Transient Plane Source (TPS) and Modified TPS (MTPS) methods.? Transient Line Source (TLS) method.? Laser Flash Method (LFM).? Time-domain thermo-reflectance method41In this research we use Transient Plane Source (TPS) technique for calculation of the thermal con-ductivity of NBSK foam-papers at three different temperatures. The samples were measured atfreezing point of water (0?C), boiling point of water (100?C) and at the standard atmospheric tem-perature (20?C) in different porosities. The ThermTest- TPS-2500-S Thermal Constants Analyseris the instrument chosen for all bulk thermal conductivity measurements. The TPS-2500-S meetsthe ISO Standard (ISO/DIS 22007-2.2)11.Transient Plane Source Method, utilizing a plane sensor and a mathematical model describingthe heat conductivity, combined with electronics, enables the method to be used to measure thermaltransport properties. It covers thermal conductivity range of at least 0.01-500 WmK (in accordancewith ISO 22007-2) and can be used for measuring various kind of materials, such as solids, liquids,pastes and thin films. In 2008, it was approved as an ISO standard for measuring thermal transportproperties of polymers. This TPS standard also covers the use of this method to test both isotropicand anisotropic materials. The TPS technique typically employs two samples halves, in betweenwhich the sensor is sandwiched. Normally the samples should be homogeneous, but extendeduse of transient plane source testing of heterogeneous material is possible with proper selection ofsensor size to maximize sample penetration. This method can also be used in a single-sided config-uration, with the introduction of a known insulation material used as sensor support. By recordingtemperature versus time response in the sensor, the thermal conductivity, thermal diffusivity andspecific heat capacity of the material can be calculated. Today, many companies like Thermo TestInc. work on producing very accurate machines based on TPS method. One of the very accuratemachines for this purpose is ThermTest-TPS-2500-S, the product of Thermo Test Inc..The Transient Plane Source TPS-2500-S constitutes an industry benchmark thermal conductiv-ity test system. The TPS-2500-S thermal conductivity instrument is designed for analysing thermaltransport properties of solids, liquids, paste and powders, including various types of geometries anddimensions. Table 3.4 shows the specifications of TPS-2500-S.11http://www.thermtest.com/42Table 3.4: Specifications of the apparatus(TPS-2500-S) which is used to measure thermalconductivity of foam-papers.Materials Solids, liquids, powders and pasteThermal Conductivity Range 0.001 to 1000 W/mKThermal Diffusivity Range 0.1 to 100 mm2/sSpecific Heat Capacity Range Up to 5 MJ/m3KMeasurement Time 1 to 1280 secondsReproducibility Typically better than 1%Accuracy Better than 5%Temperature Range ?160?C to 1000?CISO Standard ISO/DIS 22007-2.2The basic principle of the TPS System is that the sample surrounds the TPS sensor in all di-rections and the heat generated in the sensor freely diffuses in all directions. The solution to thethermal conductivity equation assumes the sensor is in an infinite medium, so the measurementand analysis of data must account for the limitation created by sample boundaries. Preliminarymeasurements are on each sample. This included selection of the TPS sensor to be used, test timeand output of power to the TPS sensor. A standard wait period of 10-15 minutes is implementedbetween successive measurements on each sample to ensure the samples are isothermal prior mea-surement. Samples for this test are in circular shape and the diameter of each sample is 35mm.3.2.2.3 Results and discussionSamples of the same grammage (100g/m2) of standard foam-paper made up of NBSK pulp andBiotergeAS-40 surfactant are used for this test. Table 3.5 shows some physical and morphologicalcharacteristics of the samples used for this test. Different porosities are because of the different air-content(manufacturing condition) of the foam. The bulk thermal conductivity of the samples madein 21%, 34%, 43% and 54% foam air-contents are measured at three different temperatures (0, 20and 100 degrees of Celsius) as mentioned before. Bulk Thermal Conductivity was measured usingthe TPS Standard Analysis Method. This method is used to determine the bulk thermal conductivityof porous materials in the wide range of thermal conductivities (0.001W/mK to 1000W/mK) andtemperatures as shown in Table 3.4.43Table 3.5: Characteristics of NBSK foam-paper samples for thermal conductivity test.Foam air-content (%) Bulk (cm3/g) Porosity (%)21 20.15 96.634 46.54 98.143 50.90 98.754 95.21 99.5The results of thermal conductivity versus foam air-content at different temperatures and ther-mal conductivity versus temperature for different porosities are illustrated in Figure3.25 and Fig-ure3.26 respectively. The results show that by 3% increasing in porosity of samples, the values ofthermal conductivity decrease by about 24% in low temperature and about 15% in high tempera-ture range. This is an interesting result which shows the considerable effect of porosity in thermalconductivity reduction in order to make light thermal insulators out of pulp. As the porosity ofthe samples increases, the volume fraction of air in the equi-size samples increases. This leads todecreasing the thermal conductivity of the final product since the thermal conductivity of air (about0.026 at 25?C) is much lower than thermal conductivity of NBSK fibres (about 0.056 at 25?C).The results also show that by increasing the working temperature from 0?C to 100?C the valuesof thermal conductivity increase since the air convection through the pores (small scale channels)considerably increases. The error bars of the experimental data using TPS-2500-S show outstand-ing repeatability and reproducibility in measurements. Figure 3.5 also shows that by increasing theworking temperature, the thermal conductivity of foam-papers increasing. This is reasonable sincefoam-papers are porous materials i.e. they consist of a solid phase (fibres) and fluid phase (air). Thevalues of the thermal conductivity of the porous materials depends on both thermal conductivityof solid phase and fluid phase. The thermal conductivity of fluid phase increases by increasing theworking temperature which incurs increasing the thermal conductivity of the final material.4420 25 30 35 40 45 50 550.0350.040.0450.050.0550.060.065Air Content (%)ThermalConductivity(W/m.K)T = 0 CelsiusT = 20 CelsiusT = 100 CelsiusFigure 3.25: Thermal conductivity of NBSK foam-paper vs. foam air-content in three work-ing temperatures.0 20 40 60 80 1000.030.0350.040.0450.050.0550.060.065Temperature (C)ThermalConductivity(W/m.K)Porosity = 96.6%Porosity = 99.5%Figure 3.26: Thermal conductivity of NBSK foam-paper vs. working temperature for differ-ent porosities.45In the next step the results of the thermal conductivity of our samples are compared to thevalues of thermal conductivity of commercial heat insulators. The results show that as the porosityof foam-paper increases, the thermal conductivity values of foam-papers get closer to commercialheat insulators. Figure3.27 also shows that the values of the thermal conductivity of the foam-paper made up of 54% foam air-content are in the range of commercial insulators. To summarize,we developed porous papers demonstrated thermal conductivity close to commercial insulatorshowever more research must be done for the commercializing purposes.0 20 40 60 80 1000.030.0350.040.0450.050.0550.060.0650.07Temperature (C)ThermalConductivity(W/m.K)Polyurethane Foam (PUR)Expanded Polystyrene (EPS)Extruded Polystyrene Foam (XPS)NBSK (54% Air Content) NBSK (21% Air Content)Figure 3.27: Thermal conductivity of NBSK foam-papers vs. temperature comparing to com-mercial heat insulators.3.2.3 AcousticsIn this section, attention is focused on measuring the sound absorption (acoustic) properties offoam-papers. Air dried NBSK and freeze-dried Nanofibrillated Lyocell Fibre (NLF) are used tomake foam-papers for acoustic tests. The samples are made in both 50% and 70% air-contents inorder to study the effect of manufacturing condition on their final acoustic performance. To studythe acoustic properties of foam-papers, sound absorption and sound transmission are measured forfoam-papers with two different porosities (two different foam air-contents) in the range of low-frequencies. To improve the performance of air-dried NBSK samples, the effect of using multi-46layer foam-papers is investigated. Other non-acoustical parameters such as air-flow resistivity,tortuosity, viscous characteristic length and thermal characteristic length are also measured. Atthe end of this chapter, the results of foam-papers are compared to the results of some availablecommercial acoustic materials and also compared to some available mathematical models.3.2.3.1 BackgroundAll materials have the capability of absorbing sound to some extent but an ?acoustical? materialor ?sound absorber? is a material which is able to absorb most of the sound energy striking itssurface. This characteristic can be explained by sound absorption coefficient in a specific frequencyrange which is defined as the percentage of the incident sound energy which has been absorbed.Acoustical materials have vast industrial applications such as sound absorbers in buildings. Theycan be categorized in three major groups:? Porous absorbers: A thick solid that contains cavities, interstices and channels so that soundwaves can propagate through them and be converted to heat. Common porous absorbers arefibreglass and generally foams like Polyurethane foam.? Panel absorbers: Non-porous, thin materials which are usually placed over an airspace, thatvibrate in a flexural mode in response to sound pressure exerted by adjacent air molecules.Best examples are thin wood panelling traps at orchestra platforms.? Resonator absorbers: Perforated materials which have openings capable of absorbing soundin a limited frequency range. A classic example is the Helmholtz resonator used to identifythe various frequencies or musical pitches.Porous absorbers are the most popular and efficient materials in acoustical applications. They gen-erally have a high sound absorption coefficient compared to other materials. In porous materials,essentially closed pores do not contribute in absorbing sound since they do not have a continuousinteraction path for friction and viscous forces which are substantially responsible for momentumloss of sound waves. This phenomenon eventually result in sound attenuation as reported by Are-nas and Crocker[49]. They explain the mechanism of sound absorption as following. When aporous material is exposed to incident sound waves, the air molecules on the surface of the mate-rial and inside the pores, are forced to vibrate. During the time that the vibration is happening, airmolecules lose some of their initial energy and simultaneously convert it to heat (viscous losses onthe walls of the interior pores). This mechanism is different at low and high frequencies. At lowfrequencies, the mechanism is isothermal but at high frequencies the same mechanism is adiabatic.47In fibrous materials, much of the energy can also be absorbed by scattering from the fibres and bythe vibration caused in the individual fibres. However in granular materials or rigid porous ones thesound absorption mechanism is mainly due to the viscosity of the air trapped among the granules.Many experimental works and mathematical modelings have been done by researchers on mea-suring and modeling the acoustic properties of porous materials. A standard technique to charac-terize these absorbing materials is the impedance tube. In 1976 Seybert and Ross[50] proposeda setup based on random excitation in a two microphone tube to measure the acoustic propertiesof some available sound absorbers. They claim that the two-microphone random-excitation tech-nique can be used to evaluate the acoustic properties very rapidly since traversing is necessaryand random excitation is used. Lauriks et al.[51] observe that only free-field technique for theperpendicular surface impedance can give reliable results. They prove that the friction betweenthe screen of the porous sound absorber and the tube causes extra sound attenuation. Allard etal.[52] study on the acoustic properties of resonant materials. A very interesting research is doneby Vigran et al.[53] on the effect of constraints at the tube wall on the absorption coefficient of anelastic porous material. They compare the results of standing wave tube and free field method andfigure out that in the case of those porous materials which are influenced by the wall constraints ofthe tube, the standing tube method should be used carefully since both impedance and absorptioncoefficient could only be measured inaccurately. Champoux and Stinson[54] carry out an experi-mental research on the acoustical properties of some simple porous materials. They also comparedtheir results for some available theoretical models between 500 Hz and 4500 Hz. They showedthat their results are in good agreement with available theoretical models. Moreover, Champouxand Stinson[54] also study the effect of pore shape factor on the acoustical properties of porousmedia both in low and high frequencies. Their results show that a single shape factor leads to verygood agreement with the exact solutions for all frequencies. Champoux and Stinson[55] investi-gated on the acoustical properties of porous ceramics both experimentally and theoretically. TheAttenborough[56] and the Biot-Allard[52, 54] models are used for the theoretical modeling partof their research. They show that the experimental results of porous ceramics can totally be pre-dicted by Biot-Allard model. The method of three-microphone impedance-tube is first proposedand developed by Izumi et al.[57] to improve the two-microphone impedance-tube. By using thismethod, more acoustical parameters can be measured. This method is improved by Salisuo andPanneton[58] later. Recently, Doutres et al.[59] proposed an experimental technique based on athree-microphone impedance tube setup to measure acoustic and non-acoustic properties of poroussound absorbing materials. Based on their new approach, they measure the acoustic properties of48some porous materials such as sound absorption coefficient, sound transmission loss, effective den-sity and effective bulk modulus. In addition to that, they calculate the non-acoustic parameters suchas static airflow resistivity, tortuosity, viscous and thermal characteristic lengths using the measureproperties by an indirect approach. For calculation of the non-acoustic properties, they used thevalues of the open porosities of the porous materials.Besides the experimental studies, lots of mathematical models are proposed by many authorsto predict the acoustical properties of porous materials. The final purpose of these mathematicalmodels is to derive the characteristic wave impedance (Zc) and characteristic propagation constant(k) of a plane wave in a porous material (generally in any sound absorbing media) as functions ofnon-acoustical parameters. Two approaches can generally be used to model the sound propagationphenomenon in porous absorbers. The first one is empirical approach. The most widely used workfor this approach is the work done by Delany and Bazely[60]. However this model has some weakpoints in predicting the acoustic properties of some porous materials with high flow resistivity inthe range of low frequencies. In order to fill this gap, many other researchers such as Miki[61]contributes to Delany and Bazely model. Miki model compensates for the weak points of Delanyand Bazely model. The second approach is purely analytical approach based on micro-scale studyof the wave equation under specific assumptions for sound propagation in porous media. Biot[62]and Attenborough[56] propose two different models based on the second approach to predict theacoustical properties of a rigid porous sound absorber. The only problem in this sort of modelsarises when one deals with very complicated internal micro structures. In these structures the finalresults are highly dependent on how one defines the geometrical parameter called ?shape factor?.Besides these two major approaches, some simple phenomenological approaches are proposed toignore some complicated details which do not violate the essential physics behind the problem.The models proposed by Johnson et al.[63] and Allard-Champoux[64] are in this category.Some parameters have significant effects on the acoustical results of porous materials. Theseparameters are investigated by numerous authors. Koizumi et al.[65] finds that smaller fibre di-ameters in the internal structure of a porous material, leads to higher value of sound absorptioncoefficient. This effect is because using finer fibres in the structure of a porous medium resultsin having more tortuous medium with the same porosity. The effect of fibre diameter is also in-vestigated by Lee and Joo[66] in a separate research. Ren and Jacobsen[67] studied the effect ofair flow resistivity on the acoustical properties of porous medium. They figure out that as the airflow resistivity per unit thickness of the porous material increases, the characteristic impedance49and the coefficient of sound absorption increases. Conrad[68] shows that the current effect hap-pens because when sound wave propagates inside a material with higher resistivity against flow,its amplitude diminishes by friction and its energy of propagation is converted to heat. Anotherimportant parameter is open porosity. Shoshani et al.[69] show that by increasing the number ofpores in the direction of sound propagation, the coefficient of sound absorption of the porous ma-terial increases. They also state that the size and the shape of the pore have noticeable effect onsound absorption. Ibrahim and Melik[70] carry out a research about the effect of material thicknesson sound absorption coefficient. Coates and Kierzkowski[71] demonstrate that the effective soundabsorption in porous media happens when the thickness of the porous material is at least one tenthof the wavelength of the incident sound. Density is another important parameter and its importanceis more than other factors as it is directly related to the price of porous material. Koizumi et al.[65]show that porous materials with lower density are very good sound absorbers in low frequenciesand materials with higher densities are good sound absorbers at high frequencies.3.2.3.2 MethodologyTo measure the acoustic properties of foam-papers, an experimental setup based on the three-microphone impedance-tube method is used in the centre for acoustics researches at the Universityof British Columbia. Figure 3.28 shows the experimental setup that is used for these tests.This method is carried out by first measuring two different transfer functions. The first one is thetransfer function between the first microphone and the rigid end (the location of third microphone)and the second transfer function is that between the second microphone and the rigid end. By usingthese two transfer functions, the one between the first and the second microphones can be calculatedand eventually the material?s impedance and characteristic propagation constant will be calculated.Based on the criteria specified by ISO CD 10534-2, the results from the impedance tube that weused to calculate the acoustic properties in low frequency ranges, are valid from approximately 250Hz to 1750 Hz. The samples for this purpose is made in circular shape with diameter of 10 cm.50Figure 3.28: Three microphone impedance tube for acoustics tests.The air flow resistivity of the samples is measured using a simple experimental setup similar towhat is used before to measure pressure drop in filtration tests. For this purpose, ASTM C522-03standard is used to design the experimental setup for the test. The experimental setup induces air-flow through the sample so that the extra pressure drop can be measured due to the presence of thesample. The slope of the air flow-rate versus pressure drop graph is a function of air flow resistiv-ity. Air flow resistivity can also indirectly be calculated based on the results of the impedance-tube.Sometimes, a great difference between the results of these two methods are reported for somematerials.3.2.3.3 Results and discussionIn this section the effect of thickness, air-content and crowding number12 on the acoustic propertiesof foam-papers are investigated. The results are also compared with commercial sound absorbers.Figure 3.29a shows the variation of absorption coefficient in the low frequency range for threedifferent thickness of NBSK samples in 50% air-content. By increasing frequency, the values ofabsorption coefficient increases similar to other typical sound absorbers. It is also shown that at aparticular frequency, the absorption coefficient of foam-papers is dependant on its thickness. Thelarger the thickness of the foam-paper is, the greater the absorption coefficient in that particularfrequency will be. In fact, the amount of wave absorption is strongly dependant on the particlevelocity of that wave inside the porous material. The higher particle velocity leads to the higher12Crowding number is defined as the expected number of fibres in a sphere of diameter one.51amount of friction and energy loss and consequently higher absorption coefficient. Hence, thehigher frequency in the particular thickness provides higher particle velocity and as a result, higherabsorption coefficient. The greater thickness at a particular frequency also leads to more energywhich can be damped and consequently higher absorption coefficient. That?s the reason that multi-layer samples of foam-papers shows better sound absorption performance as shown in Figure 3.29b.Figure 3.29c shows the effect of air-content on the sound absorption performance of foam-papers.Higher air-content leads to having lower absorption coefficient at the same thickness. Higher air-content means lower number of fibres per unit volume and consequently less capability of dampingwave energy.10300.20.40.60.81FrequencylLHz=AbsorptionlCoefficienttttNBSKlL50PlACl?lL = 14mm)NBSKlL50PlACl?lL = 24 mm)NBSKlL50PlACl?lL = 30 mm) (a)10300.20.40.60.81FrequencylLHz=AbsorptionlCoefficient Compressedthree?layers Compressedtwo?layerstttNBSKlL50PlACl?lL = 28 mm)NBSKlL50PlACl?lL = 37 mm)NBSKlL50PlACl?lL = 14 mm) (b)10300.20.40.60.81FrequencylLHz=AbsorptionlCoefficient tttNBSKlL50PlACl?lL = 24 mm)tNBSKlL50PlACl?lL = 30 mm) NBSK L70PlACl?lL = 24 mm)NBSKlL70PlACl?lL = 30 mm) (c)Figure 3.29: Acoustic properties: (a) NBSK (50% AC) in different thickness (b) Effect ofmulti-layer and compressed samples (c) Effect of air-content.5210300.10.20.30.40.50.60.70.80.91FrequencyCPHzMAbsorptionCCoefficientNLFCP70=CACC?CLt tt=C60CmmMNBSKCP70=CACC- L tNBSKCP70=CACC- L = 30 mm) NBSKCP70=CAC -CL = 60 mm) =C24CmmMFigure 3.30: Acoustic properties: results of freeze-dried fibrillated fibres (NLF).Figure 3.30 compares the results of freeze-dried NLF with air-dried NBSK ones at the sameair-content and different thickness. Since sample at the thickness of 60 mm is not available, Miki?smodel[61] applied to predict its absorption coefficient. The results show much better sound ab-sorption performance for NLF samples. This happens because of large crowding number and veryhigh tortuous medium which contributes to high wave energy damping.Figure 3.31 shows the comparison between the results of absorption coefficients of some foam-paper products with commercial sound absorbers. The results show that the absorption coefficientsof foam-paper products in low frequency are in the range of commercial acoustic products. The airflow resistivity of the samples are also in the range of commercial sound absorbers at almost thesame thickness range. Table 3.6 shows the results of air flow resistivity of foam-paper samples andcompared them with commercial sound absorbers.5310300G20G40G60G81FrequencyL5HzBAbsorptionLCoefficientttCottonLBaffleSemi?RigidLGlassLFibreAcousticLFoamNBSKL550SLACL?LL = 30 mm) tNBSKL5CompL50SLACL?LL = 37 mm) NLFL570SLACL?LL = 60 mm)Figure 3.31: Acoustic properties of foam-paper products compared to commercial acousticmaterials.Table 3.6: Air flow resistivity of air-dried NBSK and freeze-dried NLF foam-papers com-pared to commercial sound absorbers.Sample Thickness (mm) Air flow resistivity ( Pa.sm2 )NBSK (50% AC) 14 2285NBSK (50% AC) 30 2936NBSK (50% AC - 2 stacks) 28 5738NBSK (50% AC - 3 stacks) 37 6913NBSK (70% AC) 30 2175NLF (70% AC) 60 28669Cotton baffle 25.4 52000Semi-Rigid Glass Fibre 25.5 21000Acoustic foam 50 2200054Chapter 4Foam-paper: Lattice-BoltzmannsimulationOn of the most usual and interesting phenomena in nature is fluid flow inside grossly disorderedgeometries. The most common example of this transport phenomenon is the flow of water insidesoil. Studying this problem is also important for many industries such as oil industry, membranetechnology, fuel cell technology, aerosol filtration technology, tissue engineering and various stagesof technical processes in different parts of process industry.For instance in oil industry, in order to extract the trapped oil in sedimentary rock, another fluidis pumped through the media to force the oil out. In aerosol technology, the goal is to capturesub-micron particles by letting them pass through a fibrous porous media called aerosol filter. Theimportance of improving our information and understanding of such processes arises because ofthese type of applications. In the first application (oil industry), the amount of energy consumedto reach the final goal satisfactorily, is high. In the second application (aerosol filtration), filtrationefficiency is strongly depend on the dynamics of fluid in filter media.One of the industries in which comprehending the dynamics of flow in porous structures hasa tremendous importance, is pulp and paper industry. The final product of this industry is mainlypaper which is made up of compact random orientation of cellulose fibres. The dynamics of fluidin different stages of paper-making process which control the final properties and the quality ofpaper, is related to the flow characteristics of porous media. For example, fibre deposition on theforming fabric or de-watering the paper in the last stage of paper-making process are good repre-sentatives of a very complex fluid flow through a porous structure. Recently, novel developments55in pulp and paper industry along the line of producing porous 3D papers (foam-paper), furtherincreased the significance of this subject. As we know by now, foam papers have applications indifferent technologies such as aerosol filtration, heat insulation and acoustic technologies. Thougha lot of experimental research has been done on characteristics and applications of foam-papers,the mysterious physics of fluid flow through grossly disordered structures of foam-papers has roomfor research research.In both theoretical and experimental research works which have been done by now on thephysics of flow through porous media, the first and the trickiest part of the research is to find amathematical correlation between the macroscopic characteristics of the porous material such asporosity and the main macroscopic flow properties. The first study which resulted in finding anempirical relation, is done by Darcy. Darcy, found a linear relation between the flow-rate and thepressure-drop of a porous medium. The slope of the flow-rate versus pressure-drop is a function ofa parameter called permeability which can be considered as the indicator of the fluid?s conductivitythrough porous structure. In laminar flow regimes, permeability is defined as the coefficient oflinear response of the fluid to the non-zero pressure gradient if the flux is induced[72, 73]. Darcy?slaw is more generalized for the turbulent flow regimes by contribution of Forchheimer[73, 74].Another important correlation is the relation between permeability and porosity. The well-knownequation for that is empirical Kozeny-Carman equation which relates permeability of the media toporosity, specific surface and tortuosity as represented in Equation 4.1.K =(1??)2? 3S2T (4.1)Other mathematical models are also proposed by many authors[75?77] for this purpose whichwill be discussed in the next section. Another important correlation is between friction factor andReynolds number called Ergun equation (Equation 4.2).(?PL)(DP?V 2s)? 3(1??) =150ReP+1.75 (4.2)where ReP =DPVs?(1??)? is defined as Reynolds number.Though the mathematical models are accurate to some extent to predict the behaviour of fluid flowin porous media but the physics is incomplete due to variation in internal geometry, the pore size,shape and scale.56Improvement in computers in terms of both CPU and memory, computational algorithms suchas parallel and asynchronous programming architectures and numerical techniques in computa-tional fluid dynamics (CFD) have made it possible to solve very complicated fluid dynamics andheat transfer problems. Recently very huge and geometrically complex fluid dynamics problemsare solved utilizing GPU1 programming. Beside improving the computational facilities and pro-gramming techniques, the vital role of fast numerical algorithms for computational fluid dynamicspurposes should not be forgotten.One approach for solving the flow equations is using conventional CFD techniques such asFVM2 to solve the Navier-Stokes equations. Fluid flow is determined by its macroscopic parame-ters like velocity and pressure and traditional CFD techniques normally solves them in a discretedomain. For instance in FVM, the domain (which is a complex geometry) should be discretized infinite volumes (which are much simpler geometries) and finally solve the continuum media equa-tions using a numerical algorithm. Though, traditional CFD techniques are considered as a goodapproach for the vast range of fluid dynamics problems but some serious limitations make onenot to consider them for solving the flow problem in random porous structures. The first reasonis that the results of CFD techniques are very sensitive to both grid resolution and the method ofmesh generation. Practically, one of the most challenging tasks of the traditional CFD approachis generating an appropriate mesh in order to reach a correct solution. Sometimes the geometryis very complicated and randomly structured and it is very time consuming to generate a propermesh to simulate the flow correctly and sometimes it is totally impractical. The other limitation ofFVM methods is that in the rarefied flow regimes, in cases that the mean free path of gas moleculesis considerable in comparison to the characteristic length of the control volume, the NavierStokesequations fail to simulate the correct physics of the flow. This phenomenon happens in simulationof flow through some porous structures such as filters made up of micro and nano-fibres. In thiscase, continuity assumption breaks down and Navier-Stokes equations are not valid.Another approach for flow simulation is from the microscopic-level point of view. In principle,the modeled flow behaviour can be modeled and simulated based on the direct modelling of thedynamics of fluid molecules. This method is called direct molecular dynamics (MD) and is com-putationally very expensive. Beside the MD approach, the Lattice Boltzmann Method (LBM) is1Graphical Processing Unit2Finite Volume Method57very promising in terms of ease to simulate the dynamics of fluid and effectiveness is being dealtwith distributed computational domains (parallel processing architecture and platform). The mod-els based on LBM are quite different from both accurate microscopic description (MD method) andfrom a macroscopic description of Navier-Stokes equations in the continuum hydrodynamics limit.The key point is that the macroscopic description is generally quite insensitive to the underlyingexact interactions among particles. Because of this reason, a very simplified and idealized micro-scopic model can be taken into account as long as the basic laws of physics such as conservationlaws, are satisfied. For this reason, an ?artificial grid-land? can be designed in a way to satisfy allfundamental conditions of flow equations like conservations of mass and momentum, rotationalinvariance and so on. Indeed, though the LBM known as a finite difference discretized version ofthe Boltzmann equation, historically it has evolved after realizing that fluid?s motion can be simu-lated with very simple discrete models which were developed based on physics laws (Lattice Gastechniques). Generally, in lattice gas simulation techniques, each node represents probability of thepresence of the group of particles in a velocity interval. This probability at each node is evolved bypropagation of molecules to their adjacent lattice sites and exchanging their momenta in the nextcollision.What has been explained tries to elucidate that there are possibilities to switch from micro-scaleto meso-scale simulations since computational cost of micro-scale simulations are much higher.Meso-scale simulations (specifically LBM) also have some advantages compared to the regularcomputational fluid dynamics techniques even in the range of continuum fluid dynamics. Firstly,LBM is known as a mesh-less technique and so flow problems especially multiphase ones in irreg-ular geometries can be handled with ease. Broadly speaking, easy handling of irregular boundaryconditions and mesoscopic forces which drive phase transition, simpler settling the hydrodynam-ics equations in the territory of soft matter research works when the hydrodynamic equations arenot well settled and easier coping with other related complexities hard to describe within a con-tinuum approach, are the advantages of using Lattice Boltzmann method instead of regular CFDtechniques. Moreover, lattice Boltzmann method allows solving the problem by massively parallelcomputing algorithms.The main objective of this section is to figure out the flow behaviour in uniform and non-uniform (like foam-paper internal structure) porous structures. In order to reach this goal, it isnecessary to import the foam-paper geometry to computer and compute the fibre distribution inthe thickness of a foam-paper. For this purpose, micro X-ray tomography technique is employed58to take high resolution images from cross section of some samples. The result of the X-ray to-mography device is thousands of 2D images containing salt and pepper noise. Image processingtechniques for the purpose of noise detection is applied to 2D tomographic images to make the datausable for lattice Boltzmann simulation. All the data are binarized (converted to lattice) using areasonable threshold. After preparing the 2D data, a 3D image reconstruction code is developed tomake the 3D geometry.Since the size of the 3D matrix is very large, it is decided to make a medium using 3D random cylin-der geometry. The cylinders are sampled out of the fibre distribution function which is computedby X-ray tomography technique. Though this is not the exact alternative of the actual geometry,it is very useful to understand the difference between the flow pattern in a non-uniform geometry(such as foam-paper geometry) and a uniform one. Very useful information can be extracted aboutthe difference between the permeability values in above mentioned medium.4.1 Kinetic theory and lattice gas cellular automata4.1.1 Kinetic theoryIn theoretical physics, intricate dynamics of non-equilibrium systems and their relaxation towardthermodynamic equilibrium was always a matter of question and interest. Kinetic theory is a branchof statistical physics dealing with this topic. The quantitative investigation on the dynamics of theflowing substance can be elaborated in two major levels of description: continuum level (macro-scopic) and molecular level (microscopic). Flows in a continuum framework are described in termsof field parameters such as pressure, density and velocity which are both space and time dependant.For most applications, the current view point and description are quite sufficient since the physicsof the flow is totally identified. Nonetheless, in some special cases this approach is not able todefine the whole picture. For these cases, micro-scale point of view becomes important. In thisapproach, the motion of individual molecules of the fluid is studied. The final flow will be the re-sult of the collection of molecular motion and interaction. Though this approach directly simulatesthe flow, it is computationally very expensive. To overcome this problem, statistical mechanicsproposed a third intermediate level of description, in which fluid is represented by a distributionfunction which acts as a representative for collection of particles. In other words, flows are repre-sented in terms of the probability f (t;~r,~v) of the presence of given fictitious ensemble of particlesat a given position in space,~r and time t and velocity~v. The current view point is called mesoscopicapproach and is described by kinetic theory.59In 1872 the Austrian physicist, Ludwig Boltzmann derived his celebrated equation for de-scribing the evolution of distribution function ( f (t;~r,~v)) in terms of micro-dynamics interactions.Although this equation is not an easy approach to a fluid dynamics problem mathematically andcomputationally, it contains more information due to this fact that f exists in a six-dimensionalphase space. The mathematical representation of the Boltzmann Equation (BE) is for the probabil-ity distribution function (PDF) f (t;~r,~v) and is as following[78, 79]:[?t +~v.?~r +~F .?~v] f (t;~r,~v) = Q12 (4.3)The left hand side of Equation4.3 describes streaming motion of the molecules with the pres-ence of the external force-field (~F) and the left hand side of the equation explains the effect ofinter-molecular collisions which diverts fluid?s molecules from streaming trajectory. Equation4.3shows that the collision operator for the evolution of one-particle PDF ( f ) depends on two-particlesPDF ( f12) which is concealed inside Q12. It is proved that the evolution of f12 depends on thethree-particles PDF ( f123) and generally, dynamics and evolution of the s-particles PDF dependson (s+1)-particles PDF. This chain of equations which imply the dependency of the evolution off1?s to f1?(s+1) is called BBGKY3 hierarchy. BBGKY-hierarchy is first derived from Liouvilleequation for the evolution of an N-particle system[80]. In order to terminate this hierarchy and toclose up the mathematical description of the collision operator in Equation4.3, some simplifyingassumptions are taken into account[80]. The most important assumption is that the fluid is consid-ered as a rarefied system so that the assumption of binary collision is valid and the probability ofternary (and/or higher order) collisions are practically nil[80]. The final Boltzmann equation for ararefied fluid regime is written as following:[?t +~v.?~r +~F .?~v] f =?( f1? f2?? f1 f2) ~vrel?(| ~vrel|,?s)d~v2 (4.4)As shown in Equation4.4 , the mathematical description of collision operator is very complicatedeven for a simplified binary collision in rarefied regimes.It should be noted that if all gains and losses in the system are totally balanced the collision termwill be annihilated. Explicitly, this means that both rate of change and gradient of change of theprobability distribution function is zero and all ensembles are constant. However, constant en-sembles does not mean that there is no intermolecular interaction between molecules. This means3Bogoliubov-Born-Green-Kirkwood-Yvon60that every direct intermolecular collisions is dynamically balanced by an inverse collision and viceversa. This condition is called local equilibrium. The probability distribution function for the lo-cal equilibrium is first proposed by Maxwell and Boltzmann and is known as Maxwell-Boltzmannequilibrium distribution. Equation4.5 demonstrates the celebrated Maxwell-Boltzmann equilib-rium distribution. Since Q( f e, f e) = 0 , the Maxwell-Boltzmann distribution for local equilibriumcan shown to be:f e = ?(2piv2T )?32 exp(??3?i=1(ci?ui)22v2T) and i = 1? 3 (4.5)where vT =?KBTm is the thermal speed which is related to the fluid temperature. Another profoundcontribution of Boltzmann was quantitative discovery of irreversibility. He defined a functionalcalled H-function as represented in Equation 4.6:H(t) =??f ln f d~vd~r (4.6)and proved that this function is a monotonically increasing function of the time i.e. dHdt ? 0 . Theequality sign happens when the system is in global equilibrium and the entropy of the systemis maximum.4 H-theorem is similar to the second law of thermodynamics. By considering thegeneralized Gibbs entropy into account, it can be easily found that H = SV KB which shows that fora system of interest with a fixed volume V , entropy never diminishes.4.1.2 A brief review of lattice gas cellular automata (LGCA)Before talking about the idea of lattice gas and application of Cellular Automata for solving thehydrodynamics equations, it will be appropriate to discuss briefly its history and to understand theidea behind it. CA5 is a discrete novel model and a creative approach for solving physical phe-nomena in nature and a new smart approach to solve differential equations. The method is studiedin vast range of sciences such as statistical physics, computational biology, social science, financeand economics. For the first time, in 1950, CA was developed and introduced by Stanislas Ulamand John von Neumann[81] at Los Alamos National Laboratories while Ulam was working on the4Global equilibrium is a condition that both mean flow speed and the flow temperature is constant. Global equilibriumtakes requires much longer time to happen comparing to the local equilibrium. In fact, being in the local equilibriumcondition, Q( f e, f e) = 0, does not put any restriction on the space dependence of temperature and flow speed.5Cellular Automata61growth of crystals, using a simple lattice network and Neumann was working on self-replicatingsystems. In 1969 a German mathematician and computer scientist, Konrad Zuse proposes his cel-ebrated theory, the Zuse?s Theory that the physical rules of the universe are naturally discrete andthe entire nature is the result of a deterministic computation on a cellular automaton[81]. Besidesapplying the method to hydrodynamics, Zuse applied this model for solving electrodynamics andquantum mechanics problems. Between 1970 and 1985, several scientists such as John HortonConway[82], Martin Gardner[83?85] and Stephen Wolfram[86] work on the idea of CA and it?sdifferent applications. Cellular Automata consists of regular arrangements of single grid of cells ofthe same kind so that each one of them has a finite number of different states. To start a simulationwith CA, a particular state is given to the cells as an initial state at t = 0. The new state of eachcell is updated synchronously in terms of the current state of each cell and the states of the localadjacent cells according to a determined and fixed set of rules at discrete time steps. The currentsystematic survey on the algorithm of the one dimensional CA is studied by Wolfram[86].Lattice Gas Cellular Automata (LGCA) is the first specific Cellular Automata technique devel-oped for the sake of fluid flow simulation in 1973 by Hardy, Pomeau and Pazzis[81]. This method isalso called HPP-LGCA which is derived from initials of its authors. HPP technique is able to con-serve mass and momentum. However, it is not suitable for the actual macroscopic flow simulationsince it is not able to satisfy the other characteristics of Navier-Stokes equations in macroscopiclevel such as rotational invariance[79, 81]. In 1986, Frisch, Hasslacher and Pomeau find that if aCA is run over a lattice with hexagonal symmetry, it leads to satisfy the hydrodynamics equationsin macroscopic level. This model is called FHP-LGCA[79, 81].Though, remarkable progress is achieved in order to make LGCA suitable for simulating flow equa-tions, the interest to use them for this purpose levelled off at early 90s. The major problem washigh amount of statistical noise of simulation of high Reynolds number flows. There are also otherproblems which can be found in the literature[79]. These problems are overcome by the idea ofLattice Boltzmann method.4.2 Lattice-Boltzmann method: hydrodynamicsIn 1988, the idea of Lattice Boltzmann technique to solve the discrete Boltzmann equation for flowsimulation purposes is first proposed by G. McNamara and G. Zanetti[79, 87] and then followedby Higuera and Jimenez[79, 88]. They try to minimize and/or escape the major problems in theLGCA approach which are statistical noise and exponential complexity of collision algorithm[79].62Though the main problems of LGCA are lifted by the previous proposed lattice Boltzmann models,the third noticeable limitation is still a matter of concern which is high momentum diffusivity. Thislimitation is due to the collision operator of the LGCA method which supports very low Reynoldsnumber flows and is also dominated by the works done by Higuera, Succi and Benzi[79, 88].One of the most serious problems of lattice Boltzmann techniques is raised because of the complexnature of the collision integral in the Boltzmann equation. Several authors propose different modelsto circumvent this issue. The simplest model which converges to the solution of the Navier-Stokesequations is the model proposed by Bhatnagar, Gross and Krook[79, 89, 90] in 1954 and calledBGK model. The idea behind their model is to simplify the collision operator which contains hugeamount of information about the binary interactions, into a much simpler operator which can pre-serve the essential physics of the problem. It is proved that most of the of details in the binarycollision integral of the Boltzmann equation does not significantly influence the final macroscopicresults of the many fluid dynamic problems in the moderate Reynolds numbers[79]. The BGKapproximation conserves the collision invariants of Q( f , f ) and also expresses the tendency toMaxwellian local equilibrium. These two together are enough to make sure that this model is accu-rate enough to model the collision integral for a wide flow regimes. The BGK model is representedas:Q( f , f ) =?f ? f e? (4.7)Parameter ? is called relaxation time an it defines the time required for the molecules to relax towardlocal Maxwellian equilibrium. The combination of lattice Boltzmann method with BGK model iscalled LBGK technique and it is used by many authors[79, 90]. In the lattice Boltzmann technique,the kinetic equation is solved for the probability distribution function, f , in the phase space and themacroscopic parameters such as density, macroscopic velocity and specific internal energy fieldscan be calculated by evaluating the hydrodynamic moment equations of the distribution functionsas following[79]:?(~r; t) = m?f (~r,~v; t)d~v (4.8)?(~r; t)~u(~r; t) = m?~v f (~r,~v; t)d~v (4.9)?(~r; t)e(~r; t) = m?(~v?~u)22f (~r,~v; t)d~v (4.10)63The kinetic model that is used in this survey is the Boltzmann equation with the so called singlerelaxation-time (SRT) approximation, the BGK model.[?t +~v.?~r] f =?1? ( f ? fe) (4.11)In order to solve the Equation 4.11 numerically, first it should be discritized in the velocity spaceusing the finite numbers of velocity vectors which must obey conservative laws and also preserveother properties of the flow equations.? fk? t + ck.? fk?~r =?1? ( fk? fek ) (4.12)where ck = ?x/? t is the kth discrete velocity. The current equation must also be discretized in timeand space. The final discrete lattice BGK equation is:fk(r+ ck? t, t +? t)? fk(r, t) =?? t? ( fk(r, t)? fek (r, t)) (4.13)The chosen discrete velocity space depicts the type of the lattice model which is selected. A 2Dnine velocity square lattice (D2Q9) and 3D nineteen velocity lattice (D3Q19) models are used forsimulating two dimensional and three dimensional flows respectively. The discrete velocity set forD2Q9 lattice model is as following:~ck =[0 1 0 ?1 0 1 ?1 ?1 10 0 1 0 ?1 1 1 ?1 ?1]The discrete velocity set for D3Q19 lattice model is as following[91]:~ck =???0 1 ?1 0 0 0 0 1 1 1 1 ?1 ?1 ?1 ?1 0 0 0 00 0 0 1 ?1 0 0 1 ?1 0 0 1 ?1 0 0 1 1 ?1 ?10 0 0 0 0 1 ?1 0 0 1 ?1 0 0 1 ?1 1 ?1 1 ?1???The equilibrium distribution function for both D2Q9 and D3Q19 models are:f eqi = wk?[1+~ck.~uc2s+(~ck.~u)22c4s?~u.~u2c2s] (4.14)where cs = 1?3 is the lattice speed of sound for both lattices, and pressure can be derived from the64ideal gas equation of state which is P = ?c2s . The values of wk are lattice weighting functions. Thelattice weighting functions for D2Q9 is:wk =?????4/9 k = 01/9 k = 1,3,5,71/36 k = 2,4,6,8(4.15)and for D3Q19 lattice is:wk =?????1/3 k = 01/18 k = 1? 61/36 k = 7? 18(4.16)In order to derive the incompressible flow equations (Ma =|~u | /cs 1) (the so-called Navier-Stokes equations) from discrete Boltzmann equations, Chapman-Enskog expansion[92] is applied.The corresponding equation for the dynamic viscosity in the Navier-Stokes equations can be de-rived as following[79]:? = (??0.5)c2s? t (4.17)Since the value for dynamic viscosity should be positive to be physically meaningful, the condition? > 12 also must be satisfied.The lattice BGK equation, Equation 4.13 can be computationally solved in a specific manner. If thephysical time step for computation (? t) is sufficiently less than the mean collision time of particles,the streaming and inter-molecular collision parts of Equation 4.13 can be decomposed into twoindependent local equations as following:f?k(~ri, t) = fk(~ri, t)?1? ( fk(~ri, t)? fek (~ri, t)) Collision (4.18)fk(~ri +~ck? t, t +? t) = f?k(~ri, t) Streaming (4.19)where f?k depicts the post collision probability distribution function.4.2.1 LBM versus Navier-Stokes solversOne of the most important questions which is required to be answered is what are the major differ-ences of a Lattice-Boltzmann based solver and Navier-Stokes solvers based on the regular contin-uum CFD techniques. In other words, it should be answered that who really needs LBM to solve a65fluid dynamics problem. In this section, some advantages and disadvantages of LBM are listed indetail.? One of the most challenging problems of regular NS6 incompressible solvers, is to solve aPoisson equation in order to compute the pressure field. This results in global data commu-nication and parallelization is relatively difficult. The condition for the LBM is reverse. Thecollision step in LBM is completely local and no calculation is happened in the streamingstep. The pressure field is also solved using an equation of state (for an ideal gas this equationis a linear function of density field, P = ?c2s ) [79]. This shows that the data communicationsin the whole LB algorithm is local making it easy to incorporate massive parallelization.? Uniform data shifting is another advantage of LB method compared to the conventional NS-solvers. Non-linear convective term in the Navier-Stokes equations (~u.?~u) is hard to treat.However, in the LB method, due to the linear operator in the Boltzmann equation(Equation4.4), the non-linear convective term is replaced by linear advective term (Equations 4.18 and4.19) and hence it is very straightforward to model.? The value for the lattice time-step and the lattice grid size are both equal to 1. So that theCFL7 number which is proportional to ? t?x in lattice Boltzmann method is 1. This implies thatthe rate of convergence in LB is totally imposed by the acoustic propagation which must bevery low in order to take care of the destructive effect of high lattice Mach number. Conse-quently, LB is considered as a low-speed convergent method for steady-state simulations.? Handling complex geometries and boundary conditions are matter of concern for both NS-solvers and LB method. In NS-solvers, the boundaries of the objects must be preciselydetected and the domain should be discretized with special and careful treatment. Howeverit is not always possible to mesh the domain properly due to the gross disorderedness andcomplexity (such as random disordered porous structures like the geometry of paper). More-over, proper handling of normal and tangential shear stress components on the boundaries isrequired. LB methods also suffer from problems related to curved boundaries and lack ofa specific counterpart in micro and meso-scale for some boundary conditions (for exampleno-slip boundary condition does not have a very accurate representation in micro-scale).6Navier-Stokes7Courant - Friedrichs - Lewy664.2.2 Applications of LBMNowadays the vast applications of Lattice-Boltzmann method in solving fluid dynamics problemshas become more common because of the availability of huge computational resources, espe-cially in the domain of parallel computing architectures. Lattice-Boltzmann incorporates complexphysics in a very simple and natural way by playing the details of particle-particle interactions orlevel of expansion of the distribution function in velocity space. In other words, complex physicscan easily be added to the fundamental equations of the Lattice-Boltzmann method. What ismeant by complex physics are two-phase flow problems(liquid-gas or solid-liquid transitions)[93],physics of the single-component or multi-component flows through random disordered geometriesand porous media[94], particle suspensions in fluid (such as sub-micron aerosol filtration)[95],chemically reacting flows and combustions[96], MHD8 flow[97], fully turbulent non-reactive andreactive flows[94, 98], crystallization and so on.For simulation of fluid flow through porous media, a lot of reasonable work has been done bothpurely on understanding the physics of the flow and the applications such as aerosol filtration, oilindustry and so on. The works which are done in FP-Innovation by Drolet[95] on the numericalsimulation of aerosol filtration in fibrous media are good examples. Rebai et. al.[95] carried out acomplete analysis on the effect of different computer-generated fibrous filters. They calculate theeffect of filtration efficiency of fibrous filters by changing different parameters such as aspect-ratioof fibres, crowding number and the Reynolds number using Lattice-Boltzmann method. Thomaset. al.[99] study the Stokes flow through uni-modal fibrous porous media using LBM. They usedcoarse numerical lattice technique to predict the permeability of the media in different solid frac-tions and different fibre diameters. Przekop and Gradon[100] simulate fluid flow and unsteadyparticle deposition in nano-structured fibrous media using the combination of LBM and Browniandynamics. They finally investigate the efficiency and the dynamics of deposition on the nano-fibrous media. There are several other papers similar to the presented works investigate differentparameters of fibrous filters and their effect on the filter efficiency.4.2.3 Source of errors in LBMA matter of concern of every numerical method to solve hydrodynamics equations is about un-avoidable errors. In this section, the sources of errors for the lattice Boltzmann technique areintroduced.8Magnetohydrodynamics67? Domain discretization errors: One of the main error sources in numerical methods is dueto the problem?s domain discretization into finite number of cells or nodes. The level foraccuracy needed in LB simulations determine the largest lattice size that can be used. It isvery common to determine the size of the lattice spacing based on the minimum scale lengthof the geometry. For instance the lattice spacing sizes in the flow simulation inside a 2Dchannel and around a 2D cylindrical obstacle, are the width of the channel and the radiusof the cylinder respectively. In the flow simulation in porous media it is very common touse the minimum pore size as the scale length which sometimes results in very fine latticenodes spacing near percolation threshold i.e. many pores have very small size leading tohuge simulation time. Several works are done to circumvent this problem by applying localgrid refinement methods[101, 102].It has been shown that in LBGK method, smaller values of the BGK relaxation time, de-creases the finite size effect. On the other hand very small values of relaxation time (? closeto 0.5) might result in some errors and numerical instabilities[79, 102].? Compressibility errors: Compressibility errors are caused and propagated due to the depen-dency of the pressure field to density in Lattice-Boltzmann simulation. This is in contra-diction with the assumption of incompressible fluid. In order to simulate the incompress-ible flow equations with LB method which is inherently quasi-compressible, the simulationmust be done in low Ma number. The error in the LB due to the compressibility effect isO(Ma2)[79]. The value of lattice Ma must not locally be greater than 0.3 [79].4.2.4 Boundary conditionsOne of the main strong points of lattice-Boltzmann method is its capability of handling different(even very complex) boundary conditions, especially handling boundary conditions according togrossly irregular geometries which can be illustrated by the case of no-slip boundary conditions atsolid walls. Two different approaches are applied to model boundary conditions for LBM. One iscalled ?Distribution Modification? and the other is ?Distribution Reconstruction?[103, 104]. Thedistribution modification approach tries to imply that some physical rules such as conservation lawsand bounce-back rules are able to obtain the unknown distribution functions. The distribution re-construction approach tries to recreate the distribution function on the boundaries by relating themto the macroscopic properties and their derivatives. The half-way bounce-back[79, 102], full-waybounce-back[79] boundary conditions for wall no-slip condition in solid-fluid interactions and theso-called Zou and He boundary condition[105] and Inamuro boundary condition[106] for mod-68elling of the open boundaries are examples of the distribution modification approach. Methodssuch as finite-difference boundary condition[107], regularized boundary condition[108] and theone proposed by Guo et al.[109] for modeling of the open boundaries are categorized in the distri-bution reconstruction approach.It should be noted that all the mentioned boundary conditions are developed for the boundariesand geometries which can be described by horizontal and vertical lines in Cartesian coordinates.Special treatment based on the interpolation and/or extrapolation techniques must be taken intoaccount in order to calculate much more precise approximation of those values that must be com-puted on the nodes away from the walls of the solid geometry in the Cartesian/lattice coordinatesystem. The interpolation and/or extrapolation techniques must be utilized since the boundarynodes can not locate on the curved boundaries (lattices configuration) as it can for methods likefinite difference in regular CFD. The method developed by Filippova and Hanel [110] (FH bound-ary condition), the modified method developed by Mei et al.[111] in continuation of the work doneby Filippova and Hanel (FHM boundary condition), the method proposed by Bouzidi et al.[112](BFL boundary condition) and the method proposed by Yu et al.[102] (Yu boundary condition) arethe examples of distribution modification approach for the curved wall boundaries. However, themethod proposed by Verschaeve et al.[108] which is based on the regularized boundary conditionsfor the open curved boundaries is a good example of distribution reconstruction.4.2.4.1 No-slip boundary conditionIn fluid dynamics the no-slip boundary condition is used in the solid-fluid interaction and states thatat a solid boundary, the velocity of the fluid is zero relative to the boundary interface. In lattice-Boltzmann method, this can be achieved by implementing an algorithm based on particle?s bounce-back on the boundaries[79]. Though the algorithm is simple, it can preserve the particles? momenta.Other techniques for simulating no-slip boundary condition is based on the idea of bounce-backalgorithm but they provide more detail in order to diminish the errors on the boundaries.4.2.4.1.1 Bounce-back algorithmDuring the streaming process, the components of the distribution function which are propagatinginto the solid nodes9 are bounced-back exactly in the opposite direction to its adjacent fluid-node.This means that the portion of the particle propagated into the solid boundary is reflected back to9solid-nodes are those nodes that are inside the solid geometry and are in contact with at least one node out of thesolid geometry (fluid-node)69the adjacent fluid-node in the opposite direction. In this process, the momentum of the particleis reversed. The applied surface force and applied torque are also calculated from the momentumtransfer at each boundary node and summed over the whole number of solid-nodes. The full-waybounce-back method can be described as following:f?k?(~rb, t) = fk(~rb, t) (4.20)where fk(~rb, t) is the set of the incoming populations from computational domain into solid bound-ary before collision and f?k?(~rb, t) is the set of the outgoing populations after collision from solidboundary to the computational domain. The full-way bounce back algorithm is a first-order accu-rate algorithm for no-slip boundary condition.In order to increase the accuracy of algorithm to second-order algorithm, the idea of half-waybounce-back can be taken into account. In this method, the solid boundary locates halfway thelink (between a solid node and a fluid node). In this algorithm solid nodes are practically act asghost nodes so particles bounce back before reaching the actual solid node in halfway the link.From computational point of view, in the half-way bounce-back algorithm, at a given time t thepost-collision values which are given by the population in the opposite direction at time t?1, areassigned. The collision step takes place after this step. Computationally, this is totally differentfrom the full-way bounce-back method since the full-way algorithm takes place after the collisionstep. This implies that the half-way bounce-back algorithm keeps the normal collision step butmodifies the streaming step quite differently in comparison to the full-way algorithm which mod-ifies the collision step and keeps the normal streaming step. The mathematical description of thehalf-way bounce-back algorithm is as following:fk?(~r f , t +? t) = f?k(~r f , t) (4.21)4.2.4.1.2 Filippova-Hanel-Mei (FHM) algorithmThe current algorithm is developed to take care of the solid-fluid interactions on the curved bound-aries by applying a interpolation/extrapolation technique. Two major problems are involved inmodeling of the curved boundaries in lattice Boltzmann method. The first problem is that thesolid boundary nodes (wall nodes) are not exactly located on the lattice coordinates (or even atthe halfway the link between sold nodes and fluid nodes). The second problem is that the fieldmacroscopic parameters cannot be computed by only the conservation laws since the number ofunknown populations are more (or sometimes less) than the equations. In this section, the algo-70rithm introduced by Filippova and Hanel and the modification that is done by Mei et al.[110, 111]is elaborated. As it can be seen in Figure 4.1, the vertical distance between the wall node and itsadjacent fluid node is ?.?x where ?= |~r f?~rw||~r f?~rb| and obviously 0? ?? 1.ckck~?.?x?xfluid nodessolid nodeswall nodesFigure 4.1: Curved boundary condition and lattice orientation.In this method, after the collision step, the distribution function for the solid nodes whichare linked at least with a fluid node should be updated based on an interpolation technique. Thesolid-to-fluid post collision distribution function can be interpolated by the values of post collisiondistribution function which is computed in the collision step (referring to Equation 4.18) and afictitious equilibrium distribution function ( f ?k )[110, 111]:f?k?(t,~rb) = (1??) f?k(t,~r f )+? f?k (t,~rb)+2wk?(t,~r f )~u(t,~rw).~ckc2s(4.22)where~u(t,~rw) is the velocity of the moving wall. The last term will vanish if the walls are stationary.The fictitious equilibrium distribution function is calculated as:f ?k (t,~rb) = wk?(t,~r f )(1+~ck.~u(t,~rb)c2s+(~ck.~u(t,~r f ))22c4s+~u(t,~r f ).~u(t,~r f )2c2s) (4.23)71where~u(t,~rb) is calculated as following:{~u(t,~rb) = ~u(t,~r f ) ?< 12~u(t,~rb) = (??1? )~u(t,~r f )+(1?)~u(t,~rw) ??12(4.24)In Filippova and Hanel algorithm the interpolation parameter, ? is calculated as:{? = 2??1??1 ?<12? = 2??1? ??12(4.25)The method proposed by Filippova and Hanel is improved by Mei et al.[110, 111] since the stabilityof their method is weak as ? ? 1. The following formula for ? is proposed by Mei et al.:{? = 2??1??2 ?<12? = 2??1? ??12(4.26)4.2.4.2 Open boundary conditionsBoundary conditions such as inlet/outlet boundaries, boundaries with free-slip condition (lines andplanes of symmetry) and periodic boundaries are categorized in the scope of the open boundaries.In this section, algorithms for open boundaries in lattice Boltzmann method are briefly discussed.4.2.4.2.1 Inlet/Outlet boundary conditionSeveral methods are proposed to model the inlet/outlet conditions for lattice Boltzmann method.One of the most successful and straight-forward methods to implement the Dirichlet boundaryconditions at inlet/outlet is the method proposed by Zou and He[105] based on the general rules ofconservation of mass and momentum at boundaries. The idea is derived by applying the bounce-back rule for the non-equilibrium part of the distribution function in the normal direction wherethe other two unknown functions are derived by the conservation rules. Zou and He boundarydynamics can be considered as extended bounce-back rule for open-boundaries. The probabilitydistribution function can be decomposed into the equilibrium part and the non-equilibrium part asrepresented in Equation 4.27 .fk = fek + fnek k = 0? 8 (4.27)7226 537 4 81Wall nodesInlet/OutletFluid nodesKnown PDFUnknown PDFFigure 4.2: Unknown and known probability distribution functions in D2Q9 lattice.For instance the x-component of velocity at inlet (west boundary) is uinx and the y-componentof velocity is uiny as it is shown in Figure 4.2 for D2Q9 lattice. Recalling from Equations 4.8, 4.9and ?? macroscopic field variables can be re-written as:? in =8?k=0fk (4.28)? inuinx = ( f1 + f5 + f8)? ( f6 + f3 + f7) (4.29)? inuiny = ( f2 + f5 + f6)? ( f4 + f7 + f8) (4.30)The equilibrium condition normal to the boundary (Equation 4.27) results in:f1? fe1 = f3? fe3 (4.31)where the equilibrium functions in Equation 4.31 can be calculated by Equation 4.14 . SolvingEquations 4.31 , 4.28 , 4.29 and 4.30 the four unknowns can be found as[105]:f1 = f3 +23? inuinx (4.32)f5 = f7?12( f2? f4)+16? inuinx +12? inuiny (4.33)f8 = f6 +12( f2? f4)+16? inuinx ?12? inuiny (4.34)73? in = 11?uinx( f0 + f2 + f4 +2( f3 + f6 + f7)) (4.35)The equations for other boundaries (North, East and South) and other lattices (such as D3Q19 andso on) can be derived using the same technique[105].4.2.4.2.2 Symmetry boundary conditionIn many practical problems some geometrical symmetry can be identified. Hence, in order to re-duce the computational time, the solution can be found only in one part of the problem and assumethat the identical solution is repeated beyond the lines (and/or surfaces) of symmetry. Mathemat-ically, symmetric boundary is usually understood by setting the component of velocity normal tothe boundary to be zero:~u(~r, t).n? = 0 (4.36)In Lattice-Boltzmann method this condition can be satisfied by setting the unknown distributionfunctions equal to their mirror images[79] . The mathematical representation of such a condition,for example, for the south boundary( f2, f5and f6 are unknowns) of Figure4.2 is as following:?????f2 = f4f5 = f8f6 = f7(4.37)4.2.5 Numerical stabilityOne of the main issues of different algorithms in fluid dynamics are problems related to the sta-bility of that algorithm. Lattice Boltzmann technique is also not immune to the error dispersionsrelated to the stability problems. Numerical instability in Single Relaxation LBM (SRT-method)are mainly arises due to the poor trade off between the relaxation parameter (?) and the latticeMach number (Ma). The problem arises because of inadequate-definition of the initial conditionwhich results in high values of the lattice Mach number and spurious pressure waves. A very com-prehensive research is done by Sterling and Chen[113] at Los-Alamos National Laboratories on theeffect of different parameters on the SRT-method?s stability. The distribution of the mass at a sitebetween the different discrete speeds, the relaxation time in BGK model, the mean velocity and thewave number of perturbations are introduced as parameters on which the stability of SRT-methodare depends on them[113]. Another important reason for numerical instabilities in SRT Lattice-74Boltzmann method is singular points in geometries. For instance, in a lid driven cavity problem,the two upper corners are singular points in the geometry which may result in a very noisy pressurefield [113]. This effect is reduced by applying multi-block method[102] but, in the case of transientflow at high Reynolds number it is very hard (sometimes impossible by increasing the Reynoldsnumber) to damp the acoustic wave propagated due to the singularity in geometry[102]. This ismainly because in SRT lattice Boltzmann, both bulk viscosity and shear viscosity are relaxed bya single value[102]. It should also be noted that the stability of boundary conditions (specificallyopen boundaries) is very important. The stability of Zou and He boundary condition is guaranteedin the range of low and moderate Reynolds numbers[102, 105]. However this boundary conditionbecomes vulnerable (and sometimes totally unstable) in high Reynolds numbers. In the case of no-slip boundary conditions in solid-fluid interaction, Yu et al.[102] compare the stability of differentboundary conditions in different flow regimes.In order to circumvent these problems, different algorithms are proposes which are more stablethan the SRT-method. D?Humieres[114] proposes an algorithm based on the general collisionmatrix for collision operator that is called Multi Relaxation Technique (MRT). The idea of MRT-method is to relax different modes with different relaxation times while all modes in the SRT-method are relaxed with a single value of relaxation. A lot of scientific investigations is done on theMRT Lattice-Boltzmann method. Mei et al. compare the performance of SRT and MRT Lattice-Boltzmann methods in detail. A comprehensive study on stability characteristics and dispersion onMRT Lattice-Boltzmann method is carried out by Lallemand and Luo [115].4.3 Results and discussionsThis section is divided into two parts. In the first part, some 2D and 3D benchmark problems aresolved for code validation purpose. Error and stability analysis of the code are also presented in thissection. The results and discussions about the flow simulation inside 3D model geometry of foam-papers will be presented in the second part of this chapter. Codes for all 2D and 3D simulations aredeveloped in C++ language in Linux platform from scratch. In both 2D and 3D codes the memoryis dynamically allocated for handling huge data structure. Moreover, both 2D and 3D codes areparallelized for multi-core machines using OpenMP computer architecture. Complementary codesfor geometry generation (for example generating a geometry for 3D random fibre medium) is alsowritten from scratch for the simulations.754.3.1 Benchmark studiesIn this section, for the purpose of code validation, flow in a 2D channel at low Re number, flowaround a circular cylinder at low and moderate Re numbers and pressure driven flow around bundleof 2D cylinders are solved. For the validation of 3D code, pressure driven flow inside a uniformrandom orientated cylinders (a model geometry for uniform fibrous porous media) are solved.4.3.1.1 Flow inside a 2D channelThe first benchmark study is done on a 2D flow inside a channel geometry at a low Reynolds num-ber. Zou and He boundary condition is applied for the inlet and outlet open boundaries. A uniformvelocity profile is considered as an inlet boundary condition and a constant pressure is applied tothe outlet boundary. No-slip boundary condition is applied to the top and bottom wall boundariesusing full-way bounce-back algorithm. The channel aspect ratio is 25, the grid resolution in thewidth of channel is 80 and the flow solved for Re = 4.0.0 0.2 0.4 0.6 0.8 100.511.52y/HU/U inLBM (? = 0.6)Analytical(a)0 0.5 1 1.5 2x 105?12?10?8?6?4?2Iteration numberResidualsRes?URes?VRes??? = 0.6(b)Figure 4.3: (a) Lattice-Boltzmann simulation of a 2D channel flow at Re = 4.0. (b) Residuals.Figure4.3a shows the results of x-velocity in the thickness of the channel at fully-developedregion. Comparing the LBM results to the analytical solution of 2D channel flow shows an excel-lent agreement. The residual plots are shown in Figure4.3b. The residuals are defined based on themaximum difference of a macroscopic parameter between the current and the previous time-step.In order to be more precise about the simulation parameters and source of errors in LB simulations,domain discretization and compressibility errors are studied. Figure4.4 shows how accuracy ofthe solution changes by increasing the grid resolution in the domain. Beside having a error in by76applying first-order accurate full-way bounce-back method on the channel walls, the other sourceof error causes using Zou and He algorithm in the simulation of open boundaries as shown in Fig-ure4.5. This algorithm becomes unstable at high Reynolds number flow regimes due to increasingthe error overshoot on boundaries by increasing Re number.20 40 60 8010?410?310?2NyL2 Norm Error? = 0.6Figure 4.4: L2 norm errors of x-velocity for five different grid resolutions.3300 3400 3500 3600 3700?6.5?6?5.5?5Number of time?steps2000 4000 6000 8000 10000?7?6.5?6?5.5?5?4.5?4Number of time?stepsResidualsRes?URes?VRes??Figure 4.5: Errors made by applying Zou-He boundary conditionsFigure4.6a shows how the solution of the channel flow changes by increasing the lattice Manumber at the same number of iterations, same grid resolution and same Re number. Figure4.6bshows the effect of Ma number on the RMS error of the solution. By increasing the Ma numberin the simulation domain the solution loses its accuracy. By increasing the value of Ma number toMa = 1.0 the simulation becomes unstable and divergence occurs as is shown in Figure4.7.770 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 100.20.40.60.811.21.41.61.82y/Hin Ny =u40,uReu=u4, ? =u1.70,uMau=u0.100Ny =u40,uReu=u4, ? =u1.00,uMau=u0.045Ny =u40,uReu=u4, ? =u0.80,uMau=u0.026Ny =u40,uReu=u4, ? =u0.70,uMau=u0.014Ny =u40,uReu=u4, ? =u0.60,uMau=u0.006AnalyticaluSolutionSimulations at 20000 iterations.Channel aspect ratio (L/H): 25.0U/U(a)10?2 10?110?810?610?4MaRMS Error(b)Figure 4.6: Compressibility errors (a) Solution (b) RMS error vs. lattice Ma number.0 20 40 60 80?2024681012Number of time stepsResidualsRes?URes?VRes??(a)980 1000 1020 1040 1060?2024681012Number of time stepsResidualsRes?URes?VRes??(b)Figure 4.7: Compressibility instability plots at Ma = 1.4.3.1.2 Flow around a 2D circular cylinderThe second benchmark study is done on a 2D flow around a circular cylinder geometry at low andmoderate Reynolds numbers. Zou and He boundary condition is applied for the inlet and outletopen boundaries. A uniform velocity profile is considered as an inlet boundary condition and aconstant pressure is applied to the outlet boundary. No-slip boundary condition is applied to thetop and bottom wall boundaries using full-way bounce-back algorithm. The channel aspect ratiois 6.25, the grid resolution in the width of channel is 160, the blockage ratio (Ratio of the channel78width to cylinder diameter) is 8 and the flow solved for Red = 1,5,10,20 and 40. Second order ac-curate FHM boundary condition is implemented to simulate the no-slip condition on the surface ofthe cylinder more accurately. Figure4.8 a shows how those solid-nodes which are in contact withat least one fluid-node are identified by the code. The values of geometric parameter, ?, is alsoplotted in Figure4.8b. This figure shows that the value of ? is between 0 to 1 as it is predicted inFHM theory. Figure4.9 shows the results of the Drag coefficient in different Red numbers for bothLB method and CFD. The CFD simulation is done using commercial software ANSYS-FLUENT14.5. The results show excellent agreement compared to the CFD results. The centre of the cylin-der is located at xc = Nx/4 and yc = Ny/2.240 245 250 255 260707274767880828486889020 40 60 80 100 120 140 160 18000.20.40.60.81Number of solid nodes?(a)(b)Figure 4.8: (a) Solid-nodes which have a connection at least to one fluid-node (b) Value ofgeometric parameter ?.790 10 20 30 4005101520RedC DFV methodLBMFigure 4.9: Drag coefficient vs. Red number.In order to be more precise about the number of grid resolution which is used to simulate theflow around the cylinder, a sensitivity test study is done. Figure4.10 shows the minimum numberof nodes which should be considered for a cylinder in order to have accurate results. Minimumnumber of 8 nodes per radius of cylinder is required to have accurate results. This number willchange by changing the Red number (Accurate resolution for Re = 1 is r = 4).2 4 6 8 10 12 142.662.712.762.8CylinderVradiusV(latticeVunits)C DLBMFVVmethodReV=V10Figure 4.10: Drag coefficient vs. cylinder resolution at Red = 10.80Figure4.11 shows the contour plot of velocity magnitude around the circular cylinder in Red =1,5,20 and 40. The result depicts how the flow pattern changes by increasing Red number asexpected.(a)(c)(b)(d)Figure 4.11: Velocity magnitude around circular cylinder at low and moderate Red .4.3.1.3 Pressure driven flow around bundle of 2D cylindersThe third benchmark study is a pressure driven flow around bundle of 2D cylinders in low Reynoldsnumber regime. A periodic arrays of identical cylinders which are oriented perpendicular to theflow are considered. Since the geometry is periodic in this particular case, the solution of theproblem is reduced to one cylinder but the boundary conditions is changed to periodic boundaries.The porosity of the geometry can be changed by changing the diameter of the cylinder in the flowdomain. Zou and He boundary condition is applied for the inlet and outlet open boundaries. Inorder to apply the open boundaries, a buffer zone is taken into account before the inlet and afterthe outlet. This buffer zone must be large enough (in this simulation the size of the buffer zone isNx10 ) to minimize the interaction between inlet and outlet of the simulation domain while periodicboundary condition is applied. Second order accurate FHM boundary condition is implemented to81simulate the no-slip condition on the surface of the cylinder more accurately.0.7 0.75 0.8 0.85 0.9 0.9510?210?1100101102?K/d2 fSanganisLandLAcrivosLModelDrummondLandLTahirLModelLBMLSimulationL(RebaiLetLal.)LBMLSimulationFigure 4.12: Normalized permeability vs. porosityResults of normalized permeability vs. porosity are presented in Figure4.12. The results arecompared to the models presented by[76, 77] and the LBM simulations done by Rebaiet al.[95].The discrepancy from analytical models in lower porosities are expected since the models are validfor high porosity media[75].4.3.1.4 Flow inside a uniform 3D random fibrous structureThe fourth benchmark study is for validation of the 3D code. A uniform distribution of randomentangled cylinders are generated to simulate the geometry of fibrous porous media and filter-like structures in different porosities. The geometry made in a way that non of the cylinders canintersect with each other in the simulation domain. Resolution of r = 4 is used since the flow issolved in Red < 2. The aspect ratio of all the cylinders are the same and equal to 5 (L/d = 5). Forsake of diminishing the huge amount of CPU time and large amount of memory, periodic boundarycondition is applied to the side walls of the channel. Zou and He boundary condition is applied forthe inlet and outlet open boundaries. In order to apply the open boundaries, a buffer zone is takeninto account before the inlet and after the outlet. This buffer zone must be large enough (in thissimulation the size of the buffer zone is Nx6 ) to minimize the interaction between inlet and outlet of82the simulation domain while periodic boundary condition is applied. Second order accurate FHMboundary condition is implemented to simulate the no-slip condition on the surface of each cylindermore accurately. The other details of the simulations is reported in Table4.1. The sample geometryis shown in Figure4.13.0 0.5 11.5 201200.511.52XYZFigure 4.13: Sample uniform 3D fibrous geometry.Table 4.1: 3D Benchmark test case.Porosity No. of cylinders Nx Ny Nz93.5 340 330 220 22095.5 234 330 220 22097.5 132 330 220 220Results of normalized permeability vs. porosity are presented in Figure4.14. The results arecompared to Davies empirical model[75] and the LBM simulations done by Rebai et al.[95]. Thesmall discrepancy from empirical model of Davie is expected du to the random nature of themedium and small numerical errors in the LB algorithm and numerical implementation of boundaryconditions.830.65 0.7 0.75 0.8 0.85 0.9 0.9510?210?1100101102?K/d2 fDaviesbModelLBMbSimulationb(Rebaibetbal.)LBMbSimulation0.92 0.93 0.94 0.95 0.96 0.97 0.98100?K/d2 fDaviesbModelLBMbSimulationb(Rebaibetbal.)LBMbSimulationFigure 4.14: Normalized permeability vs. porosity (3D case).Figure 4.15: X-velocity: 3D uniform random fibrous media.Figure4.15 shows the x-velocity in three different planes perpendicular to each other in orderto show how the pattern of flow changes inside a grossly disordered geometry.844.3.2 Simulation of flow inside model geometry of a foam-paperTo simplify the complexity of geometry and reduce the computational time of the simulations,a model based on random oriented cylinders are proposed to be used instead of simulating themesh came out of the tomographic results. Since the geometry of foam-paper is a non-uniformporous geometry, the data for fibre concentration in the thickness of foam-papers is required. Thisdata is provided by Micro X-ray tomography from the cross-sections of the samples as shown inFigure4.16. Figure4.16 actually shows the fibre concentration distribution in the thickness of afoam paper.Figure 4.16: Tomographic images of side and top view of a foam-paper.0 500 1000 1500 200005101520Foam?paper thickness (?m)Fibre Concentration (%)NBSK (21% AC)Figure 4.17: Fibre concentration frequency in the thickness of foam-paper.The result of Figure4.17 is used to generate random numbers from this distribution function85instead of using the uniform random number generator in the thickness direction of the model ge-ometry domain. Figure4.18 shows the random numbers which are generated from the distributionof fibre frequency in the thickness of a foam-paper (results of x-ray tomograph, Figure4.16). Thegeometry for the simulation is made from these random numbers in two different porosities asreported in Table4.2.Figure 4.18: Random numbers generated from the distribution of fibre frequency in the thick-ness of a foam-paper.Table 4.2: Foam-paper simulation cases.Porosity No. of cylinders Nx Ny Nz93.5 338 330 220 22095.5 231 330 220 220Figure 4.19a shows the residual plots for hydrodynamic parameters to illustrate the convergenceof the simulations. Since the simulations are done by applying a pressure drop between the inletand the outlet of the channel, the flowrate will evolve in the simulation. Since the simulations arein low Reynolds regime, the value for the average velocity at the inlet must converge to a constant86value. This condition is another criteria for the convergence of low Re number flows as shown inFigure4.19b. In order to guarantee the stability of simulations, the value of maximum Ma numberis presented at each time-step as shown in Figure4.19c (Final value of Mamax = 0.024). This valueis lower than the critical Ma number value (Mamax < 0.3) [79].2 4 6 8 10x 104?11?10?9?8?7?6?5?4Number of time?stepsResidualRes??Res?URes?VRes?W(a)0 2 4 6 8 10x 10405101520253035Number of time?stepsRe inlet(b)0 2 4 6 8 10x 10400.0050.010.0150.020.025Number of time?stepsMa max(c)Figure 4.19: Convergence criteria (a) Residuals (b) Reinlet (c) Mamax.(a)(b)Figure 4.20: Comparing the x-velocity simulation of a foam-paper with a uniform 3D randomfibrous media at same porosity.87Figure 4.20 compares the x-velocity simulation of a foam-paper with a uniform 3D randomfibrous media at the same porosity. This qualitative result show the difference between the fibresfrequency distribution of a foam paper and a uniform medium. The difference between the flowpattern in both geometries are also shown in Figure4.20.0.93 0.94 0.95 0.96 0.97100?K/d2 fDavies)ModelLBM)(uniform)medium)LBM)(foam?paper)model)(a)0.93 0.935 0.94 0.945 0.95 0.95510?1100101?K/d f2DaviescModelLBMc(uniformcmedium)LBMc(foam?papercmodel)NBSK (Experimental)(b)Figure 4.21: Convergence criteria (a) Residuals (b) Reinlet (c) Mamax.Results of normalized permeability vs. porosity are presented in Figure4.21a. The resultsare compared to Davies empirical model[75] and the LBM simulations of a uniform 3D randomfibrous media which is done by author. The results of foam-paper simulations are compared to theexperimental results of actual NBSK foam-paper as shown in Figure4.21b.Cylinders that are used here for the simulations are rigid bodies and have an aspect ratio of 5 whichis much different from the actual foam-paper fibre strands which are flexible and have an aspectratio range of 60 - 130. It may be said that since the cylinders are stiff as well as random in theirarrangement, there are void spaces in between which cannot be occupied. The situation in thereal case is much less stringent as the fibres can bend which, is assumed, also increases blockage inflow leading to high tortuosity and low permeability. Also, in the simulations, random generation offibres is carried out in a way that no two cylinders intersect or are in tangential contact. This is notthe case with foam-paper micro-structure as fibre strands are not completely spatially segregatedand there may be a large number of points where one single strand is in contact with many otherstrands. This fact suggests that the flow around a real fibre strand is different in the vicinity ofthese contacts. A cylinder in simulation on the other hand is completely surrounded by a flow field.Hence there may be some flow leakage in simulations as the fibre-fibre contacts are not taken intoaccount.Now if there are two different foam-paper blocks both with the same length, cross-section and88average porosity, but one is uniform and the other has a distribution, that is, the local porosityvaries along the length, the extent of difference in permeability values for the two different blockswill depend on the porosity distribution of the non-uniform block. If the minimum and maximumlocal porosities for the non-uniform block have a small difference, the resulting permeability willnot be significantly different from that of the uniform block. It is assumed that this will not continueto hold true if there is a huge overall variation in the distribution or steep gradients in porosity alongthe length.89Chapter 5Conclusions and final remarksSince the total amount of paper consumption in the north America was diminished approximately25% in the past five years, it is necessary to make new products with different applications fromcellulose fibres. Foam-forming is one the potential ideas to make new products with various appli-cations out of cellulose fibres. Foam-forming process can be simply implemented and attached toan automated paper-machine. The 3D, honeycomb structure of foam leads to making high porousmaterials out of cellulose fibres called foam-paper. The experimental tests show pretty reasonableresults in aerosol filtration of sub-micron particles, heat insulation and sound absorption.In the sub-micron air-filtration tests, foam papers made by NBSK and CTMP pulp can notbe used as efficient filters because of their poor filtration properties. However, the results showthat the key parameters for making efficient filters from cellulose fibres are crowding number andfibre specific surface. The experimental results show that having 10% to 30% weight ratio ofNBSK valley beaten fibres in the structure of foam-papers improves both their mechanical proper-ties and filtration characteristics compared to the commercial filters. Moreover, foam papers madeby nanofibrillated fibres show excellent filtration properties at high foam air-contents. The onlyproblem is that the process of making them (Freeze-drying) is not economically affordable. More-over, at very high foam air-contents they show very poor mechanical properties. For the next step,the author recommends to do a research on filtration properties of air-dried foam-papers by addingdifferent weight ratio of nanofibrillated fibres to their structure.Comparing the results of thermal conductivity of NBSK foam-papers with commercial heatinsulators show that the new products made in high foam air-content can be used as heat insula-90tors. Multilayer NBSK foam-papers also show reasonable acoustic properties. However, the bestresults are obtained from samples which are made by nanofibrillated fibres. Overall, based on theexperimental data in the range of the experiments of this dissertation, there is a possibility to thinkabout commercialization, however doing much more research is required.Understanding the physics of the flow inside the disordered geometries of foam-papers is an-other aim of this research. Practically, it is hard to use regular Navier-Stokes solvers for these kindof geometries because of so many issues which has already been discussed in this thesis. An Open-MP parallel C/C++ code for handling large data structures is developed for multi-core machines inLinux platform. Since lattice Boltzmann method is a local numerical algorithm, it can be massivelyparallelized for parallel machines. In the modeling of foam-papers using the model geometries ofcylinder, it should be noted that due to computational constraints, the number of fibres that is usedin the simulations is less than what it actually is in the foam-paper geometry. Nevertheless, thesimulations may be considered to be a preliminary step to numerically studying the flow throughfoam-paper media as they do provide ample amount of qualitative information about nature of flowthrough this kind of domains which have highly randomized geometries. More precise results maybe obtained by carrying out the simulations for a larger ensemble of fibres as that is expected to bea better statistical approximation of the porous medium. Also the random fibre generator can beimproved to mimic the real geometry. Aspect ratio of cylinders used in the simulations is kept at 5whereas in reality it is in the rage of 60 - 130. Also the real fibres are present in bent and twistedshapes and are in contact with other fibres. This geometrical details should be analysed and takeninto consideration for randomly generating the fibres.91Bibliography[1] Jim Ford; Susan Kinsella; Gerard Gleason; Andrew Goldberg; Robin Averbeck. The stateof the paper industry. 2011. ? pages ix, 1, 2[2] Timo Lappalainen and Jani Lehmonen. Determinations of bubble size distribution offoam-fibre mixture using circular hough transform. NORDIC Pulp and Paper Research, 27(5):930?939, 2012. ? pages 6, 10[3] D. Weairea and R. Phelana. A counter-example to kelvin?s conjecture on minimal surfaces.Philosophical Magazine Letters, 69(2):107?110, 1994. ? pages 7[4] Christopher James William Breward. The Mathematics of Foam. PhD thesis, TheUniversity of Oxford, 1999. ? pages 7[5] J.C.Isarin; A.D.J.Kaasjager and R.B.M.Holweg. Bubble size distribution during theapplication of foam to fabrics and its effects on product quality. Philosophical MagazineLetters, 65(2):61?69, 1985. ? pages 8[6] Samuel P. Gidoa; Douglas E. Hirta; Susan M. Montgomerya; Robert K. Prud?hommea andLudwig Rebenfelda. Journal of Dispersion Science and Technology, (6):785?793. ? pages8[7] A.B.J. Kroezen; J. Groot Wassink and C.A.C. Schipper. Journal of the Society of Dyersand Colourists, pages 393?400. ? pages 8[8] A.B.J. Kroezen and J. Groot Wassink. Journal of the Society of Dyers and Colourists,pages 397?401. ? pages 8[9] D. E. Hirt; R. K. Prud?homme and Ludwig Rebenfeld. AIChE Journal, (6):326?328. ?pages 8[10] R. C. Chang; H. M. Schoen and C. S. Grove. Industrial and Engineering Chemistry, (11):4a?6a. ? pages 8[11] K. A. Brakke. Experimental mathematics 1, pages 141?165. ? pages 892[12] Yi Jiang. Extended large- q potts model simulation of foam drainage. PhilosophicalMagazine Letters, 74(2):119?128, 1996. ? pages 8[13] J.J.Bikerman. Foams. Springer-Verlag, 1973. ? pages 9[14] John P. Heller and Murty S. Kuntamukkula. Critical review of the foam rheology literature.Industrial and Engineering Chemistry Research, 26(2):318?325, 1987. ? pages 9, 12, 14[15] M.J. Rosen and J.T. Kunjappu. Surfactants and interfacial phenomena. John Wiley andSons., 2012. ? pages 10[16] Liguang Wang and Roe-Hoan Yoon. Effects of surface forces and film elasticity on foamstability. International Journal of Mineral Processing, 85(4):101?110, 2008. ? pages 10[17] W. B. Russel; D. A. Saville and W. R. Schowalter. Colloidal Dispersions. Cambridge,1989. ? pages 11[18] J. A. Krug. Oil gas journal, 73, 1975. ? pages 11[19] C. L. Wendorff and R. B. Ainley. Presented at the Annual Fall Meeting of the Society ofPetroleum Engineers, 73, 1981. ? pages 11[20] J. D. Lincicome. World oil, 28, 1984. ? pages 11[21] J. O. Sibree. The viscosity of froth. Trans. Faraday Soc, 28:325?331, 1934. ? pages 11[22] C.S.Grove; G.E.Wise; W.C.March and J.B.Gray. Viscosity of fire fighting foam. Ind. Eng.Chem., 43:1120?1122, 1951. ? pages 11[23] S.S.Marsden and S.A.Khan. The flow of foam through short porous media and apparentviscosity measurements. Soc. Pertrol Eng., 6:17?25, 1966. ? pages 12[24] R. J.; Stelson T. E Wenzel, H. G.; Brungraber. The viscosity of high expansion foam. J.Mater, 5:396?400, 1970. ? pages 12[25] J.T. Patton; H. Belkin; S. Holbrook and M. Kuntamukkula. Rheology of mobility-controlfoams. Soc. Petrol Eng., 23:456?460, 1983. ? pages 12[26] G.J.Hirasaki and J.B.Lawson. Mechanisms of foam flow in porous media: apparentviscosity in smooth capillaries. Annual Technical Conference and Exhibition of the Societyof Petroleum Engineers, San Francisco, CA, 1983. ? pages 12, 13, 14[27] A. M. Kraynik. Annual Meeting of the Society of Rheology, Evanston, IL, 1982. ? pages1393[28] H.M.Princen. Rheology of foams and highly concentrated emulsions. i. elastic propertiesand yield stress of a cylindrical model system. Journal of colloidal interface science, 1983.? pages 13[29] B Radvan and A.P.J Gatward. The formation of wet-laid webs by a foaming process.TAPPI, 55(5):748?751, 1972. ? pages 15, 18[30] M.K Smith and V.W Punton. The role of the forming process in determining the structureand properties of paper. Pulp and Paper Canada, 76(2):114?117, 1975. ? pages 15, 18[31] R.W. Tringham and M.A. Manager. New developments in Radfoam process. PaperTechnology, pages 288?294, 1974. ? pages 15[32] V.W Smith, M.K; Punton and A.G Rixson. The structure and properties of paper formed bya foaming process. TAPPI, 57(1):107?111, 1974. ? pages 15, 18[33] T. Lappalainen and J. Lehmonen. Determinations of bubble size distribution of foam-fibremixture using circular hough transform. NORDIC Pulp and Paper Research, 27(5):930?939, 2012. ? pages 15[34] H.; Salmela J.; Poranen, J.; Kiiskinen and J. Asikainen. Breakthrough in papermakingresource efficiency with foam. TAPPI journal, pages 1688?1695, 2013. ? pages 16[35] W.C.Hinds. Aerosol Technology: Properties, Behavior and Measurement of AirborneParticles. J. Wiley and Sons. New York. USA., 1982. ? pages 27[36] R.C.Brown. Air Filtration. Pergamon Press Inc., New York, USA., 1993. ? pages 27, 31[37] J. Seinfeld and P. Spyros. Atmospheric Chemistry and Physics: From Air Pollution toClimate Change. John Wiley and Sons, Inc., 2006. ? pages 27, 28[38] KW Lee and B.Y.H. Liu. On the minimum efficiency and the most penetrating particle sizefor fibrous filters. Air Pollution Control Association Journal, 30(4):337?381, 1972. ?pages 28[39] J.A. Pritchard. A guide to industrial respiratory protection. National Institute forOccupational Safety and Health. Cincinnati, OH., pages 5?8, 1976. ? pages 29[40] Clarke A. Rodman and Richard C. Lessmann. Automotive nonwoven filter media: theirconstructions and filter mechanisms. TAPPI, pages 161?168, 1988. ? pages 30[41] Jingliang Mao, Biljana Grgic, Warren H. Finlay, John F. Kadla, and Richard J. Kerekes.wood pulp based filters for removal of submicrometer aerosol particles. NORDIC Pulp andPaper Research, 23(4):420?425, 2008. ? pages 3094[42] Alan L. Macfarlane, John F. Kadla, and Richard J. Kerekes. High Performance Air FiltersProduced from Freeze-Dried Fibrillated Wood Pulp: Fiber Network Compression Due tothe Freezing Process. Industrial and Engineering Chemistry Research, 51(32):10702?10711, 2012. ? pages 30[43] C. Tseng, M. Yamaguchi, and T. Ohmori. Thermal conductivity of polyurethane foamsfrom room temperature to 20 K. Cryogenics, 37(6):305?312, 1997. ? pages 40[44] A.J. Loeb. Thermal conductivity: Viii, a theory of thermal conductivity of porousmaterials. American Ceramic Soc., 37(2):96, 1954. ? pages 40[45] J.C. Maxwell. A Treatise on Electricity and Magnetisim. Dover, Newyork, 1954. ? pages40[46] H.W.Russel. Principles of heat flow in porous insulation. American Ceramic Soc., 18(1):96, 1935. ? pages 40[47] G. Reichenauer; U. Heinemann and H.P. Ebert. Relationship between pore size and the gaspressure dependence of the gaseous thermal conductivity. Colloids Surf., pages 204?207,2007. ? pages 40[48] David S. Smith; Arnaud Alzina; Julie Bourret; Benoit Nait-Ali; Fabienne Pennec andNicolas Tessier-Doyen. Thermal conductivity of porous materials. Journal of MaterialsResearch, 28(17):2260?2272, 2013. ? pages 41[49] J.P. Arenas and M.J. Crocker. Recent trends in porous sound-absorbing materials. Soundand Vibration, Materials Reference Issue, 2010. ? pages 47[50] A.F. Seybert and D.F. Ross. Experimental determination of acoustic properties using a twomicrophone random excitation technique. Journal of the Acoustical Society of America,pages 1362?1370, 1977. ? pages 48[51] W. Lauriks; J. F. Allard and A. Cops. Inhomogeneous plane waves in layered media.Physical Acoustics, pages 425?431, 1991. ? pages 48[52] J.F. Allard. Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials.Chapman and Hall, 1993. ? pages 48[53] Vigran T.E.; Kelders L.; Lauriks W.; Leclaire P. and Johansen T.F. Prediction andmeasurements on the in?uence of boundary conditions in a standing wave tube.Acustica/Acta Acustica, pages 419?423, 1997. ? pages 48[54] Yvan Champoux and Michael R. Stinson. Journal of accoustical society of America, (2), .? pages 4895[55] Yvan Champoux and Michael R. Stinson. Journal of accoustical society of America, (4), .? pages 48[56] K. Attenborough. Acoustical characteristics of rigid fibrous absorbents and granularmaterials. Journal of accoustical society of America, 73:785?799, 1983. ? pages 48, 49[57] Y. Izumi; T. Iwase and R. Kawabata. A new measuring method for sound propagationconstant by using sound tube without any air spaces back of a test material. J. Internoise,98, 1998. ? pages 48[58] Y. Salissou and R. Panneton. Wideband characterization of the complex wave number andcharacteristic impedance of sound absorbers. Journal of accoustical society of America,128, 2010. ? pages 48[59] O. Doutres; Y. Salissou; N. Atalla and R. Panneton. Evaluation of the acoustic andnon-acoustic properties of sound absorbing materials using a three-microphone impedancetube. Journal of Applied Acoustics, 71(6):506?509, 2010. ? pages 48[60] M. Delany and E. Bazley. Acoustical properties of fibrous absorbent materials. Journal ofApplied Acoustics, 24(1):63?70, 1988. ? pages 49[61] Y. Miki. Acoustical properties of porous materials: modifications of delany-bazley models.J.Acoust. Soc. Jp., 11(1):19?28, 1990. ? pages 49, 53[62] M.A. Biot. Theory of propagation of elastic waves in a fluid-saturated porous solid. i. lowfrequency range, ii. higher frequency range. Journal of accoustical society of America, 28,1956. ? pages 49[63] D.L. Johnson; J. Koplik and R. Dashen. Theory of dynamic permeability and tortuosity influidsaturated porous media. Journal of Fluid Mechanics, 176:379?402, 1987. ? pages 49[64] J.F. Allard and Y. Champoux. New empirical equations for sound propagation in rigidframe fibrous materials. Journal of accoustical society of America, 9:3346?3353, 1992. ?pages 49[65] T.; N. Tsujiuchi Koizumi and A. Adachi. The development of sound absorbing materialsusing natural bamboo fibers. WIT Press, 2002. ? pages 49, 50[66] Y.E. Lee and C.W. Joo. Sound absorption properties of recycled polyester fibrous assemblyabsorbers. AUTEX Research Journal, 3(2), 2003. ? pages 49[67] M. Ren and F. Jacobsen. A method of measuring the dynamic flow resistance and reactanceof porous materials. Applied Acoustics, 39(4):265?276, 1993. ? pages 49[68] J. Conrad. Engineering Acoustics and Noise Control. Prentice-Hall, 1983. ? pages 5096[69] Y. Shoshani and Y. Yakubov. Use of nonwovens of variable porosity as noise controlelements. I.N.J, 2003. ? pages 50[70] M.A. Ibrahim and R.W. Melik. Physical parameters affecting acoustic absorptioncharacteristics of fibrous materials. Proceedings of the mathematical and physical societyof Egypt, 1978. ? pages 50[71] M. Coates and M. Kierzkowski. Technical-textiles-international. Acoustic Textiles -Lighter, Thinner and More Absorbent, 2002. ? pages 50[72] A. Tamayol and M. Bahrami. Analytical determination of viscous permeability of fibrousporous media. International Journal of Heat and Mass Transfer, 52(9):2407?2414, 2009.? pages 56[73] F. E. Teruel. Characterization of a porous medium employing numerical tools:Permeability and pressure-drop from darcy to turbulence. International Journal of Heatand Mass Transfer, 52(25):5878?5888, 2009. ? pages 56[74] D. Shou; J. Fan and F. Ding. Hydraulic permeability of fibrous porous media. InternationalJournal of Heat and Mass Transfer, pages 4009?4018, 2011. ? pages 56[75] C.N. Davies. The separation of airborne dust and particles. Proceedings of Institute ofMechanical Engineers, London, 1952. ? pages 56, 82, 83, 88[76] J. Drummond and M. Tahir. Laminar viscous flow through regular arrays of parallel solidcylinders. International Journal of Multiphysics Flow, 1984. ? pages 82[77] A. Sangani and A. Acrivos. Laminar viscous flow through regular arrays of parallel solidcylinders. International Journal of Multiphysics Flow, pages 193?206, 1982. ? pages 56,82[78] C. Cercignani. The Boltzmann equation and its applications. Springer-Verlag, New York,1988. ? pages 60[79] S. Succi. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. OxfordScience Publications. ? pages 60, 62, 63, 65, 66, 68, 69, 74, 87[80] M. Kardar. Statistical Physics of Particles. Cambridge University Press, 2007. ? pages 60[81] D.A. Wolf-Gladrow. Lattice-gas cellular automata and lattice Boltzmann models: anintroduction. 2000. ? pages 61, 62[82] E. R. Berlekamp; J. H. Conway and R. K. Guy. Winning ways for your mathematical plays.Academic Press, New York, 1984. ? pages 6297[83] M. Garder. The fantastic combinations of conway?s new solitaire game of ?life?. ScienceAmerica, 1970. ? pages 62[84] M. Garder. Geometric fallacies: hidden errors pave the road to absurd conclusions. ScienceAmerica, 1971. ? pages[85] M. Garder. On cellular automata, self-reproduction, the garden of eden and the game ?life?.Science America, 1971. ? pages 62[86] G.D. Doolen; U. Frisch; B. Hasslacher; S. Orszag and S. Wolfram. Lattice Gas Methodsfor Partial Differential Equations. Academic Press, New York. ? pages 62[87] G. McNamara and G. Zanetti. Use of the boltzmann equation to simulate lattice gasautomata. Physical Review Letter, 61:2332, 1988. ? pages 62[88] F. Higuera and J. Jimenez. Boltzmann approach to lattice gas simulations. Europhys. Lett.,9:663, 1989. ? pages 62, 63[89] P. Bhatnagar; E. Gross and M. Krook. A model for collision process in gases i: Smallamplitude processes in charged and neutral one component system. Phys. Rev., 94:511,1954. ? pages 63[90] Y. Qian; D. d?Humieres and P. Lallemand. Lattice bgk models for the navier-stokesequation. Europhys. Lett., 17:479, 1992. ? pages 63[91] Martin Hecht and Jens Harting. Implementation of on-site velocity boundary conditions ford3q19 lattice boltzmann simulations. Journal of Statistical Mechanics: Theory andExperiment, 10:1018, 2010. ? pages 64[92] S. Chapman and T. Cowling. Mathematical Theory of Non-Uniform Gases. CambridgeUniversity Press. ? pages 65[93] X. Shan and H. Chen. Lattice boltzmann model for simulating flows with multiple phasesand components. Phys. Rev. E, 47:1815?1819, 1993. ? pages 67[94] Martin Hecht and Jens Harting. Lattice boltzmann method for fluid flows. Annual Reviewof Fluid Mechanics, 30:329?364, 1998. ? pages 67[95] M Rebai, F Drolet, and D Vidal. A Lattice Boltzmann approach for predicting the captureefficiency of random fibrous media. Asia Pacific Journal of Chemical Enginering, pages29?37, 2011. ? pages 67, 82, 83[96] O. Filippova and D. Hanel. Lattice-bgk model for low mach number combustion. Int. J.Mod. Phys. C, 9:1439, 1998. ? pages 6798[97] S. Chen; H. Chen; D. Martinez and W. Matthaeus. Lattice boltzmann model for simulationof magnetohydrodynamics. Phys. Rev. Lett., 67:3776?3779, 1991. ? pages 67[98] O. Filippova; S. Succi; F. Mazzocco; C. Arrighetti; G. Bella and D. Hanel. Multiscalelattice boltzmann schemes with turbulence modeling. Int. J. Mod. Phys. C, 170:812829,2001. ? pages 67[99] M.L.R; D.B. Bingham Thomas and M. Pourkashanian. Prediction of the permeability offibrous porous media using the lattice boltzmann method in conjuction with coarsenumerical lattices. Open Transport, pages 80?89, 2010. ? pages 67[100] R. Przekop and L. Gradon. Non-steady-state aerosol filtration in nanostructured fibrousmedia. Philosophical transactions. Series A, Mathematical, physical, and engineeringsciences, 369:2476?2484, 2011. ? pages 67[101] O. Filippova and D. Hanel. Grid refinement for lattice-bgk models. Journal ofComputational Physics, 14:219?228, 1998. ? pages 68[102] D. Yu; R. Mei; L. Luo and W. Shyy. Viscous flow computations with the method of latticeboltzmann equation. Progress in Aerospace Sciences, 39:329?367, 2003. ? pages 68, 69,75[103] Urpo Aaltosalmi. Fluid flow in porous media with the Lattice-Boltzmann method. PhDthesis, University of Jyvaskyla, 2005. ? pages 68[104] Chen Peng. The lattice boltzmann method for fluid dynamics: Theory and applications.,2013. ? pages 68[105] Q. Zou and X. He. On pressure and velocity boundary conditions for the lattice boltzmannbgk model. Physics of Fluids, 9:1591, 1997. ? pages 68, 72, 73, 74, 75[106] T. Inamuro; M. Yoshina and F. Ogino. A non-slip boundary condition for lattice boltzmannsimulations. Physics of Fluids, 7:2928?2930, 1995. ? pages 68[107] P. A. Skordos. Initial and boundary conditions for the lattice boltzmann method. Phys. Rev.E, 48:4823?4842, 1993. ? pages 69[108] J. Verschaeve and B. Muller. A curved no-slip boundary condition for the lattice boltzmannmethod. Journal of Computational Physics, 2010. ? pages 69[109] C. Zheng Z. Guo and B. Shi. An extrapolation method for boundary conditions in latticeboltzmann method. Phys. Fluids, 14:2007?2010, 2002. ? pages 69[110] D. Hanel O. Filippova. Boundary setting and local grid refinement for lattice-bgk models.International Journal of Modern Physics C, 9:1271?1279, 1998. ? pages 69, 71, 7299[111] R. Mei; L.S. Luo and W. Shyy. An accurate curved boundary treatment in the latticeboltzmann method. Journal of Computational Physics, 155:307?330, 1999. ? pages 69,71, 72[112] M. Bouzidi; M. Firdaouss and P. Lallemand. Momentum transfer of a lattice boltzmann uidwith boundaries. Phys Fluids, 2001. ? pages 69[113] James D. Sterling and Shiyi Chen. Stability analysis of lattice boltzmann methods. Journalof Computational Physics, 123(1):196?206, 1996. ? pages 74, 75[114] D. DHumieres; I. Ginzburg; M. Krafczyk; P. Lallemand and L. Luo. Multiple relaxationtime lattice boltzmann models in three dimensions. Philosophical transactions, 360:437?51, 2002. ? pages 75[115] P. Lallemand and L.S. Luo. Theory of the lattice boltzmann method: Dispersion,dissipation, isotropy, galilean invariance, and stability. Physical Review E, 61:6546?6562,2000. ? pages 75100