w , .q." 24...... A") firfi§;.Lm , 11:531....» ‘3 vi? a. 3,5...Hu J..: . a 1. 4. . J. :3- Qfi arui 45m}? ”an”. 3a. b 3.36“}: 21:. 12...: 1.9.1.... 2K1: $ 2 : (8;: \\ .L, :3 1.123.21‘ i: L." .1. {21 .w... :39. 5 1: 3. .l...auvnu_. u...uw;_.u.ah“~nl.... ..,,..‘.n_¥av.a5_(.l_I. 9 (it ;H g. .32. 1 . . .3, . . ‘ V , .91; m... x :1: IE. . l , cg”! . :vh‘w : nu g“! ”if .K ? ‘ M . . .. . :t. 9 J x .r ,n in! 3...... . . “HIV...“ .i V . . ...~§ v. » xii a. ; This is to certify that the .7 dissertation entitled Numerical Simulation of Cellular Natural Convection in a Wide Porous Box to Evaluate Thermodynamic Wavenumber Prediction Theories for an Infinite Porous Layer presented by Joseph Bayne Schroeder has been accepted towards fulfillment of the requirements for the Ph. D. degree in Mechanical Engineering @WJMJ Major Professor’s Signature ZMG [y 5 Q 006 Date MSU is an Affirmative Action/Equal Opportunity Institution __ *M_ _-__~ .——~-———————-- LIBRARY .. 3:24“; State University PLACE IN RETURN BOX to remove this checkout from your record. To AVOID FINES return on or before date due. MAY BE RECALLED with earlier due date if requested. DATE DUE DATE DUE DATE DUE 2/05 p:/C|RC/DateDue.indd-p.1 NUMERICAL SIMULATION OF CELLULAR NATURAL CONVECTION IN A WIDE POROUS Box TO EVALUATE THERMOOYNAMIC WAVENUMBER PREDICTION THEORIES FOR AN INFINITE HORIZONTAL POROUS LAYER By Joseph Bayne Schroeder A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mechanical Engineering 2006 IN} Wher that ii move whicr Theor the we unpro few. T waver entIOp throng exDies Gauilib For thi: IBCIang BIInkm ABSTRACT NUMERICAL SIMULATION OF CELLULAR NATURAL CONVECTION IN A WIDE POROUS BOX TO EVALUATE THERMODYNAMIC WAVENUMBER PREDICTION THEORIES FOR AN INFINITE HORIZONTAL POROUS LAYER By Joseph Bayne Schroeder When a layer of fluid or fluid-saturated porous media is heated from below so that the temperature gradient exceeds a critical value, the fluid begins to move in a cellular convective pattern, known as Rayleigh-Bénard convection, which significantly increases the rate of heat transferred across the layer. Theories to predict the size of the convective cells, which are characterized by the wavenumber of convection because of their periodic nature, are still unproven; and experimental measurements of the wavenumber have been few. This work seeks to verify two thermodynamic theories to predict the wavenumber of convection: one theory based on extremum of the rate of entropy generation caused by both heat transfer and viscous fluid flow through the medium, and another theory employing a stability functional expressing the ”generalized excess entropy production," based on the non- equilibrium thermodynamic theory of stability and fluctuations. For this work numerical simulations are performed for a wide, horizontal, rectangular enclosure filled with fluid and porous media. The modified Darcy- Brinkman-Forchheimer momentum and continuum-model thermal energy equafi Inethc thch mode vekaci entrot funciii FouRe vaRaU Nume ”Umbl thata; aSpect Brinkn Nussel anub‘ enUop‘ fUllctio Compa DOrOUS i”dICate equations are solved with a mixed Fourier-Galerkin and finite-difference method. The nonlinearities in the momentum equation are transformed into the Fourier space; and the equations are quasi-linearized about each diagonal mode. The discretized equations are then solved by iteration, alternating the velocity and temperature solutions. The stream function, Nusselt number, entropy generation rate, and generalized excess entropy production functional are calculated from the converged solutions in the discretized Fourier space. The wavenumber of convection is then measured from the variations in the flow velocity field. Numerical solutions with successively wider boxes demonstrate Nusselt number and wavenumber dependence on the porous media Rayleigh number that approaches the infinite layer behavior in the central region of boxes with aspect ratios of about 15 for Rayleigh numbers between 40 and 200. The Brinkman and Forchheimer momentum terms are seen to only affect the Nusselt number and wavenumber for extreme values. Similar numerical simulations of an infinite horizontal layer were performed that employed the entropy generation rate and generalized excess entropy production stability functional theories to determine the wavenumber for the system. Comparisons with the Nusselt number and wavenumber dependence on the porous media Rayleigh number from the porous box simulations of this work indicate the best agreement with the entropy generation rate theory. Copyfightby Joseph Bayne Schroeder 2006 To the E Preparir an0ppc To the Engineering students of Olivet Nazarene University who, while preparing for lives of service to God, humanity, and their profession, give me an opportunity to fulfill my purpose. Thesuc assistan PmfieSSI acadeni compar contlnw Theren‘ Bénard, mscussi Mmhga academy MY Colle “mhedp Reams, I Brow”, I e”COUrag SDTRUaI‘ Acknowledgements The success of this work owes much to many who provided support and assistance, and to whom I owe a debt of great gratitude. Professor Craig Somerton has been a mentor and advisor on matters academic, technical, and educational, as well as a friend and occasional golf companion, for close to fourteen years. I thank him for his patient and continued support and encouragement. The remainder of the Doctoral Committee: Professors James Beck, André Bénard, and Charles MacCluer, provided useful advice, review, and discussions. I also thank the Department of Mechanical Engineering at Michigan State University for providing opportunities and support early in my academic career. My colleagues at and around Olivet Nazarene University, including (but not limited to) Mike Morgan, lvor Newsham, Rod Korthals, Eric Erickson, Max Reams, Fran Reed, Gary Streit, Brock Schroeder, Daniel Green, Stephen Brown, Larry Vail, and Frank Garton; have been a great source of encouragement, listening ears, and sound advice. I thank them all for their spiritual and scholarly support. I would also like to thank the GNU Department vi cfi E ancl Rte 2 Daci than! seeir Ifinah Jesus (Nder loan Creatii of Engineering and the Office of Academic Affairs, who provided resources and support enabling me to complete this work. My wife, Barb, and my dear children: Lucy, Adam, and Megan, have stood by me and smiled upon me, while often suffering the absence of husband and Dad for the seemingly endless hours it takes to complete a dissertation. I thank them and state my unending devotion. I hope they can now enjoy seeing Daddy wear his funny hat. Finally, and most of all, I must give thanks and praise to God the Father and to Jesus the Christ, who created all things in His perfect wisdom, who brings order out of the chaos, who continues to shape the world by His gracious, loving, and redeeming hand, and who reveals Himself and the works of His creation to us in beautiful and amazing ways. Amen vii Ufldh URMR Nomenc 3 l o g o I u I I o A a _ - - - _ _ - N 9? A! A\A\A\‘\ [\J M 0‘ UT A (A) (.3 r.\ r.‘ r.\ NN C301 Table of Contents List of Tables .................................................. xi List of Figures ................................................. xi Nomenclature ................................................ xiii 1 Introduction ............................................. 1 1.1 Description of Problem ................................ 1 1.2 Literature Review ..................................... 3 1.2.1 Porous Media Flow Analysis 1.2.2 Cellular Convection in Fluid Layers 1.2.3 Cellular Convection in Porous Layers 1.2.4 Entropy Generation Theories 1.3 Scope of Dissertation ................................. 20 1.4 Outline of the Dissertation ............................. 21 2 Governing Equations ...................................... 22 2.1 Assumptions ........................................ 22 2.2 Governing Conservation Equations ...................... 24 2.2.1 Mass Conservation (Continuity) 2.2.2 Momentum Conservation 2.2.3 Conservation of Thermal Energy 2.2.4 Boundary Conditions 2.3 Non-Dimensionalization and Scaling ..................... 28 2.3.1 Mass Conservation 2.3.2 Conservation of Momentum 2.3.3 Energy Equation 2.3.4 Boundary Conditions 2.4 Reduction of Mass and Momentum Equations ............. 35 2.5 Entropy Generation Rate Equation ...................... 37 2.5.1 Entropy Generation Rate Calculation for Porous Media Flow 2.5.2 Non-Dimensionalization 2.6 Generalized Excess Entropy Production Rate Functional ..... 43 viii 3 4 5 Numerical Solution Method ................................ 48 3.1 Galerkin Forms of Governing Equations .................. 49 3.1.1 Galerkin Ansatz 3.1.2 Mass Continuity Equation 3.1.3 Galerkin Form of the Momentum Equation 3.1.4 Expansion of the Nonlinear Forchheimer Momentum Term 3.1.5 Galerkin Form of the Thermal Energy Equation 3.1.5.1 Mode Zero Equation 3.1.5.2 Higher Modes Equation 3.2 Modal Ouasi-Linearization ............................. 59 3.3 Finite Difference Equations ............................ 62 3.3.1 Momentum Equation 3.3.2 Zeroth Mode Energy Equation 3.3.3 Modal Energy Equations 3.3.4 Expansion of the Nonlinear Term 3.4 Calculation to Recover Horizontal Velocity ................ 70 3.5 Post-Processing ..................................... 70 3.5.1 Computing the Pressure Field 3.5.2 Stream Function Calculation 3.5.3 Nusselt Number Calculation 3.5.4 Entropy Generation Rate 3.5.5 Generalized Excess Entropy Production Rate 3.6 Numerical Solution Algorithm .......................... 79 Results ................................................. 82 4.1 Effect of Box Width .................................. 82 4.1.1 Onset of Convection 4.1.2 Nusselt Number 4.1.3 Wavenumber Measurement Method 4.1.4 Wavenumber Results 4.2 Non-Linear Results ................................... 99 4.2.1 Forchhiemer Effect 4.2.2 Brinkman Effect 4.3 Comparison with Wavenumber Prediction Theories ....... 102 4.4 Numerical Validation ................................ 107 Conclusions and Recommendations ......................... 114 5.1 Summary of Contributions ........................... 114 5.2 Conclusions ....................................... 115 5.3 Recommendations for Future Study .................... 116 Appendh Appendh Appendh Appendi lfibfiogr Appendix A: Galerkin Inner Products of Sine and Cosine Bases ......... 118 Appendix B: Triple-Inner Products ................................ 121 Appendix C: Galerkin Integrals for Entropy Generation Rate ........... 124 Appendix D: Galerkin Integrals for Excess Entropy Production Rate Functional ................... 128 Bibliography ................................................. 136 Tabl Tabl Note Figur Flgur Figur‘ Figur Flgun Figuri Figun Figure Figure Figure Figure Figure Figure Figure I:lgure Figure Figure Flgure I:Igure FIgure Figure Flgure Flgure Figlire FIQUIe FiQure Flgure Flgure Figure FlSjure Table 4.1 Table 4.2 Note: Some Figure 2.1 Figure 2.2 Figure 4.1 Figure 4.2 Figure 4.3 Figure 4.4 Figure 4.5 Figure 4.6 Figure 4.7 Figure 4.8 Figure 4.9 Figure 4.10 Figure 4.11 Figure 4.12 Figure 4.13 Figure 4.14 Figure 4.15 Figure 4.16 Figure 4.17 Figure 4.18 Figure 4.19 Figure 4.20 Figure 4.21 Figure 4.22 Figure 4.23 Figure 4.24 Figure 4.25 Figure 4.26 Figure 4.27 Figure 4.28 List of Tables Ranges for Non-Linear Parameters ...................... 99 Results from Porous Layer and Porous Box Simulations . . . . 103 List of Figures images in this dissertation are presented in color. Diagram of Physical System ........................... 24 Non-Dimensionalized Porous Box Heated From Below ...... 34 Streamlines for Ram=50, y=3 .......................... 84 Streamlines for Ram=75, y=3 .......................... 84 Streamlines for Ram=100, y=3 ......................... 84 Isotherms for Ram=50, y=3 ............................ 85 Isotherms for Ram=75, y=3 ............................ 85 Isotherms for Ram=100, y=3 ........................... 85 Streamlines for Ram=50, y=8 .......................... 86 Streamlines for Ram=75, y=8 .......................... 86 Streamlines for Ram=100, y=8 ......................... 86 Isotherms for Ram=50, y=8 ............................ 87 Isotherms for Ram=75, y=8 ............................ 87 Isotherms for Ram=100, y=8 ........................... 87 Streamlines for Ram=50, y=15 ......................... 88 Streamlines for Ram=75, y=15 ......................... 88 Streamlines for Ram=100, y=15 ........................ 88 Streamlines for Ram=120, y=15 ........................ 88 Isotherms for Ram=50, y=15 ........................... 89 Isotherms for Ram=75, y=15 ........................... 89 Isotherms for Ram=100, y=15 .......................... 89 Isotherms for Ram=120, y=15 .......................... 89 Bounds of Critical Rayleigh Number for Convective Motion . . 91 Heat Transfer for Linear Simulations ..................... 92 Central Average Wavenumber vs. Rayleigh Number ........ 97 Overall Average Wavenumber vs. Rayleigh Number ........ 97 Comparison of Overall and Central Average Wavenumbers . . 98 Effect of the Forchheimer Parameter, M ................. 100 Effect of the Brinkman Parameter, B .................... 101 Nusselt Number Comparison for Wavenumber Prediction Schemes .......................................... 104 xi Fr Fig Fig Figure 4.29 Figure 4.30 Figure 4.31 Figure 4.32 Figure 4.33 Figure 4.34 Figure 4.35 Wavenumber vs. Rayleigh Number Under Various Prediction Schemes .......................................... 106 Nusselt Number Convergence with Number of Galerkin Modes ................................................. 111 Nusselt Number Convergence with Finite Difference Mesh Size ................................................. 111 Wavenumber Convergence with Number of Galerkin Modes 112 Wavenumber Convergence with Finite Difference Mesh Size 112 Maximum Vertical Velocity Component Convergence with Number of Galerkin Modes ........................... 113 Maximum Vertical Velocity Component Convergence with Finite Difference Mesh Size ................................ 113 xii 4 him ’I . .1 ‘ldmt‘ C. 00 . I... rum 3%... B .I o m m WI.» 9 II m m I .m. 11*: v. 1' C ~IABCCCI€DDI DD FFFFF.Fm.ggGGhHHvAIKmLmUIJn\.nDIDIP.qRu: Bl.i9 BZJ ’ BJJ 3*”. 'I ~b| ‘01 '71: 551» ‘7 5 U ZDSBBSw-pgsgk‘cgcam oqla‘q \ -- .9 Nomenclature unit vectors in x-, y-, and 2- directions wavenumber of convection wavenumber coefficient of k-th galerkin mode basis function inner product of cosine and sine trigonometric triple inner products nonlinear momentum term coefficient of Wm collected Brinkman momentum term coefficients specific heat effective specific heat of porous medium finite difference coefficients for discretized momentum equations mass-specific internal energy finite difference coefficients for zeroth mode thermal energy equafions finite difference coefficients for higher mode thermal energy equafions collected nonlinear advective thermal energy transport terms m-th mode inner product of nonlinear advective terms remainder of quasi-linearized advective thermal energy terms RHS for discretized momentum equations gravitational acceleration vector magnitude of gravitational acceleration total entropy generation rate components for entropy generation rate integral finite-difference grid spacing k-th Fourier expansion mode of velocity magnitude term i-th finite difference node of k-th mode of velocity magnitude number of Galerkin modes effective thermal conductivity of porous medium vertical thickness of porous enclosure Forchheimer coefficient collected Forchheimer momentum term coefficients Forchheimer exponent number of vertical finite difference mesh nodes non-dimensional pressure (dimensional) pressure k-th Galerkin mode of pressure heat flux vector average heat flux at wall xiii u I O ‘ W .0 m m. m. . a . ~ ~ we E. . . . .T. .I Q0. 0. 0. DA SS 5 S 1..in T.n§/CT.HAT;IT.LH.U .vwl £62 §Q>© 330’”!st * 0 to OS” gen CI.) gen XS N.“ collected nonlinear Forchheimer momentum terms m-th mode inner product of nonlinear momentum terms remainder of quasi-linearized nonlinear momentum terms discretized quasi-linearized nonlinear momentum terms porous media Rayleigh number mass-specific local entropy components for ‘I’ integral calculation volumetric rate of entropy generation dimensionless entropy generation rate dimensionless time dimensional time non-dimensional scaled temperature (OsTs1) dimensional temperature reference temperature for Boussinesq approximation temperature at cold horizontal boundary (top wall) temperature at heated horizontal boundary (bottom wall) temperature difference TH' — Tc' k-th Galerkin mode of temperature i-th finite difference node of k-th Galerkin mode of temperature Darcian velocity component in depth dimension Darcian velocity vector (interstitial) fluid velocity vector horizontal Darcian velocity component, or specific volume k-th Galerkin mode of horizontal Darcian velocity i-th finite difference node of k-th mode of horizontal velocity vertical Darcian velocity component horizontal width of porous enclosure i-th finite difference node of k-th Galerkin mode of vertical velocity k-th Galerkin mode of vertical Darcian velocity depthwise space dimension horizontal space dimension RHS for discretized thermal energy equations vertical space dimension xiv Greek L- R ”‘i §owovv--<‘CD 1: '»x m PO31:§2>8 _-3 up 3 P Subscript 1U, m fuperscril Greek Letters 3 fl PPR: region-<11”: on >3 thermal diffusivity of porous medium volumetric coefficient of thermal expansion aspect ratio: W/L Newton iteration damping factor for temperature Newton iteration damping factor for velocity porosity of porous medium velocity magnitude nonlinear momentum term permeability of porous medium wavelength of convection cells, measured in the central rolls wavelength of convection cells, measured across width of box dynamic viscosity of fluid Brinkman-effective viscosity of porous medium mass density of fluid at Boussinesq reference temperature To mass density of fluid mass density of porous medium porous media viscosity ratio lJm/“f irreversibility ratio stream function generalized entropys tability functional viscous dissipation function porous media thermal storage capacity ratio generalized excess entropy production functional Non-Dimensional Numbers Subscripts k, l, m 1 m f Superscripts *- Darcy number exponent-modified Forchheimer number buoyant heat transfer number Nusselt number porous media Prandtl number Rayleigh number porous media Rayleigh number Galerkin mode number finite difference node number porous media property fluid property dimensional value, or solution value prior to Newton step damping XV a M «C i d D 0 V fly 5 WITCIJ Mathematical Symbols 44006::- ~cu ‘HM: III partial derivative operator ordinary derivative operator total derivative operator material derivative operator variation operator gradient or divergence operator Laplacian operator definition vector magnitude (Euclidian norm) summation integral xvi 1.1 Thksvv inapo fimnb chfical tmnper the hea howeve theory 1 charact hequen Th€ihel inseekh mennod l- TI he and 2 As Du) flu Chapter 1: Introduction 1.1 Description of Problem This work is a numerical study of the behavior of Rayleigh-Bénard convection in a porous medium. When a layer of porous media filled with fluid is heated from below (and cooled above) with a temperature difference exceeding a critical value, the fluid begins a cellular convective motion caused by temperature-induced density gradients. This motion significantly enhances the heat transfer rate across the layer. The characteristics of this motion, however, are still not fully understood. In particular, there remains no proven theory to predict the size of the convection cells. This size is typically characterized by the wavenumber of convection: a non-dimensional frequency, because of the periodic nature of the convection cells. The theory employed in this work will use the second law of thermodynamics in seeking to provide closure to the wavenumber question. Two thermodynamic theories will examined: 1. That the rate of entropy generation from the combined effects of heat transfer and of viscous fluid flow is minimized in the steady state, and 2. A stability functional, ‘1’, related to the generalized excess entropy production of the non-equilibrium thermodynamic theory of fluctuations, is minimized at the stable steady state. In t hor mor will the I conc the c layer ofa i deter the S) existii The in 3. T, CC kn In this work, non-linear numerical simulations will be performed for a wide horizontal, rectangular enclosure filled with fluid-saturated porous media. The momentum and energy equations describing the flow and temperature fields will be solved, and from these solutions the heat transfer behavior and size of the convective planform will be determined. Subsequent solutions will be conducted by increasing the width of the box, until the behavior of the fluid in the central region of the porous box approaches that of the infinite horizontal layer. These results will then be compared with similar numerical simulations of a horizontal porous layer, which require these thermodynamic theories to determine the wavenumber to adequately mathematically describe and solve the system. The results of both simulations will also be compared with existing published experimental and numerical data. The intended results of this study are: 1. To determine whether the thermodynamic theories adequately predict the wavenumber of convection for Rayleigh-Bénard convection in a porous medium. 2. To compare, and possibly reconcile the results of, two seemingly distinct thermodynamic theories. 3. To generate additional wavenumber and heat transfer data for natural convection in porous layers and enclosures, adding to the body of knowledge in this area of heat transfer. 4, To no 1.2 Li The beh. medium been on researcl undergr nuclear This se analysi are av; Cheng reviev. VWfila body: IBIEVE COurs this 0 4. To develop and employ novel numerical solutions for heat transfer and non-linear flows in porous media. 1.2 Literature Review The behavior of a layer of pure fluid or of a layer of fluid saturated porous medium which undergoes convective motion when heated from below has been one of the classic problems in fluid mechanics and heat transfer research. The knowledge of the behavior of porous layers can be applied to underground geothermal beds, horizontal layers of thermal insulation, and nuclear reactor cores, to list just a few applications. This section will begin with a review of principles of porous media flow analysis relevant to this dissertation. Several review articles, handbooks, etc. are available which more thoroughly cover flows in porous media, such as Cheng (1978), Bejan (1987), Kaviany (1991), and Nield and Bejan (1998). A review of the Rayleigh-Bénard problem in fluid layers will follow, and continue with a review of research in the porous layer problem. There exists a vast body of research in these areas, and the review here is only a sample of the relevant work. A number of additional references which were gathered in the course of this research are also listed in the Bibliography at the conclusion of this document. 1.2.1 For The stuch have beg between Law, for Darcy's flow ur BCCOUI Flow ' by th. Law. conc 10m bee A v Ch; ex Ini 1.2.1 Porous Media Flow Analysis The study of flow through a porous medium is traditionally considered to have begun with Darcy, who determined the empirical linear relationship between the pressure gradient and the fluid flow rate, now called Darcy's Law, for public fountains fed by underground pressure (Darcy 1856). Darcy’s Law has its limitations. It fails to accurately describe the physics of the flow under many situations, including rapid fluid flow rates, and does not account for observed flow dependence on several additional parameters. Flow variations depending on the particle size in the medium, characterized by the Darcy number, are described by the Brinkman extension to Darcy’s Law. This term also allows for a characterization of no-slip boundary conditions at the solid boundaries, which would over-specify the Darcy Law formulation. (Neale and Nader, 1974) The use of the Brinkman extension has been often debated in the literature by Nield (1983, 1985, 1991, 1995). A variation with fluid properties and inertial effects at higher flow velocities is characterized by the Prandtl number, and is described by the Forchheimer extension to Darcy's Law. A Prandtl number effect which accounts for the inertial properties of the flowing fluid was presented by Somerton (1983), and later modified by Georgiadis and Catton (1986), derived by order of magnitude analysis from Navier-Stokes flow over parallel cylinders. This leads exten ofthe Jonss medii of Pre condi meas devia ofthe scattr "effec Carm Forcf theni abou1 SOme pUbHs BOlle: nUmb infinlt. leads to inclusion of the Kozeny-Carman term, or Forchheimer inertial extension to the momentum equations, to mathematically capture the effect of the Prandtl number variations. Jonsson and Catton (1987) performed experiments with various packing media and fluids in a cylinder with variable heights to demonstrate the effect of Prandtl number on natural convection in porous media. Effective thermal conductivities and Forchheimer equation permeability values were also measured. The use of values calculated from correlations was shown to have deviated greatly from actual measurements, and it is suggested that the use of these calculations to establish media properties could be cause for much scatter of the results in the literature. Jonsson and Catton defined an ”effective Prandtl number" as the product of porous medium and Kozeny- Carman numbers, based on the appearance of those parameters in the Forchheimer-extended momentum equation. Their Nusselt number results then exhibited good correlation along Prandtl number curves. Some scatter about the curves was attributed to a Darcy number dependence, as well as some experimental uncertainties. The data also compared favorably with published results of Georgiadis and Catton (1986) and Combarnous and Bories (1975). For effective Prandtl numbers greater than 0.1, the Nusselt number results exhibited no Prandtl number dependence, and matched the infinite Prandtl number case (exhibiting no inertial effects). Ber: Jon Prar that doe how num ITIUC The l veloc magr inves; to det flows, eXDOn ShKfles apDTO) 12.2 ( One of l beenihi Because the effective Prandtl number varies inversely with the layer thickness, Jonsson and Catton determined that “...no matter how large the medium Prandtl number is, there will always be a layer thickness that is large enough that inertial effects will become important.” This effect of the layer thickness does not apply to the case of Rayleigh-Bénard convection in a pure fluid, however. This dependence on plate spacing, as well as medium Prandtl number, is also suggested as a cause for scatter in published data which is so much more pronounced in the porous media than in the pure fluids. The Forchheimer extension has evolved in several forms. Often used as a velocity-squared term, the more favored vector representation is the velocity magnitude, multiplied by the velocity vector component. Genik (1998) investigated both linear and non-linear flow in several different porous media to determine the permeability and Forchheimer coefficients for non-linear flows. The Forchheimer term was further modified to include a variable exponent for the velocity magnitude factor. Genik performed parameter studies for this exponent, and determined that the exponent should be approximately 1.8, depending on the flow conditions. 1.2.2 Cellular Convection in Fluid Layers One of the fundamental problems in the study of convection heat transfer has been the characteristics of a layer of fluid which undergoes natural convective motion when under a temperature gradient. Thest gradie with Li ceHula Howev shown additio Sternlir phenor heat tra Much VI. configm Charactg become fluid lay. alld also The part. CO’Il’é‘Ct/i Patterns! numbets. The study of the stability of a layer of fluid under an adverse temperature gradient — to determine the conditions which will initiate fluid motion —- began with Lord Rayleigh (1916). Bénard (1900) had previously discovered a natural cellular convective phenomenon in a fluid layer similarly heated from below. However, the convective cells observed in these experiments were later shown to have been caused by surface tension effects at the free surface, in addition to the thermal gradient, as studied by Pearson (1958), Scriven and Sternling (1964), Nield (1964), and Scanlon and Segel (1967). Yet the phenomenon, which has become an essential part of the study of convective heat transfer, continues to bear the names of both Rayleigh and Bénard. Much work has been done on the issue of the stability of fluid in this configuration, to determine the degree of the adverse temperature gradient, characterized by the Rayleigh number, required to cause the fluid layer to become unstable and undergo continuous motions. This issue of stability in a fluid layer is explored in a full section of the book by Chandrasekhar (1961) and also the two volumes by Joseph (1976). The pattern configuration of the resulting pattern motion, termed the convective planform, has been an issue of study, and the size of repeated cell patterns, characterized by the wavenumber of convection, at Rayleigh numbers above the critical value has been a fertile area of research. Willis cells Béna diam may whicl Seeki meas silicor gmph to-heig compa infinhe The res dimens MMed Cmnan e)(Derim Rome) With inc” With II'ICrE Willis, Deardorff, and Somerville (1972) measured the diameter of convection cells from photographs of air, water, and silicone oil undergoing Rayleigh- Bénard convection at supercritical Rayleigh numbers. They found that the diameter "depends uniquely on [the Rayleigh number],” and that this variation may explain the disagreements found with other heat transfer predictions which use a fixed wavenumber of convection. Seeking wavenumber data, Luijkx and Platten (1982) took experimental measurements of Rayleigh-Bénard convection cells for Dow Corning 200 silicone oil in a thin fluid layer heated from below using a simple shadow graph method. The wavelength of the cells were measured at several width- to-height ratios near the critical Rayleigh numbers. Their results were compared to three theories predicting behavior of the Rayleigh-Bénard cells: infinite longitudinal rolls, finite transverse cells, and three-dimensional cells The results of the experiments agree well with theory based on three- dimensional perturbations, with differences attributed to measurement error in the depth and the fact that the layer is actually finite, and therefore must contain an integer number of cells. The measured temperature gradient in the experiment also must have been slightly supercritical, so that the cells could actually be observed to develop. The wavenumber was observed to decrease with increased temperature gradient, and to decrease and then increase again with increasing width (or width to height ratio). Based on the wavenumber re th Thr WBV MCDOn thermOI in DIEUII layer of l thermod EntrOpy I Was Used onetell (3 Shown to prOdlictio results, the authors assumed that the cellular flow systems were actually three-dimensional. Three theories seem to have been used to attempt to predict the proper wavenumber of convection. 1. Using the critical wavenumber which establishes at the onset of convection. 2. Selecting the wavenumber that maximizes the heat transfer rate or Nusselt number. 3. Using a stability functional associated with the ”generalized excess entropy production rate” of the system, which seems to exhibit a linear characteristic at preferred wavenumbers. McDonough (1980), in his dissertation work, used the non-equilibrium thermodynamic theory of Glansdorff and Prigogine (1971) to provide closure in predicting the wavenumber for natural convection in an infinite horizontal layer of fluid. This theory accounted the convection of the fluid with the usual thermodynamic entropy production in a so-termed ”generalized excess entropy production." A mixed Galerkin-finite difference numerical scheme was used to solve the Navier-Stokes and thermal energy equations over a one-cell period. Although a theoretical solution for wavenumber selection was shown to be possible by minimizing the generalized excess entropy production functional, McDonough instead approximated using an ”inflection point pri emperica number the cons Data hav the horiz number Still an 3 out a fine turbulent isotherm. COmputai The inabil theoretica accuretely layer of fit 1'2-3 Cell‘ Like the R, point principle of generalized excess entropy production," based on the emperically observed linearity of the generalized functional with the Rayleigh number at stable states, in order to carry out the numerical calculations within the constraints of the computers available in 1980. Data have been published to describe the heat transfer behavior of the fluid in the horizontal layer, typically reduced to a Nusselt number versus Rayleigh number correlation. Still an area of active research, Tian and Karayiannis (2000) recently carried out a finely detailed experimental study of two-dimensional, low level turbulence, natural convection in an air-filled vertical square cavity driven by isothermal vertical walls at different temperatures, to be used for computational code validation. The inability to accurately prescribe the wavenumber of convection in theoretical flow models has been suggested as one major difficulty in accurately modeling the behavior of the natural convection in a horizontal layer of fluid. 1.2.3 Cellular Convection in Porous Layers Like the Rayleigh-Bénard problem in the fluid layer, cellular natural convection in a fluid saturated porous medium has also been a topic of much research 10 fortht fiukllz Beck( rnedka vannng predic: stabilit thatva contini determ across existen the mo. only wl Confirnt ”Mshah OihErs. JonSSOr Spheregl Ravieigh damfeh for the last 40 years or so. Many of the same issues surrounding the study of fluid layers are material for horizontal layers of porous media. Beck (1972) published a chart showing the flow planforms for flow in porous media-filled enclosures heated from below, based on Darcy law theory at varying aspect ratios. Beck showed that the energy theory for stability predicts a critical Rayleigh number, Raw of 41:2 or greater, and the linear stability theory predicts 3 R3,, of 4a2 or less, thus bounding and confirming that value as necessary and sufficient for convective motion to begin and continue for an enclosure or layer with solid, impermeable walls. Beck also determined the stable mode of cellular convection: that the number of cells across the length of the box and the direction of two-dimensional rolls or the existence of three-dimensional cells, determining that the cells take shape in the most square cross sections possible, and that two-dimensional rolls form only when the box is ”thin" in the direction of the roll axis. Beck's results confirm the results of Lapwood (1948) for the infinite horizontal layer, while illustrating flaws in the use of the nonlinear momentum term by Lapwood and others. Jonsson and Carton (1987) also obtained wavenumber data for water in steel spheres. The results showed two distinct wavenumber possibilities at each Rayleigh number, though both produced the same Nusselt number. Most data fell outside the Strauss stability envelope. It was suggested that the 11 Prandtl quesfic toroida compa wide dI Beck’s enCIOSL shapes differer critical Several Uhderc Saatjiar case of pr ESSurI demons nUmber conilam and SDEI are all 31 number, Prandtl number effect should shift this envelope. However, It is also questionable whether the wavenumbers in the supposed two-dimensional toroidal convective rolls in their experimental test cylinder would be directly comparable with long cylindrical rolls in rectangularly bounded or infinitely wide domains. Beck's analysis has recently been extended by Wang (1999) for the case of an enclosure with isothermal top and constant flux bottom. The pattern of mode shapes are similar to the constant temperature bottom wall case, but with different sizes: there are generally fewer rolls for a given box width, and the critical Rayleigh numbers are generally much lower. Several researchers have shown that the critical Rayleigh number can shift under certain circumstances. Nield (1982), correcting an earlier work of Saatjian (1980), derived the proper (non-Boussinesq) equations for the special case of flow of an ideal gas through a porous medium, accounting the pressure work done during volume changes in the fluid. This work demonstrates conditions which shift the critical (porous media) Rayleigh number above the usual value of 41:2. Other non-Boussinesq effects are contrasted with the behavior of a Boussinesq liquid: the thermal expansivity and specific heat decrease and viscosity increases with temperature, which are all stabilizing effects for the ideal gas, also increasing the critical Rayleigh number. 12 Rees (1997) showed that the critical Rayleigh number for onset of cellular convection could also be increased above 411:2 when inertial effects, modeled by the Forchheimer extension, are made more important by a pressure-driven horizontal flow of the fluid through the layer. It was also suggested that this effect could be used to experimentally evaluate the relative importance of the Forchheimer extension under various flow conditions, by determining the relative shift in the critical Rayleigh number. Studies similar to those of the fluid layer have been performed to attempt to determine the wavenumber of convection for the porous layer. Somerton, McDonough, and Catton (1982) used the mixed finite difference- Galerkin method to solve the Darcy momentum and energy equations for the horizontal layer heated from below. To specify the wavenumber a stability functional based on the Glansdorff and Prigogine (1971) theory of nonequilibrium thermodynamics was used. It was found that at the preferred wavenumbers this functional is linear with Rayleigh number, with good correspondence to existing experimental data. There has been much greater disparity in studies of the heat transfer characteristic for the horizontal porous layer, including experimental data, theoretical, and numerical predictions. 13 Sc lve Br dil llll (i! th« 8C (1S na co the 9H Ra ant Ulll Ah the KEy Some early numerical flow and heat transfer resutls were obatined by Chan, lvey, and Barry (1970) who solved the two-dimensional steady Darcy- Brinkman equations with the Boussinesq-Oberbeck approximation using finite differences for flow in a porous rectangular enclosure. A possible angle of inclination was also included in the equations. Palm, Weber, and Kvernvold (1972) used a power series expansion of the Rayleigh number to determine the Nusselt number dependence for porous Rayleigh-Bénard convection according to the Darcy flow model. A study of the effect of the aspect ratio was conducted by Prasad and Kulacki (1984), who used a Darcy’s Law numerical simulation of two-dimensional natural convection cells in a wide and shallow rectangular enclosure with the temperature gradient along the vertical walls. They found that multiple convection cells occurred when the aspect ratio (depth to width) was less than 0.9, but that only a single cell was established when the aspect ratio was greater than unity. A further finding was that the various Nusselt number- Rayleigh number curves could cross, depending on both the Rayleigh number and the aspect ratio, and that the Nusselt number is not always maximized under the established flow pattern. An approach used to try and reconcile the divergence between numerical- theoretical and experimental data was put forth by Prasad, Kulacki, and Keyhani (1985). They proposed a model using an ”effective thermal 14 condUCi This effE conduct the path when m thermal empirica flow reg predictiv effective Kladias a Combinai cavity he isolating ' Forchheir f“HI/Chara experime, results usi energl’Eq Manmjc conductivity," which would depend on the flow regime or Grashoff number. This effective conductivity would tend toward the stagnant thermal conductivity at low flow rates, when conduction through the solid matrix is the path for more of the heat transfer. Under higher flow rates, however, when most of the thermal energy is carried by fluid conduction, the effective thermal conductivity would approach that of the pure fluid phase. This semi- empirical approach succeeded in collapsing much of the data in the Darcian flow regime, but requires iteration between the experimental results and a predictive model to analyze the results and determine the proper value for the effective thermal conductivity model. Kladias and Prasad (1991) performed experiments using 16 different combinations of fluid, solid materials, and particle sizes, in a rectangular cavity heated from below. This variety allowed comparisons of trends by isolating individual parameters: Rayleigh number, Darcy number, Forchheimer number, and thermal conductivity ratio, in order to ”isolate and fully characterize the contribution of the solid matrix properties." The experimental data were compared to the authors' theoretical-numerical results using a Darcy-Brinkman-Forchheimer (DBF) flow model and single energy equation, with investigations including effects of varying porosity and thermal conductivity caused by flow channeling near the walls and another trial including dispersive effective thermal conductivity. 15 A clear Darcy number dependence was shown by Kladias and Prasad. As the particle size was decreased, increasing the packing density and thereby decreasing the permeability and the Darcy number, the Nusselt number was decreased for the same Rayleigh number. A Prandtl number dependence was also demonstrated using various fluids with the same solid materials and particle sizes. Heat transfer increased with increasing (fluid-)Prandtl number. This increase was greater with larger particle sizes, because of the reduced contribution of the Forchhiemer inertial effects. Some increase in heat transfer was attributable to the accompanying increase in thermal conductivity ratio, which also explains the stronger Prandtl number effect observed in porous media compared with convection in a pure fluid. Kladias and Prasad emphasized the use of the Nusselt number and Rayleigh number based on the fluid properties, rather than the commonly used ”effective" properties of the fluid-medium combination, in order to more clearly show the divergence in the heat transfer results, and to illustrate the dependence upon the porous matrix structure and its thermal properties. Kladias and Prasad further reported their numerical results for comparison. The authors found that multiple steady-state solutions could be achieved, depending in the initial conditions supplied. To resolve their uncertainty in determining the proper wavenumber (or number of convective rolls), sinusoidal initial conditions were supplied based on the experimental 16 measure ‘zeroin' fimDB theexp expenn porosh‘ chosen doseri Mghco Uansfei undera dISpers hansfei 1.2.4 5 Thesec Wavenu SYSIEm j measurements of temperature variations along the bottom and top walls, with ”zero initial conditions in the interior of the cavity.” The DBF solutions (without wall channeling or dispersion) compared well with the experiments only for low Da and Ra,, and were over 20% lower than experimental Nusselt numbers at higher values. Wall channeling effects on porosity and thermal conductivity were modeled with empirical correlations chosen to match near-wall velocity measurements. Their results were much closer in agreement with the experimental data, except for the results with high conductivity solid matrix, which considerably under-predicted the heat transfer. All theoretical DBF predictions with and without wall channeling under-predicted the Nusselt number to varying degrees. The studies using a dispersive thermal conductivity model seriously over-predicted the heat transfer. 1.2.4 Entropy Generation Theories The second law of thermodynamics holds the most promise for closure to the wavenumber prediction problem. Glansdorff and Prigogine (1971) extended the concept of a measure of entropy to a ”generalized excess entropy,” which captures the motion of the system in a description of the state of the non-equilibrium system undergoing convective motion. The need for a wavenumber prediction to obtain closure 17 of the governing equations for the Rayleigh-Bénard problem was addressed by considering the non-equilibrium thermodynamic theory of Glansdorff and Prigogine, who defined a stability functional, termed the ”generalized excess entropy production,” which could be used to predict the stable stationary state of the system. This theory was used by Roberts (1966) to numerically study the planform shape and behavior of Rayleigh-Bénard convection in a pure fluid layer. McDonough (1980) also followed the analysis of the theory of Glansdorff and Prigogine for the fluid layer, but determined that the formal calculation to minimize the stability functional was too computationally expensive. Instead, a semi-empirical method was employed by determining that the functional should be linear with the Rayleigh number when evaluated at the preferred wavenumbers. This did result in wavenumbers which decrease with increasing Rayleigh number, which was consistent with experimental results and had not been achieved by previous numerical-theoretical schemes. Somerton, McDonough, and Catton (1982) used McDonough’s approximation of the Glansdorff and Progogine theory to predict the wavenumber for Rayleigh-Bénard convection in the porous layer using a Darcy flow model. The dissertation of Somerton (1982) also followed this wavenumber prediction scheme for the two cases of porous layers heated from below and those internally heated. Somerton also included the Brinkman term in the 18 momen bmmm the por terms iI Bejan ( for non led to t tool to It is fun same fa system Several Convect Demirel heattrai rate lhro MOSI rec from nati Theinco momentum equations to properly model the no-slip condition at the top and bottom horizontal walls. Georgiadis and Catton (1986) continued the work in the porous layer to include both the Brinkman and Forchheimer momentum terms in the numerical solutions. Bejan (1982) demonstrated another calculation of the entropy generation rate for non-equilibrium systems, including fluids in viscous flows. That work has led to the growing field of entropy generation minimization (EGM), a design tool to optimize thermo-fluid systems by minimizing the amount of lost work. It is further shown how natural systems tend to be ”self-optimizing" in the same fashion, leading to a ”constructal theory" for self-organization in natural systems. Several recent studies have studied the entropy generation rate in natural convection in porous media and in similar systems. Demirel, et. al. (1997) analyzed the entropy generation rate for convection heat transfer in a “packed duct," including maps of the entropy generation rate throughout the duct. Most recently, Baytas (2000) has analyzed the entropy generation resulting from natural convection in a porous cavity tilted to various inclination angles. The incompressible Darcy-Boussinesq equations (neglecting inertial and other 19 nonlinl genera and vis the cav minimi. 1.3 5 The wo numeric Darcy-B thermal behavio wide po l- c. w th w; 2' Ex De 3' E) COl nonlinear effects) were solved numerically for a square cavity. Entropy generation results, calculated by Bejan’s method to include both heat transfer and viscous contributions, show the distribution of entropy generation within the cavity. These results are used to determine the ”optimal” angle to minimize the entropy generation. 1.3 Scope of Dissertation The work to be presented in this dissertation is the development of a numerical simulation of Rayleigh-Bénard convection in a porous box, using a Darcy-Brinkman-Forchheimer momentum model and single porous medium thermal energy equation. When the box is made sufficiently wide, the behavior in the central region is assumed to be the same as in an infinitely wide porous layer. Results from this simulation model will be used to: 1. Compare the wavenumber of convection and resulting heat transfer with similar results of a simulation of the porous layer which uses the thermodynamic entropy generation theories to predict the wavenumber. 2. Examine the effect of the Brinkman and Forchheimer extensions to Darcy’s law on the heat transfer and wavenumber. 3. Examine the effect of the box width on the critical Rayleigh number for convection in the porous enclosure. 20 This has hia The devn con: gen: func Rnh The fort exce thet resu ShnL conv entrc lhod DGRO Iehfih and ti deScr 1.4 Outline of the Dissertation This dissertation is divided into five main chapters. This introductory chapter has served to introduce the problem under study (cellular natural convection in a layer of porous media), and expose the body of work which precedes it. The second chapter will document the assumptions and lay out the development of the equations used in the numerical solution: the governing conservation equations and associated boundary conditions, the entropy generation rate equation, and the excess entropy production rate stability functional. The third chapter details the implementation of the mixed Galerkin- Finite Difference method used to numerically solve the governing equations. The post-processing computations, such as calculation of the stream function for flow visualization and calculations for the entropy generation rate and excess entropy production rate stability functional, are developed, as well. In the fourth chapter the flow, temperature, heat transfer, and wavenumber results are presented, and used to evaluate results from the porous layer simulation, which are based on alternate methods used to predict the convective wavenumber. These methods include the maximization of the entropy generation rate, and the minimization of the excess entropy production rate stability functional. A final chapter summarizes the work performed, the conclusions reached, and lists directions for further research relating to this topic. Details of inner products used in the Galerkin equations and the computation of the integral functions using the Galerkin variables are described in several appendices. 21 In this i introdu their WI consen describ for entr which a predicti 2.1 A The aSSl here for follows; I. Th mc be ll. Tw dir lso .l Chapter 2: Governing Equations In this chapter the governing equations for the porous enclosure will be introduced in basic dimensional form, non-dimensionalized, and reduced to their working forms for the numerical solution. These will include the conservation equations of mass, momentum, and thermal energy, which describe the heat and fluid flow in the porous enclosure, and the equations for entropy generation and for the generalized excess entropy production, which are the measures used for the two thermodynamic wavenumber prediction theories for cellular convection in the porous layer. 2.1 Assumptions The assumptions made in the derivation of the working equations are listed here for completeness, and to serve as a framework for the analysis which follows: I. The steady state condition will be assumed in the solution of the momentum and energy equations. The time-dependent equations must be considered in the analysis of entropy production rates, however. ll. Two dimensional fluid and heat flows are assumed in the y- and 2- directions. lll. (Boussinesq-) Incompressible flow is assumed lV. Isotropic materials with uniform properties (except Boussinesq approximation for fluid density) 22 VI. VII. VIII. XI. XII. XIII. XIV. FI Lc be en effi 88 ea I VI. VII. VIII. XI. Xll. XIII. XIV. (Modified) Darcy-Brinkman-Forchheimer (DBF) momentum equation models the flow in a porous medium Boussinesq fluid approximation: p=p0[1-fl(T*-TO*)] (2.1) with reference temperature at cold wall temperature ( T0‘=TC‘ ), Velocity represented by Darcian (or filtration) velocity integrated over some local neighborhood: to: _. 1 _ e u = 2 Inf dA (2.2) No-slip and impenetrable flow conditions at solid boundaries Gravity acting in .2 direction ( g = — g1?) No phase changes (de = c dT) No internal heat generation in either the fluid or the porous matrix Stationary solid phase of the porous medium Fourier Law for heat conduction applies at interstitial level Local Thermodynamic Equilibrium (LTE), which takes two forms. One is between solid and fluid phases, implied by a single continuum-model energy equation using a local volume average temperature and an effective thermal conductivity of the fluid-solid medium pair. The second is in the specification of a local entropy state value under non- equilibrium conditions (see Sections 2.5.1 and 2.6) 23 t8 XVI. N- bl 2.2 G A diagra Figure 2. (position superscr' values w. througho z'=o Flsure 2.1 XV. Infinitely conductive top and bottom plates at constant, uniform temperatures TC' and TH'. XVI. Neglect viscous dissipation of mechanical energy into thermal energy, but consider for second law analysis 2.2 Governing Conservation Equations A diagram of the Rayleigh-Bénard system in a porous enclosure is shown in Figure 2.1. Note that the dimensional quantities for the state variables (position, velocity, pressure, temperature, and time) will be starred (with a superscript asterisk) in this section, to distinguish from the dimensionless values which will be developed in the subsequent section and will be used throughout most of the remainder of this work. =2” v',w*=o 3 z =L ,,, //j l o I ‘._.,//: .x/ ‘ l/ // ‘ 3T //, -,"_/ \ /, 6T :0 // " Ni ‘tln ‘ I HI, "I \I .=0 ay' r“ l’r’ I I n': t I '7, _ , j ‘/ = a ell \‘ // ‘ /' ',/. / \ . , //. 5"” I” //// ' ; //, zt=0 /:// / l T'=TH' v,w’=0 . ‘ y =0 y =W Figure 2.1 Diagram of Physical System 24 2.2. The GUI” wit vel the 2.: Th ea. vis BX 2.2.1 Mass Conservation (Continuity) The conservation of mass equation for the porous system, beginning with compressible and unsteady considerations, would read :6 + v (pr-2*): o (2.3) with the notation {4* = u*z° + v*]' + WT]; for the volume-averaged Darcian velocity vector. By then assuming incompressible fluid flow (in the sense of the Boussinesq assumption) the velocity becomes divergence-free. .. an" at" aw" - -O (2.4) V-z‘i— *+ *+ *_ a} dz 6:? 2.2.2 Momentum Conservation This analysis uses a modified Darcy-Brinkman-Forchheimer momentum equation. The form chosen includes the time dependent term, an effective viscosity for the porous medium in the Brinkman term, and a parametric exponent in the Forchheimer term: * l at e __ 2_:I- i_* B_* (n—l) *— -£pa*-V(: +E§+MV u — Ku —£nu uJ—O (2.5) (A) ) U (D) (E) (F) 25 Note ti operati Thougf has a u. momen (A) (B) IE) IF) T P sh De bo DIE Mc No ex; Note that the application of the Laplacian operator on a vector implies the operation of the usual (scalar) Laplacian component-wise: V2‘*- 072; é’zuI' é’zu‘l ,, u = *2 + *2 'I' *2 1 dc é) a} 52v!“ 022$ é‘zv* A. + &*2 + é)*2 + &*2 .] azw" azw“ azw“ . + ,2 + ,2 + ,2 k at Q) a.- (2.6) Though much more complex than the linear Darcy law, each additional term has a unique contribution in describing the physics of the flow. The respective momentum terms represent: (A) (B) (C) (D) (E) (F) Transient term — time rate of change of local fluid momentum. Pressure gradient. Hydrostatic gravitational body force. Brinkman term — fluid-fluid shear. Also considered a "boundary" term, since it accomplishes the effect of no-slip at a solid boundary. Darcy term - viscous shear at the (interstitial) fluid-solid phase boundaries. Also considered a diffusion law for flow caused by the pressure potential. Modified Forchheimer term — describes convective fluid inertial effects. Note that the coefficient m must take units which balance the modified exponentn. 26 2.2.3 Conservation of Thermal Energy For a Boussinesq-incompressible fluid in a porous medium without phase changes or internal heat generation, a differential energy balance produces the thermal energy energy equation (Bejan 1984, Cth): fl _* III 2 III _ (pc)m a, +(fx)fuv -VTJ- kmvl T + fife _ o (2.7) The terms in this equation represent physical processes of: (G) Local (sensible) thermal energy storage in the porous medium, considering both the fluid and the solid matrix. The effective thermal storage capacity, (pc),,,, is typically a local volume-weighted average of the properties of the solid and the fluid. (H) Net convective transport of thermal energy by the fluid flowing through the medium. (I) Net conduction. heat flux, modeled by the Fourier law and characterized by the effective km for the fluid-porous media system. The value of km, though bounded by several models, must typically be determined by experiment. (J) Internal heating (dissipation) produced by work of viscous flow. This term will be neglected for the solution of the flow and temperature fields, but must be considered for the analysis of the entropy 27 gl 2.2.4 B Implemi side wal W The tern and the I Addition; back wall 2'3 No The State DIOCBSSE generation. The calculation of this term will be revisited in Section 2.5.1. 2.2.4 Boundary Conditions Implementing the no-slip and impenetrable flow conditions at the top and side walls: 11* :I- = 0 , 11* at: = 0 , 11* at: '-' 0 , * = 0 g :0 Z = L y =0 y*:W III II! it! =0: : : O : 0 : V zit—O O , V 2*:L , V y*=0 , y*=W 0 (28) w*.=o,w*.=o,w*.=o "‘ =0 :: =0 " =L y =0 y*=W The temperatures are specified on the top and bottom plates, II: it: * it: T ZT=O:TH ; T Z*=L:TC (29) and the right and left lateral walls are assumed adiabatic ar‘ 07'" = , = 0 (2 10) t III . é) y*_0 é; y*=W Additional no-slip and impenetrable, and adiabatic conditions on the front and back walls would also be assumed for a 3-D formulation. 2.3 Non-Dimensionalization and Scaling The state variables are scaled in the parameters typical for convective processes. This includes offsets for the pressure and temperature by 28 subtracting the hydrostatic pressure gradient and the reference cold temperature: * II t x )2 z x E —" ; E _ ; Z 5 — L y L L t*a _ m t: L2 PE K (P + pogz ) #f am T“ - TC“ T" - TC" T5 :i- :I- " at: TH - TC AT The length and time scales imply the proper scaling for the non-dimensionalized velocities: * L u: z—u ; v: am MR» (2.11) (2.12) (2.13) (2.14) (2.15) (2.16) For convenience when substituting, also note the reverse of the definitions * ya :I: per P =P[ fme-pogz =P[ me]_p0ng T" = TAT" + TC“ 29 (2.17) (2.18) and the length 5 Now we equatio 2.3.1 IV SubstitL produce which re. and the non-dimensionalization of the gradient operator corresponding to the length scaling t l 1.5 .3 Afl V =—V=— '— '— k— 2.19 L L(léc+Jéz+ a] ( ) Now we apply the scaling of the variables to each of the conservation equations and boundary conditions. 2.3.1 Mass Conservation Substituting the nondimensional variables into the continuity equation produces e921. Lil, Late. 0 2 20 0‘,"ch amézL am&L- (') which reduces to V ‘-O-é—u-+é)—+—O:w— 221 u“ “at d) a (' I 2.3.2 Conservation of Momentum Next substitute the nondimensional variables into the momentum equation (2.5). Considering term-by-term : (A) (employing the Boussinesq approximation): 1 072‘ 1 (am)(am) 0‘12 [poamzjgli E grim—Z- 73:?01 30 (B) OI T (C) (2 ,1 (D) I I (E) - (F) n f i? Re-asser KL \I lflfam/I I"I l‘ 5,1 (B) on the pressure term 1 ”fan: VP" = —v P — L L [ K ] pog(z )] pfam A = VP- k i .1, j a. (C) (again employing the Boussinesq approximation): P3 = Po[1’ 47* _ To*)](- 31;) -p0g[l— fl(TAT* + T5“ - T0*)]IE —p0g[1- 5(TAT“)]IE 2 (D) #mv *zl-‘T : flffifl(l) VZEEflfi) = [PA](ELa—'"]V2fi (E) ———a* =[ (F) neglecting density variations (by Boussinesq again), p (II-1):: p0 (am)-("'l)(am)_ — u —' 7:“ u —' u m - m L fl(am)nlfii("- l)2}. _* m T Re-assembling the momentum equation and multiplying every term by KL [ ] yields: .uf am L AT" . — @% g-VP+(K pog'g ]Tk+[-&"-](—%)V2ii-a EflfL .uf am .uf L (2 22) _ ._"_"’_Q_[E.m.)n HUI-”Fr, #f'" L ( 31 Now collecting terms in the following dimensionless groups: Rayleigh number AT* L3 Ra a pogfl (2.23) qu am Darc number V Da 2 :2 (2.24) L porous media Rayleigh number p0 g [M T*x L Ram E Ra-Da E (2.25) tuf am porous media Prandtl number vf ,uf Prm a 5 (2.26) am pOam porous media viscosity ratio ,2: 0m 5 — (2.27) I“f and a modified Forchheimer number (n- 2) Fsm a £(EE—j (2.28) m yields the non-dimensional momentum equation: D d? . F JD .. ——l—-VP+RamTk+DaamV2fi—fi——€m——-q-Iiil(n ”a: o (2.29) Pr," 5t Prm 2.3.3 Energy Equation Again, substituting the non-dimensional forms into the energy equation, neglecting the viscous dissipation, considering termwise : (em-:5 = (pal-E31) : (T*AT* . TC.) (GI [(pc)mamAT* * 5f * 5i _. L2 32 Re-assembling and dividing through by L2 L2 defining the thermal storage capacity ratio (AC) (2,, a ——’-"— (2.30) (wa yields the dimensionless energy equation: 07 Qm-d—Jrii-VT-VZTzo (2.31) 2.3.4 Boundary Conditions The nondimensional variables exhibit the homogeneous no-slip and impenetrable boundary conditions for the velocity, expressed as: ulzzo = 0 ; ulz=1= O ; uly:0 = O ; ulyz}, = 0 vlz=o = 0 : vlz=1= : v|y=o = 0 ; VIM = 0 (2.32) wlz=0 = O ; w|z=1= 0 ; W|y=0 = 0 ; wlyz}, = 0 where the aspect ratio, y, has been defined W E _ 2.33 7 L l I 33 The temperature boundary conditions are also simplified to a" tfl NF0=1; flfll=0; —+ =0; -+ =0 03m 03).”0 0yy=7 An illustration of the porous box system, as described in the non-dimensional variables, appears in Figure 2.2. \\C \ . L\‘ \ : I f ,1 r.“ ‘\ . \ >1 / ’/ ‘\é- \ u ::::t I . ',¢ . \\ ' \ \ _ q I f” , Q ~ \ \\) \ I I I I ’7 \ \\ \\\ \ \ \ T‘ ; " -‘_-_,r \ \ \ \ \\\ ‘Q é II o \b. \\> \ \ \ \ /" I Q \ \\: . . l : I \c \\‘ .. ,’ , .r I / .\ \ ‘ \ \ \\" . /1 \“ I) ‘ C \‘ \ \ \Kf\\ S E II o z=0-—z’ l l T=l uw=0 l g y=0 y=7 Figure 2.2 Non-Dimensionalized Porous Box Heated From Below 2.4 Reduction of the Mass and Momentum Equations Taking the non-dimensionalized continuity and momentum equations, (2.21) and (2.29), and thermal energy equation (2.31), and assuming steady, tvvo- dimensional flow and heat transfer reduces to the four scalar equations: (con Inurty) é! & - ( - I W 0’21) 52V _n_] (y-momentum) - —+ B 7+ —2- - v - Mlul V = 0 (236) e» a} a: 23.3 dual): _ _ We) - RT (237, (z-momentum) 01' @2 012 w u w — — . ar arr azT é‘ZT (thermal energy) v— + w— - —— -— = (2.38) a} 52 azfaz making use of the definitions: K rum BE DaO' = —— m L2 W (2.39) F x/D Ma Sm " (2.40) Prm R 5 Ram (2.41) (17(5 v2+w2 (242) 35 Next cross-differentiate the y-momentum equation with respect to z and the z- momentum equation with respect to y to get: 322—5 83v 0'0"] fl éoeln-l) —@&+ @2&+&3 — u v 32' dz 0 521) B[a3w 33w] av a(|‘l"“ ) er a a“ w a —%+ @3+@&2 -——"— M— =-R—" (2.43) (2.44) Subtracting the terms of equation (2.43) from (2.44) eliminates the pressure and produces the single combined momentum equation: B{- 53v _é‘3v+0"3w+ §3wj+fl_§v_ @201, 01,3 @3 @&2 d: é, aim) The horizontal velocity can be eliminated from the linear terms of the equation by differentiating again with respect to y: B[— 54v _ 54v +§4w+ 34w ]+ é’2v_§2w @301, é’&3 @4 @2322 @3‘, @2 (92 _ a2 n. + final" 11)-:6701-2) 1w) azr =-R——2- and using the following differentiated forms of the continuity equation: 54". a4.) _O @3& @2é2- 54v 54w é»&3+§z4:0 0‘21} 52W_O afar 36 (2.45) (2.46) (2.47) (2.48) (2.49) Substituting then yields 33T+2@20~22+&4 "@2'022 0’2 (lfiln-lv)_ 5’2 (lfiln—lw)] : Kill" + (2.50) 0W 37 _ a2 To further simplify, lump the nonlinear Forchheimer terms with the definition a2 _, ,,_ a2 g n_ Q(v,w) a [352—014 lv)— g-Z-(Iu) lw)] (2.51) for the resulting form of the combined momentum equation: a4 a4 a4 a2 a2 azr B;+2@2&2+&4 --5)7-;7w+ MQ=-R—7 §.(252) which is now a linear equation for w in terms of Q and T. 2.5 Entropy Generation Rate Equation Since it is an irreversible heat transfer process, cellular natural convection will continue to generate entropy and destroy available work. The processes through which entropy is generated are both the transfer of heat across a finite temperature gradient, and the viscous flow of fluid through the medium. Although the viscous dissipation term is often neglected in the thermal analysis for fluid and porous media flows, it must certainly be included in the Second Law analysis of such flow systems. This is especially so for porous 37 media fl pores 0 2.5.1 E The ana set fortI D-B-F ni Using a) porous I to expre medium media flows, in which the viscous fluid flow is largely diffusive through the pores of the medium. 2.5.1 Entropy Generation Rate Calculation for Porous Media Flow The analysis of the entropy generation rate, Sim. will follow the development set forth by Bejan (1982, Ch. 5, and 1984, Ch. 10), with adaptations for the D-B-F nonlinear porous flow model used for this work. Using an entropy flow balance of a differential element of the homogeneous porous media model, the Second Law of Thermodynamics can be employed to express the volumetric entropy generation rate at a point in the porous medium as S gen: *V-q-Tq-VT +p I: (2.53) where (j is the heat flux vector, and s is the local entropy state. The first and second terms on the right-hand side account the net rate of entropy transfer caused by convective and conductive heat transfer at the point, and the last term measures the net accumulation of entropy at the location. Conservation of mass has been assumed, balancing the in- and out-flows at the point. Again, the starred quantities indicate the dimensional values of state variables. 38 (f) 8f A note should be made here about the use of the thermodynamic state of ”entropy" throughout a system undergoing natural convection. Such a system is obviously not at equilibrium, as required by definitions of classical thermodynamic properties. The requisite assumption is that of local thermodynamic equilibrium (LTE). (Note that this assumption is also used, with a different implication in the use of the single thermal energy equation in Section 2.2.3.) The LTE assumption simply requires that ”the local entropy, s, is the same function of the local macroscopic variables as at equilibrium state." (Glansdorff and Prigogine 1971, p. 14) Armed with the LTE assumption, the local entropy state can be related to the local state variables in the solution: temperature and velocity. First, express a First Law balance equation written at the point in terms of the internal energy, e, and in the absence of internal heat generation: De 4: ‘4: met.=-V-ci-P (V-u )+rzf¢ (2.54) where (I) is the viscous dissipation function for the porous flow. Now employing the material derivative of the Gibbs relation: Ds p De P" me ,0 “— = — e - (255) "’ Dr" T” Dt pm T“ Dt" and the mass continuity statement at: TPmV'u :0 (2-56) 39 10 COm The he To exp therme work (I Brinkm So the i to combine equations (2.53), (2.55), and (2.56) as 4: 1 _ at: #f S gen=-;2—q-VT -I- T*(D (2.57) The heat flux can be expressed using Fourier’s law, a = — ka r“ (2.58) To express the viscous dissipation function, (I) , we note that the rate of thermal energy generated by viscous dissipation must be equal to the rate of work done overcoming viscous forces, as expressed by the Darcy and Brinkman terms of the steady state D-B-F force-momentum balance. :I: I“ at: :I- )1ch = a (7ft? - any 2a ) (2.59) So the viscous dissipation can be written *2 * III =—n "’1‘ -V2fi K ,Uf e t. t. .i (2.60) 1 *2 *2 #m 452v 52V *Oflzw 52W z—V +w ‘—" 42+ *2 “V *2+ .2 K H d) a é) a- The expression for the entropy generation rate under D-B-F flow becomes 2 l . km at S gen - T‘Z VT - amt-2 .vza") (2.61) 40 2.5.2 Non-Dimensionalization Substituting the definitions of the dimensionless state variables yields 2 . k Afam 1 s = ———’" IV 712 + [Ii-2P — B ii-V 2a l (2.62) gen L2(T+ 192 [HLZATT] (T+ T) ( ) Or, expressing in the two-dimensional form, the local entropy generation rate can be calculated S... _ ___]£r1_(fl]2+(§f_)2 ge" — L2(T+ 2')2 a} 01' ,ufamfl 1 52v é‘zv 62w é‘zw + [KILZAT*)(T+ T){VZ + w2 - Blvi-é7+ g) + w[?+ Ell where the absolute temperature has been non-dimensionalized with the (2.63) parameter, 1:, as * T r: C... AT (2.64) Factoring the porous media Rayleigh number out of the coefficient of the viscous flow term ”fam2]_[ ”ram Mammal?) KLZAT“ pogflxLA T" L 1 (2.65) = (Ramil “Miogflj 41 suggests a scale factor for entropy generation rate in the remainder. Dividing through by (W) yields the non-dimensional equation j: v [ 2 — (72—) 2.66 ng HBlT+ {)2 l T|2+ Ram (T+ r) lul (u u) ( ) where the non-dimensional entropy generation rate is defined with the scale factor L . S , [—ls (2.67) g " armada/3 g8" and the buoyant heat transfer number on the conductive term is defined as km W (268’ H35 The total rate of entropy generation, G, (per unit depth) is then determined by integrating S over the domain. gen 1 7 G = I jsgen dydz (2.69) O 0 Also of interest is the irreversibility ratio: a ratio of the fluid flow irreversibility to the heat transfer irreversibility (Bejan, 1984) ”faszC* ka(A Til )2 ¢ (2.70) 42 2.6 Generalized Excess Entropy Production Rate Functional The non-equilibrium thermodynamic theory of Glansdorff and Prigogine (1971) presents a thermodynamic potential which would be minimized for convective flow in the stationary state. The development here follows Glansdorff and Progogine, and also the dissertations of McDonough (1980) and of Somerton (1982), to develop a suitable expression for the generalized excess entropy production functional, ‘I'. To begin, define the "generalized" entropy as _.*Z p_f 1 pm 27;,“ (2.71) The starred quantities are, again, the dimensional quantities, which will subsequently be replaced with non-dimensional expressions. The T0' represents the local, unperturbed temperature. The first term represents the local specific entropy (per unit mass) of the averaged porous medium. The second term represents the kinetic, or hydrodynamic contribution from the fluid flow. Again, the assumption of local thermal equilibrium (LTE) must be employed. We shall next examine the perturbation about the steady condition, so that the perturbed state equals the steady state plus the fluctuations or variation, 43 e.g. T=T0‘+0T. Taking the second local variation (so 6F6y=62) of the defined functional yields the ”excess" quantity, given as p 1 2 azr = 523- —f—|a2* (2.72) pm TO since ii is an independent variable, 6227* a O, and To‘ is constant with respect to the variations. Now under the condition of “local thermal equilibrium" (LTE), then 62s < 0 everywhere locally, assuming a stable local equilibrium (Glansdorff and Progogine, 1971). Therefore azr < o (2.76) everywhere, as well. So we may also consider the second variation as a Lyapunov stability functional, with the following condition for the steady state I Til 2r) 2 o (2.77) To evaluate this expression, we again look to the Gibbs relation. de Pdv dS = —; + T (2.78) T T Taking the second variation, 6 2s , and noting that e & v are independent with respect to variations in 3, gives 523 = 5[}17]&+5[%J6v (2.79) Applying Euler's Theorem on homogeneous functions for the independent variations 5v and be in the time derivative gives the expression 1a 1 awe) P" 5(a) 3578294?) a" 5lth at" ”'8‘” which can be substituted into the stability criterion as 1a 1 ae P" av )0 1 a __.._(52r)= (4)4.z+5(.-..)t>-_f_._eiz 2a T a T a pm 2T0 a But, by the Boussinesq incompressibility assumptions, we can assume that 2 O (2.81) 6v=0, and if no phase changes occur we have that 0e=cm 67’. Also employing theidenfifies: Xi) = - -—1—aT"' (2.82) T To"2 and .. aw): 4...)...48 d c? a”? Substituting these identities reduces the stability criterion to * * —Ti;:§ar*-€l%l—%Flea*é:—lzo (2.84) 45 Next 5 and m Next n numbe Integral entrop) ................... In order: expressic moment ~— I and ener Next apply the non-dimensionalizations as before L __:I= I am am L AT ii and multiply by To'z, so — c,,,(AT* )2(a_m)67.§£)__ fl-(TOA T“ + TC*)[9’"—3]a2. 5W) 2 0 (2.85) L2 a pm L4 o? 2 . . pm Add . . . Next multiplying by 2 .. and employing the dimenSIonless pf am flfA T numbers introduced in previous sections, yields the result 01%) at 51“ D - QmHBRama'I‘d )- a (TO + r)&- 61‘ Prm 2 0 (2.86) Integrating over the domain { y 5 (0,7), 2 6 (0,1)} results in the global excess entropy production rate functional, ‘1’, 17 51‘ D an 5 LI! 5 H -omHBRamaT§(5—)-—I—)r—a-(TO+ flan-é?) dydzz o ((2.87) 00 m - 5 In order to evaluate ‘I’, we must evaluate the time derivatives. These expressions come from the quasi-linearized first variations of the unsteady momentum equation: D & ~ E1257) = —V(éP) + Remark + IN 2(a) - m - M|a|("' "627 (2.88) m and energy balance equation: omiigT—L—avwrhvzwfl (2.89) 46 Substituting and expanding the differential operators and dot products in scalar form results in the expression IH R awn + 5(a) _ awry arm ‘ B am v @ W 01, 4’2 312 ’ asp) 5(a) ' -§v—- __ a we I? +Ram5w51' ‘I’ .=. (j ( 6v[52(5v)+32(&)] ) )dydzzo (2.90) - (To + r) 4,2 @2 + (”Vein 2(a)] \ a? 522 l _ — (1 - Mlt‘rl("' I>)(av2 + 6w2) +B From this expression, the value of ‘P can be computed from the solution of the temperature, flow, and pressure fields as determined from the mass, momentum, and energy conservation equations. 47 Chapter 3: Numerical Solution Method In this chapter the equations and routines are developed which enable the solution of the Rayleigh-Bénard convection system in the horizontal porous enclosure. The entropy measures are also computed, in order to facilitate comparisons with the results from the porous layer, and to evaluate the thermodynamic theories. The numerical method employed is a mixture of Galerkin and finite- differences. The Galerkin method with trigonometric/Fourier bases is used to solve for the horizontal dependence in the temperature and flow fields. The vertical variations are then discretized with a finite difference mesh. This method was developed by McDonough (1980) for the Rayleigh-Bénard convection problem in the pure fluid layer, and for the porous layer by Somerton, McDonough, and Catton (1982) and Georgiadis and Catton (1986). The method has been shown to be reasonably computationally efficient. It is an extension of these works that will be used to produce data with which to evaluate the thermodynamic wavenumber prediction theories .for the porous layer. Therefore the same numerical technique is chosen for this dissertation to solve in the porous box, in order to minimize differences which might arise from the numerical method. 48 3.1 G 3.1.1 G For the ( expansk safisfytl Where The pre Howev fUnctio 00ndni 3.1 Galerkin Forms of Governing Equations 3.1.1 Galerkin Ansatz For the Galerkin basis functions, the following trigonometric half-wave expansions for the velocities and temperature, have been chosen which satisfy the specified essential or natural boundary conditions. K v(y,z) = glelzl sin(aky) (3.1) K w(y,z)= Z Wk(z) sin(ak y) (3.2) k=l K T(y, z) = konk (Z) cos(aky) K (3.3) = T0(Z) + Z Tk(2) cos(aky) k=l where rrk (3 4) a E — . k 7 The pressure is not needed for the momentum and energy solutions. However, for post-calculations, is should be expanded with the cos a, y basis functions, to remain consistent with the velocity bases and boundary conditions: K P(y,z) = 2 Pk (z) cos(aky) (3.5) k: l 49 3.1.2 IllaS The Galerk substitutin dimensior evaluating Note that has been Next the over the ‘ Since m can be Then u. 3.1.2 Mass Continuity Equation The Galerkin form of the mass conservation equation is obtained by first substituting the expansions for v and w, equations (3.1) 8 (3.2), into the two- dimensional, incompressible continuity equation (2.35). After substituting and evaluating the partial derivatives termwise, the equation becomes K 2 (Vkak cosaky+ W]; sinaky)= 0 (3.6) k: 1 Note that the notation of functional dependence on 2 for the modal velocities has been suppressed. Next the inner product is formed with the pressure basis functions, cos a,,, y, over the y-dimension from 0 to y 7 K [1,210/ka cosaky+ W); sinaky)cosamydy= O (3.7) 0k : Since the V, and Wk' factors of each term are functions of 2 only, the integrals can be passed through the summation and evaluated termwise, as k}: K(Vk aklo cosakycosamyafy+ W); I; sinakycosamydy) = 0 (3.8) Then using the orthogonality of the cosines: 7 £7731 ifk=mandk¢0 (cos(aky) cos(amy)a§2 = (3.9) 0 O ifk¢mork=00rm=0 50 and defining Am as the inner product of sine and cosine over multiples of the half-wave interval: 7 Akm a (cos(aky)sin(amy)ajz 0 iy( l 1 ) if(k-m)isodd 7r k+m k-m andm¢0 (3.10) if k=m : t 0 or (k — m) is even 1— —1’" 1[——(——)—l if k=0 and m a o 71' m the continuity equation reduces to K a 431V“ ZAmkW); = 0 (3.11) k=l The Galerkin form of the continuity equation is used to recover the horizontal velocity terms, Vm , from the W,,, terms determined from the solution of the combined momentum equation. 51 3.1.3 Galerkin form of the Momentum Equation Into the combined momentum equation (2.52) the expansions for velocity and temperature are substituted, assuming the derivatives can be evaluated term- wise in the summations, which yields: sin(aky)+ MQ K z[B(ak4 Wk -2(1k2 Wk +ka')+ak2 Wk _Wk’ “1 (312) K = R2 ak2 Tk cos(aky) k=l Again, the notation of functional dependence on z for the modal velocities and temperatures has been suppressed. Now form the inner product of each side with the vertical velocity basis functions, sin am y, over the domain of y, (0, y): K {2 [3(ak4 Wk - 201:2 Wk" + Wk'm)+ ak2 Wk - Wk" ]sin(aky)+ MQ} k=l CS—HV K (3.13) R akz Tk cos(ak y) sin(am y)dy k: l - sin(am y) dy = O‘-—.V Carrying out the integrals termwise and using the orthogonality of the sines: 7 2’___.___’”" ifk=mandk¢0 . . 2 2am I s1n(a,, y) sm(a,,, y)dy = (3.14) 0 0 ifk¢mork=00rm=0 1 and the inner products of sine and cosine over multiples of the half-wave interval, equation (3.10), produces the simplified Galerkin form of the momentum equation: g-[B(am4 Wm — 2am2 Wm" + Wm"")+ am2 Wm - Wm + MQm K (3.15) where the inner product with the nonlinear term is defined: 7 7 Q... a lesinamydy= ggghalm‘v) -%(Ia|""w)]sinamydy (3.16) 3.1.4 Expansion of the Nonlinear Forchheimer Momentum Term To facilitate the linearization and computation of the nonlinear inner product term, Q," , we first identify the nonlinear component, defining n-I 775 lfiln-I : v2 + (4)2]? (3.17) The expression Qm is integrated by parts, with the result of 7 7 0., o - am lie} .-.,,,((g(.)-gtnj..w O (77w)] cosamydy I 3, elm . a a , Qm = [5070- gl'IWIISlnamy (3.18) since the sine multiplying the first bracketed term vanishes at the boundaries. 53 Next we integrate the second term by parts again to arrive at the form: 7 7’ 3( —am (507v )cosamydy + [amnwcosamyK + am2 lnwsinamydy 7° y 0 (3.19) -(r7v)(-0: )cosa ydy+ a 2 (nwsina ydy am &( m m 0 m 0 The (1) w) term in the brackets also vanishes at both boundaries because of the no-slip boundary conditions. Now we modally expand r) in terms of the velocity basis functions as: K 77(y,z)—) Z Hk(z)sinaky (3.20) k=l As will be shown, this will replace the need for quadrature integration of this nonlinear term with the relatively efficient Fourier sine series expansion in the solution procedure, and allow for improved linearization and convergence on the solution of the modal equations. We can now substitute this expansion (3.20) and the velocity expansions l3-1l 30d (32) into the reduced expression for Qm A 7 a K K Qm = —am IE (2 Hk sinaky][2 V, sina,y]]cosamydy k=1 1- 1 0 (3.21) 7 K K + am2 ([ Z Hk sinaky][2 W) sinaly] sinamydy o k=1 [:1 54 This expression may be simplified by combining the summations and differentiating the z-partial termwise 7K - ”1"“le KZI:; (HkV1)sinaky sina1y cosamy] dy 0k: II— I 7 K + am2 I 2 Z [(Hle)sinaky sinaly sinamyld)’ ok=ll=l (3.22) It is now apparent that the summations can be integrated termwise, with the modal coefficients, which depend only on 2, passed outside of the integral. This will produce the simplified form: m = Z Z [- am(HIch 'l' Hk V1,) Aklmfl) + amsz WI Aklm(0)] (3'23) : =1 after defining the two trigonometric triple inner products as 7 Aklmm) E (sinaky sinaly sinamy dy 0 (LXI 1 1 1 J if (k+l+ m) (3.24) _< 27: k- l+m k- l- m k+l+mJr k+l- m is odd T 0 if(k+l+m) L is even for k, I, and m = 1, 2, 3, and Aklmm a lsinaky sinaly cosamy dy 0 3.25 1 ifm=k-Iorm=l-k ( ) ifm=k+l for k, I, andm = 1,2,3, 4 Z 4 55 Thus, the expansion given by (3.20) has reduced need for quadrature integration of the Forchheimer term to simple evaluation the triple products (3.24) and (3.25). 3.1.5 Galerkin form of the Thermal Energy Equation The Galerkin form of the thermal energy equation is obtained in much the same fashion as the momentum equation. First the expansions (3.1) - (3.3) are substituted into the thermal energy equation (2.38). Again carrying out the partial derivatives termwise yields K’ K .K K [2 V1 sinaly][- Z ak Tk sinaky] + [Z Msinaly][Td + Z T); cosaky] k: [:1 k=l [=1 1 K K (3.26) + Z (11,2 Tk cosaky- T6’- 2 Té’cosaky = 0 k=I k-l Further, combining the summations simplifies the expression to K K F+k21(ak2 Tk— T")cosaky + T6: W1 sinaly- T6’ = 0 (3.27) k: 1:] where the double summations of the advective terms are collected as K K F a Z 2 [(- ak V1 Tk sinaky+ W) T); cosaky)sinaly] (3.28) k=U=l Next, for each equation for mode k, the inner product is formed with the temperature basis functions, cos any. 56 3.1.5.1 Mode Zero Equation The special case of m = 0 (so cos any =1), is considered first. This will produce the Galerkin equations for the mean-field temperature distribution, T0(z). The inner product formed is i’ K K I F+ Z (ak2 Tk — Tk" )cosaky+ TO' Z W, sinaly— T0" dy = 0 (3.29) 0 k= I I: l which reduces to . ,K 2 ,. F0+ToZa—W1-7To =0 1:1 I [odd (3.30) with the nonlinear advective product absorbed into the definition 7 a Odey 0 k: K 2 =1 K krr ) K 2k! = Z ( TVkaJr Z k2)WIT/t k=1 [:1 ak(12 - (Ir-l) odd l(- at V1 7k”: Sinakysinalydy + (W1 Till1i°°§akysmmydy (3.31) b 3.1.5.2 Higher Mode Equations Now we consider the inner products for the m=1, 2, 3 . . . mode cases: K 7 K ,lF +2 (alt2 Tk ‘ Tk")cosaky+ 702 W] Sinaly— 76'] cosamydy = O(3.32) 0 k=1 [=1 57 (Zn-la With the definition 7 Fm s j Fcosamydy (3.33) 0 the inner product of the cosines (3.9) and recalling the inner product of sine and cosine (3.10), the energy equation reduces to A 7! K E. (amZTm — T,;;)+ T6 2 Amle = 0 + —_ 2am I: l (l—m) odd for m=1,2,3, (3.34) To further simplify, the nonlinear double summation term can be reduced to 7 K K 13"," = [Z Z (- ak V, Tk sinakysinaly+ WIT]; cosakysinaly)cosamydy o k: 11:1 (3 35) K K ' = Z Z (' aszTk Aklm(l) + VVITIéA/c1m(2)) k: 11: l for m=1,2,3, by using the triple inner product definition of A,,,,,‘” from (3.25) and the definition 7 Au," (2) a Icosaky sinaly cosamy dy 0 7 (3.36) 1 1 1 1 ) (MN/”Hf k+l-m_ k—l+m- k—l—m if(k+l+m)isodd 0 if(k+ l+ m) is even for k, l, and m=1,2,3, 58 3.2 Modal Quasi-Linearization The fully non-linear form of the combined momentum equation, and therefore the coupled energy equation, cannot be solved directly in a single pass, and must be solved by iteration. The efficient scheme to be used is termed modal quasi-Iinearization. It is similar to the ”full" Newton quasi-linearization scheme, except that the equation corresponding to each Galerkin term, or mode, is linearized only about its own corresponding modal coefficient term. The remaining modal coefficients use the explicit value from the previous iteration. For the momentum equations, this requires that the W,,, term of Qm be extracted from the non-linear summation. For convenience, the terms Q," and ~ Am are defined so that Q... = Q... + 2me (3.37) with ~ K K K K Qm = -a,,, Z z[(H/2V1+ Hk Vz')Ak1m(l)]+ am2 2 Z [Hk W1 Aklmm) (3 38) k=ll=l k=l 1=l ' k¢m1¢m and K ( ) ~ = 2 0 Am - am kE [Hk Akmm ] (3,39) kodd 59 Now the Wm term in equation (3.38) is moved to the left-hand side of the Galerkin momentum equation, (3.15), while Q", is moved to the right-hand side, yielding the equation 2 ~ N "H % [Bam4 + am2 + ; 1m] Wm — (2am2 +1)Wm + BWm ~ K (3.40) = -MQm + R Z 01:2 Akak k=1 (k-m) odd where the definitions of Am in equation (3.10), and ak from equation (3.4), have been employed in the temperature term. In the energy equations, the mode zero equation contains no To terms, so the Iinearization is completed with a slight re-arrangement and relaxation of the velocity terms on the coefficient and higher-mode terms in [3‘0 on the right- hand side, with the result: K 2 . Z —W, Td— 7T": -F0 (3.41) 1= 60 In the higher-mode energy equations, the Tm and T3,, terms are factored out of the nonlinear term, Fm , and the remainder is defined in Fm as A ~ K K " amz VI Amlm(l)] Tm + [Z WlAmlmQ)] T62 (3'42) Fm = Fm + 1:1 [=1 where K K Fm E Z Z (‘ “k V1 Tk Aklmm + WIT/£41115”) (3.43) k=l l=l k¢ m]: m Re-arranging the energy equation with the mode m terms for temperature and its derivatives to the left-hand side and the remaining terms to the right produces the linearized equation K K (2:) amZ “((1012 VI Am1m(l)] Tm + [Z WIAmlm(2)J T61 ’ (1)7}; 2 1:1 1:1 2 K (3.44) =_Fm-Td ZAleVI 1:1 (l-m)odd The resulting finite difference equations, developed in the next sections, can then be solved for each mode successively. This de-coupling of the equation for each velocity and temperature term from each of the other modes dramatically reduces the storage and processing requirements in the matrix solutions. The coupling is maintained through repeated iterations, with a check for convergence when the equations are sufficiently well-satisfied. 61 3.3 Finite Difference Equations The Galerkin forms of the momentum and thermal energy equations are now reduced to K modal momentum equations and K +1 modal energy equations. However, the modal coefficients, Wm and Tm, are each still functions of z, and each equation represents an ordinary differential equation in 2. So next the z- variations of each modal coefficient are discretized into N grid points, and the corresponding derivatives are represented by second-order correct finite differences. 3.3.1 Momentum Equation The modal coefficients and their derivatives appearing in the quasi-linearized momentum equation, (3.40), are replaced with the following second-order correct finite difference forms: (Anderson, 1984) Wm(z) _) ij (3.45) n 1 Wm (z) —> 2%me — 2Wm,,- + Wm,,_1) (3.46) ”H 1 Wm (z) _) F(Wm,i+2 T 4KVmJH + 6Wm,i — 4KVmJ—l + Wm,i-2) (3-47) and the temperature coefficients Tm(z) —> TW- (3.48) and non-linear term expansion (which will be detailed in section 3.3.4) Hm (Z) —) Hm,i (3.49) 62 where —1—- (3.50) h N—l and the additional subscripted i on each modal coefficient corresponds to the index of the corresponding discretized node z,~ = (i- l)h (3.51) 4 Substituting into (3.40), multiplying through by (Zh/y) , and collecting nodal terms yields the finite difference equations for the central nodes, COWm,i—2 + Cle,i-l + CZWm,i + C3Wm,i+l + C4Wm,i+2 = Fm,i (3.52) with the resulting modal matrix coefficients Co = C4 5 B (3.53) CI = C3 5 _(2Ba,,,2 +1)h2 - 48 (3.54) _ 2 2 4 2 2 2 ~ 4 C2 = am (Bam + l)h + 2h (28am +1)+ 63+ 7 MAW-h (3.55) This discretization produces a symmetric penta-diagonal system for each modal momentum equation, which can be solved sequentially for each mode value. A basic Gaussian elimination method is employed, adapted from the code used by Roebers (1986), which takes advantage of the compactly banded system. The right-hand-side terms are calculated as ’ K 2 .. Fm, a -;h4 R2 (aszkJAkm)— MQW- (3.56) k=l 63 which are computed from the T“ values from the previous iteration. Now the discretizations of the linearized coefficient in the diagonal term, ZmJ- , and of the remainder on the right-hand side come directly from the discretization of the modal variables: K mZZdH H,k iAkmm(0) (357) (3.58) The first derivatives appearing for Hm and V“ are also calculated from the values of the previous iteration using second-order correct finite differences, with the central difference employed for the central locations, i=2 through i=(N-1); and three-point forward or backward differences applied at the bottom and top wall nodes. The near-boundary terms require special handling, since the first and last discretized momentum equations would refer to terms outside the domain. The impermeable wall boundary condition requires a value of zero on the boundary nodes, Wm,l = Wm,N = 0 (3.59) which also eliminates the need to solve the first and last node finite difference for each mode. For a second boundary condition on each wall, as required by the additional orders of differentiation, we look to the no-slip conditions for v 5'?"— = é"— : 0 (3.60) 6} z: 0 d) z: I combined with the continuity equation 0’12 5w 5 + E = 0 (3.61) which yields a zero-slope condition for w at each wall 0% fl: _é1: = 0 (3.62) 2:0 a" 2:1 The discretization of these conditions is applied as the ”reflected node” method for nodes adjacent to the boundaries. Wm,0 '9 Wm,2 (3 63) Wm,N+l 9 Wm,N—1 ' The first and last nodal equations, corresponding to i=2 and i=N—1, are then (C0 + C2)Wm,2 + C3Wm,3 + C4Wm,4 = Fm,2 (3.64) COWm,N-3 + Cle,N-2 +(C2+C4)Wm,N-1 = Fm,N-l (3-55) 65 l 3.3.2 Zeroth Mode Energy Equation The zeroth mode energy equation is similarly discretized with the following second-order correct finite difference forms: (Anderson, 1984) 70(2)» 70,, (3.66) 7 1 To (4')" 2(T0,1+1 - 70,1—1) (3-57) " 1 To (21)" 25(761” ‘ 2T0,i + 70.1-1) (358) which yields the tri-diagonal form: BIT0,i-l + 32T0,1+ B3TO,i+l = 4721303 (3-59) with the coefficients defined W K Ii B E -h ’ - ' E1 a, 7 (3.70) [odd 32 5 27 (3.71) K W _ l 1 B E +h ’ - 3 E a, K (3.72) [odd and the right-hand terms defined as ' 1 . K m K FOJ = Z ‘ TVkJTkJ + Z AkIWIJTIéJ (3-73) 1‘3] [:1 _ (Ir-l) odd The boundary node values are determined from the known temperature boundary conditions, To,1=1 ; To,N=0 (3.74) 66 so the first and last nodal equations to be solved correspond to those adjacent to the boundaries, expressed as 82 T02 + 3372),”) = -h2fi0,2 - B] (3.75) BIT0,N-2 + 82T0,N_] = _h2fi0,N-l (3.76) The resulting non-symmetric tri-diagonal system is then solved using a common tri-diagonal solver, in this case adapted from Roebers (1986). 3.3.3 Modal Energy Equations The energy equations in the higher Galerkin modes are discretized in the same manner as in the zeroth mode equation, with the following second- order correct finite difference forms: (Anderson, 1984) Tm(z) _) TmJ (3.77) 7 l Tm (2,)» firm,” — Tm,,-_ 1) (3.78) n 1 Tm (Zi)—) 23(Tm,i+l " 2Tm,i "l' Tm,i—l) (3.79) which yields the tri-diagonal form: Bl,iTm,i-l + 32,1 Tm,i + B3,iTm,i+l = ‘2h2Ym,i (3-80) 67 with the coefficients defined K 81,.- 412 Amzmlz’Wz. — 7 (3.81) 1:1 K 82..- 1.2,. 2amZ(V1,.-A m1... 0) +27 (382) [=1 K 83,.- +112 Am1m(2’m,.- - 7 (3.831 1:1 and the right-hand terms defined as Ym, _ -—2h2 K F...- + 16.2 (11mm) )) (3.841 [:1 and K K K =2 [" “k Tk,i[ Z V1.1Aklm(l)) + Th) 2 W1,1Ak1m(2))) (3.35) k=l 1:1 1:1 k¢m The discretized energy equations are more straightforward at the boundaries, since the first and last equations only refer to the specified boundary condition nodes. Since the applied temperature boundary conditions are satisfied in the zeroth mode, the remaining modes should vanish at the top and bottom walls, yielding the following equations corresponding to the nodes adjacent to each wall: 822 Tm,2 '1' 33,2 Tm,3 = -2h2Ym’2 (3.86) 2 31,N-1Tm.N-2 + BZ,N—1Tm,N—l = “2k Ym.N-1 (3-37) 68 3.3.4 Expansion of the Nonlinear Term As laid out in section 3.1.4, the nonlinear term (the velocity magnitude) in the Forchheimer term of the momentum equation is represented as a series of Fourier-Galerkin coefficients to facilitate the linearizations just employed. These coefficients are calculated in a manner similar to the discrete Fourier transform. First the field of velocity components and the corresponding magnitude values are computed at a number of sample points, NV, across the y-dimension. These are computed from the most recent guess or solution values in the Ansatz, (3.1), (3.2), and (3.20). To get the calculation for the expansion of 1), first form the inner product of (3.20) with the basis functions, 7 7 (77(y,z)sinamy dy = (H), (z) sinakysinamy dy 0 () (3.88) = £11142) Now solving for H, and discretizing Hk(z,-)—) H,” (3.69) yields the calculation 2 NY Hu- = N 1 Z 4(y1,z.-)sin(aky1) (3.90) Y ' 1:1 For additional computational efficiency, the sin(a,,y,) terms are computed once in the setup of the program, and recalled from storage for the field value recovery and nonlinear expansion calculations. 69 3.4 Calculation to Recover Horizontal Velocity After solving for the discretized vertical velocity modal components, Wk, the corresponding horizontal velocity can be calculated using the Galerkin form of the mass conservation equation (3.11). After applying the finite differences for the z-derivatives, the V“ components are computed as 2 K . Vm.1 = - — Z Wu Amk (3.90) mflk-l These components are calculated following the solution of the vertical velocity components, W1..- during each iteration, for use in the subsequent calculation of the right-hand inner products in the solution of the thermal energy equation. 3.5 Post-Processing After the iterations are converged, the post-processing requires calculation of several quantities, including 0 pressure coefficients - stream function - Nusselt number 0 wavenumber - entropy generation rate - generalized excess entropy production functional 70 3.5.1 Computing the Pressure The Galerkin coefficients for the pressure field are necessary to compute the global entropy generation rate, G, and the generalized excess entropy production functional, ‘1’. These pressure coefficients are determined by taking the y—momentum equation (2.36), substituting the Galerkin expansions (3.1)-(3.3), (3.5), and (3.20), and forming the inner product with the sin any basis functions, resulting in the equation . K K - 7 Z akPk sinaky+ BX (- akZVk + Vk )sinaky k=l k=l , j K K smamy dy : 0 (3.91) 0 — Z Hk sinaky) 2 V, sinaly] _ k=l 1:1 _ which reduces to mzr K K ( ) [Pm + B(- amZVm + V,;;)- VJ?) - M2 2 HleA 0 k1,, : 0 (3.92) k=11=1 This equation can be solved directly to compute the coefficients for pressure, Pm , in terms of V", and H,. K K Pm : Vm + B(am2Vm - V,;;)+(;12—-)MZ Z HkV,A(O)/dm (3.93) k=11=1 71 3.5.2 Stream Function Calculation The stream function can be calculated from the velocity solution, and is useful for visualizing the flow pattern. For the steady state flow, no fluid crosses the streamlines, which are the contours of the stream function. Most fluid mechanics and computational fluid dynamics texts describe the stream function formulation of the governing equations, but do not describe the method for computing the stream function directly from the known velocity. The procedure taken in this work follows. Beginning with the usual definition of the stream function, 1|1, such that its partial derivatives identically satisfy continuity. 11.1,.)(wsnggj Integrating the second condition (involving w) yields: y w = — )w dy+ f(z) (3.95) 0 where f is some arbitrary function of 2. Next differentiate with respect to z, with the result being equal to v, according to the definition of t|1. 01/1 0% df sz=-IEdJ)+—CIZ— (3'96) But invoking continuity again, 3’3- 9: 397 {-0} 1. 1 so v=+)flajz+g[z- é) (3.98) _ if: —.v+dz Therefore f must be a constant, which can be set arbitrarily for the convenient assignment of the stream function constants. Choosing the most convenient, f = O (3.99) leaves the value of ()1 as {/1: - wdy (3.100) Oh—-. To compute this value across the flow field, we substitute the Galerkin Ansatz and evaluate at each location fldeysinak (3.101) 1:1)(1- cosaky ’) ._. M; 73 3.5.3 Nusselt Number Calculation The average heat flux at the wall is represented by the Nusselt number, defined as at _ E- ‘7 L km Aft/Cm Nu (3.102) where the starred quantities again represent the dimensional quantities, and 1' represents the average heat flux at the wall. This can be calculated using the Fourier law applied at the top or bottom wall, and averaged across the width of the box, W. 4: q‘:fil/'J’:l ‘kmart 4: dy (3.103) wall Applying the scaling scheme from section 2.3 results in the expression _*- kmAT* 7g q ' W 0&- dy (3.104) wall Substituting the Galerkin temperature expansion yields . kmA T" W dy wall sinaky‘g] (3-105) 7 a K (0 22- To +E1Tk cosaky K 1 ' + Z —Tk wall k=lak Ql ll 7T0 wall 74 When this expression is substituted back into the Nusselt number definition and evaluated at the top and bottom walls, the elegant result obtains Nulbotmm = —T0'(0) (3.106) I Nulmp = - To (1) (3.107) 3.5.4 Entropy Generation Rate Hypothesized as a predictor of the steady-state situation in the infinite width porous layer, the global entropy generation rate, G, represents the integral of the rate of entropy generation, Sgen, over the entire domain {ye(0,y), zE(0,1)}. To calculate this value we begin by integrating equation (2.66): (Ti r)[li2|2 - B(ii-V2fi)]}dya’z \ (3.108) tdydz @2 &2 @2 &2 J 1 (.2..2-.(.(sa:v..e2_v)..(aw.6%)} Assuming the temperature fluctuations are small in the y-direction, the (T+t) coefficients are passed outside the y-integrals, with only the mode zero term, (To+1:). The integral is then broken down into six component terms, the 75 i Galerkin ansatz substituted and simplified using the inner product terms. The value of G can then be directly computed from the Galerkin coefficients of the velocity and temperature, with Simpson-rule integration across the z- dimension, as l G— __ g— (T0+BT)2G 1 l+R(T0+ T) (To + r)2 [02 - BG3]+ %[-—To— (To + 1) G6 dz where the component terms, Gl 0,, are defined and simplified as 61(2) 5 Gs (Z) (3)2131) 76 (3.109) (3.110) (3.111) (3.112) (3.113) and (3.114) K 12 12 = 2(Vk +Wk J k=l NIV Further details of the substitutions and simplifications can be found in Appendix C. 3.5.5 Generalized Excess Entropy Production Functional '- ' The generalized excess entropy production functional, ‘1’, is also calculated by integrating over the entire domain, as indicated in equation (2.90). As with the calculation of G, the appropriate Galerkin expansion is substituted for the velocity and temperature factors; the expression is simplified using the inner product coefficients; and the expression is integrated across the z-direction using a Simpson’s rule routine. Breaking the expression down into several component terms of the integrand, the ‘1’ expression can be computed as 1 )1! = “HBMSl + Sz)+ S3 - S4 + B(— 55 + S6)+ S7 - M58]dz (3.115) 0 with the component terms, SI 58 , defined as 7 am) am) S (2) a b7{v w—Jay 1 (J; a» + a} (3.116) K K K ' = Z Z 21- Aklm(l)VleTm 1' Aklm(2)Tk me k=ll=lm= 77 and _ ’ (3(a) 2 01m 2 W = ((7) l a ) dy K (3.117) = 12(ak2Tk2 + Tk ) 2k=l 7 53(2) 5 (((TO + r)[5v01;P) + 51:) 6p] - %w&}ay 0K (3.118) : - 2 {:(:0(2 + r)akaPk + Z Akak[( To + 7),,1W + Ton ]} k=l m: l )=:[(T0+ T)W (3.119) K K = (To + le Z Akaka k=1m=l 7 a? 5v a? 6v 2 5w 2 SsMEl(To+r)( 0;wa 0.121%) 1%)) )2 0 K (3.120) =-%(T0+T)kzzl[a (Vk +Wk )+Vk +Wk :l 7’ 3.1. ”(.2132 33—3)]... 0 K (3.121) = 1T Toz—lleVk + Wka I 7 S7(z)= = “To + r)(5v2 + (Swz)dy 0 (3.122) 2: and,fina"y 7 53(2) 5 ((TO + t)li2|"' ‘(M + 61v2)dy 0 (3.123) K K K : (To + I); 12: Z Aklm(0)Hk(Vle + l’VIWm) =1 =1 1 m: Further details of the derivations for each term are found in Appendix D. 3.6 Numerical Solution Algorithm The following pseudo-code lays out the structure of the calculations involved in the numerical solution routine for the velocity and temperature profiles in the porous enclosure. In this context, the superscript-starred notation on the discretized Galerkin coefficients indicates the solutions to the quasi-linearized finite difference forms of the Galerkin momentum and thermal energy conservation equations prior to updating the guess values. A Newton-iteration damping factor, (5w or (ST, is then applied, as indicated in sequence steps 3 and 4, prior to updating the new guess values. Also before updating the new guess values, convergence is evaluated by searching for the maximum absolute value difference between the solution values from the current iteration and the previous values. If any of the difference values exceeds the convergence tolerance parameter, another iteration is initiated. 79 The numerical algorithm sequence follows: 0 - Initialize Program.& Set Parameters Read values from input file calculate Galerkin basis functions calculate inner products A1m calculate triple products Ahmmh Ifihf“, Ahmu’ 1 - Make initial guesses for temperature and velocity fields Tk , Wk, or read in from file Calculate Vk from Wk [ DO j=0 (until convergence...) ] 2 - Compute zeroth.mode (mean field) of Temperature 3 - Compute velocity solution (check if Forchheimer extension is used) 0 compute Fourier sine expansion of the nonlinear term nhi(y) into Fourier sine coefficients Hki [ DO m=l to K ] ‘ 0 Compute the nonlinear inner product C2,“,i from Vk I wk I Hk,.il Aklmm)! Aklmm Compute thermal RHS terms from 1k, Assemble momentum matrix coefficients Solve for Mm; Update guess: W,“ = 6,, Wmf + (1- Saw“ check for pointwise convergence [ END DO (m) ] 0 Calculate Vhi from Whi (using Galerkin mass continuity) 4 - Compute temperature solution [ DO m=l to K ] 0 Compute energy inner product RHS Em, from Vkl Wk! Tkl Aklmm: Aklmm Assemble energy matrix coefficients. Solve for flmf Update guess: ng = éT‘flmf4 (1- 6T)flmi check for pointwise convergence [ END DO (m) ] [ Increment j and REPEAT until converged ] 80 5 Post-calculate stream function (display streamlines and temperatures) pressure field Nusselt number excess entropy production rate, PSI entropy generation rate, G 81 Chapter 4: Results The numerical simulation model solution method developed in the preceding sections was implemented and carried out on a desktop PC digital computer. The behavior of the solutions and results leading to the evaluation of the wavenumber prediction are presented first, followed by a validation study to evaluate satisfactory convergence with the numerical parameters. 4.1 Effect of Box Width In order to establish the feasibility of the code solutions, and to examine the behaviors as the box model is widened (as gamma is increased), the sumulations were first run with the Forchhiemer term coefficient, M, initially set to zero. This reduces to the linear case of Darcy flow, with the Brinkman extension. The numerical solutions for situations above the first critical Rayleigh number exhibit the anticipated cellular convection pattern, with the size of the cells increasing (corresponding to a decrease in the wavenumber) as the Rayleigh number is increased, which is consistent with previous experimental and theoretical results. The cellular patterns and the trends in cell size can be observed in the streamline plots of Figures 4.1 - 4.3, 4.7 - 4.9, and 4.13 - 4.16. The accompanying temperature fields are displaced from the uniform linear temperature profile in patters consistent with the flow field, as seen in the 82 isotherm plots in Figures 4.4 - 4.6, 4.10 - 4.12, and 4.17 - 4.20. The streamline plots were generated from contours of the stream function field, which is calculated from the Galerkin coefficients, Wk(z), according to equation (3.101). The isotherm plots are generated from the contours of the temperature field values, recovered using the Galerkin expansion (3.3). Note that these figures are not to scale, as the aspect ratio, y, varies for each set. For a given box size, one can observe the expansion of the cells as the Rayleigh number is increased. For the small box, y=3, the central roll becomes wider and more tilted with increasing Ram, pushing the small recirculating rolls into the corners. For the y=8 box, we see that one roll widens with Ram, as the adjacent roll narrows. The net effect, however, is for the cell, consisting of the roll pair, to increase in width (wavelength), with a corresponding decrease in the wavenumber. A sudden change in the roll shapes occurs around Ram of 90, and for Ram=100 we can see only 5 distinct rolls, instead of the 7 observed at Ram=50 and 75. For the widest box, with y=15, there is a similar pattern. At this aspect ratio the planform shift occurs around Ram=70, with 13 rolls observed here at Ram=50 and only 11 rolls at Ram=75, 100, and 120. The cells again show a widening with increasing Rayleigh number, as one roll widens and the adjacent rolls become narrower. The slight tilting of the rolls also becomes more pronounced as Rayleigh number increases, and the rolls against the left and right walls are squeezed into the corners. 83 41 33 25 Z index 1 101 ‘ 201 301 Q01 1 501 Y index Figure 4.1 Streamlines for Ram=50, y=3 41 33 25 Z index 201‘ l V 301; 3401 ‘ 501 Yindex Figure 4.2 Streamlines for Ram=75, y=3 101 41 33 25 Z index Y index Figure 4.3 Streamlines for Ram=100, y=3 41 1.000 33 0.875 0.750 g 25 ‘2 0.625 N 1., 0.500 0.375 9 0.250 0.125 1 0.000 Y index Figure 4.4 Isotherms for Ram=50, y=3 41 1.000 33 0.375 0.750 g 25 0.625 C “ 0.500 N 17 0.375 9 0.250 0.125 1 0.000 1 101 201 301 401 501 Y index Figure 4.5 Isotherms for Ram=75, y=3 41 1.000 33 0.875 0.750 X g 25 0.625 N 17 0.500 0.375 9 0.250 0.125 1 0.000 1 101 201 301 401 501 Y index Figure 4.6 Isotherms for Ram=100, y=3 85 33 fig 25 N 17 9 1 q .. , ‘ I ' ' 1 ' 1 1’ ' ' ‘ ‘ 1 ‘101 201 301 401 501 Yindex Figure 4.7 Streamlines for Ram=50, y=8 41 77777 33 2 * ... 8 25 1 .s , N 17 l. 9 ~ .1 fr'. ‘1 101 '201 301 ’501 Yindex Figure 4.8 Streamlines for Ram=75, y=8 41 33 :: _ “M“ g 25 * I ‘2‘ ii 17 9 1 Yindex Figure 4.9 Streamlines for Ram=100, y=8 41 1.000 0.875 0.750 0.625 0.500 0.375 0.250 0.125 0.000 33 25 Z index Y index Figure 4.10 Isotherms for Ram=50, y=8 41 33 25 Z index Y index Figure 4.11 Isotherms for Ram=75, y=8 41 1.000 0.875 0.750 0.625 0.500 0.375 0.250 0.125 0.000 33 25 Z index Y index Figure 4.12 Isotherms for Ram=100, y=8 87 Z index N VI Y index Figure 4.13 Streamlines for Ram=50, y=15 Z index N 0| Y index Figure 4.14 Streamlines for Ram=75, y=15 Z index N U! Y index Figure 4.15 Streamlines for Ram=100, y=15 Z index to UI -. . .__. . .. .-.. r ‘ ‘ " 1.. . “‘ . ‘ 1 101 201 301 401 501 Yindex Figure 4.16 Streamlines for Ram=120, y=15 88 41 33 25 Zindex 1 101 201 Y index Figure 4.17 Isotherms for Ram=50, y=15 41 33 Z index 1 1 01 201 30 1 401 Y index Figure 4.18 Isotherms for Ram=75, y=15 41 1.000 33 0.375 0.750 g 25 0.025 N 17 0.500 0.375 9 0.250 0.125 1 0.000 Y index Figure 4.19 Isotherms for Ram: 100, 11:15 41 Z index 1 101 201 301 401 Y index Figure 4.20 Isotherms for Ram=120, y=15 89 4.1.1 Onset of Convection One of the first effects of the side walls that is observed is that the side walls of the box produce stabilizing effect, especially for the narrower boxes. In the course of this study it was observed that the porous media Rayleigh number required to initiate convective motion rises above the accepted critical value for the infinite porous layer of 41tzz39.5. This trend is illustrated in Figure 4.21, which contains a plot of the Rayleigh number bounds between which convective motion began, over the range of aspect ratios simulated. The initiation of convective motion was indicated by the Nusselt number remaining at a constant 1.0 — indicating pure conduction - for Rayleigh numbers up to the lower bound, and a Nusselt number greater than 1.0 at the upper bound Rayleigh number. The values of the entropy generation rate, G, and the generalized excess entropy production functional, ‘1’, also increased by several orders of magnitude between these Rayleigh number values, indicating the sudden presence of convective motion. From these results it is observed that for simulations with an aspect ratio of 10 or wider, the critical porous media Rayleigh number remains close to the 39.5 value, and remains below 40.0 for aspect ratios wider than 12. 90 b b 50 -~ - p \ h D 48 q: 1'11. 1 o 3. \ 'ff‘ 1 I 1: . . \ 46 -- .1 1 i I porous media Rayleigh Number, Ran fi 42 '1: ’ \ . : ‘~2 40 d)- Lm-ascg‘n :3: . :01.‘ ................ _ i; *‘w..1;.m_‘ w ‘‘‘‘‘ ”Aw-"w” 38 A I 5 I T L I l L I I A T L I 0 2 4 6 8 10 1 2 14 1 6 1 8 20 aspect ratio, 37 Figure 4.21 Bounds of Critical Rayleigh Number for Convective Motion 4.1.2 Nusselt Number The effect of the side walls is also seen in the resulting flow and heat transfer behavior at Rayleigh numbers above critical. It is expected that the Nusselt number maintains a constant value of 1.0 for Rayleigh numbers below the first critical value, and that it will increase monotonically with further increases in Rayleigh number. For the purposes of this study, it is also necessary to determine the minimum aspect ratio above which the effect of the side walls exhibit negligible influence on the Nusselt number. The effect of the side walls on the heat transfer rate is illustrated in Figure 4.22, which shows the Nusselt number - Rayleigh number correlation for simulations of successively wider 91 boxes. In order to validate the results, published correlations are also included from the correlation of Somerton (1983) to the data of Combarnous, Ram(0.81) N“: 14.12 for Prm > 11.8 (4-1) and the simple result obtained by Bejan (1984) from a scaling analysis: 1 if Ram<40 Nu= Ram . (4-2) LIIIO .xblOlooe. Nusselt Number. Nu N 0| -0 20 4O 60 80 100 120 140 160 180 200 220 240 Porous Media Rayleigh Number, Ram Figure 4.22 Heat Transfer For Linear Simulations 92 The increase in Nu as Ram increases is immediately apparent. Further inspection reveals that the initial Nu vs. Ra"n slope increases with increasing aspect ratios. For a given Ram value, the Nusselt number increases significantly with an increase in the box width, and appears to approach a limiting value of an ”effectively infinite” width value in the range of y from 15 to 20. 4.1.3 Wavenumber Measurement Method The question of how to evaluate the convection cell size in order to compute the wavelength of convection required evaluating several possible methods. While the Galerkin expansion coefficients are themselves a spectral result, they occur only in integer increments, and characterize the flow in the entire width of the box. On the other hand, the wavenumber varies to a fractional degree, and characterizes the periodic cells. Several approaches were considered. 1. Searching across the box for zeroes of the horizontal velocity field component, v, which could identify the edges of a roll. 2. Searching across the box for zeroes of the vertical velocity field component, w, which could identify the centers of a roll. 3. Searching across the box for local maxima and minima of the vertical velocity field component, w, which could identify edges of a roll. 93 The first two methods were found to yield inconsistent results. The first method failed as the rolls were found to tilt to the left or right, so the anticipated verticality of the flow did not actually occur at the boundaries. In the second method, flow zeroes were often not well defined by the numerical results in the slow-moving eye of the cell. The third method was chosen as most logical and had the most consistent success in identifying roll boundaries. However, it was discovered that small oscillations near the side walls and in the centers of the rolls were also identified as local extrema, and incorrectly counted as roll boundaries. These points were filtered out by requiring the flow velocity component to be greater than 10% of the maximum value of w in the box in order to be considered as a roll boundary. This eliminated the slight perturbations in the flow found in the eyes of the cells, as well as small the recirculating regions along the sides walls and corners of the box. In effect, the 10% level established a (somewhat ad hoc) criterion for these smaller recirculations to qualify as bona fide rolls. In implementing the w-extremum search, the first differences of w at each node along the horizontal mid-plane of the box were examined from left to right. If the difference values were found to change sign, and the velocity component exceeds the 10% threshold, the location of the minimum or maximum was estimated by parabolic interpolation (Press, 1992). 94 Once the roll boundaries are identified by the search, it remains to evaluate the average wavenumber. This is accomplished by two measures, an overall average wavelength, 10A , and a central average wavelength, Ac". If only one roll was identified, the overall wavelength was assumed to be the entire width of the box. If multiple rolls were counted, the overall average wavelength is found by taking the distance between the left-most roll boundary located and the right-most boundary that completed a number of complete roll pairs, and dividing by the number of complete cells (two-roll pairs) counted. Since the rolls typically occurred in odd numbers, centered in the box, this measure was repeated between the right-most roll boundary location and the second boundary location from the left. These two values are then averaged together. This has the effect of un-weighting the measure of the left-most and right-most rolls by half. The overall wavenumber, a0A is then calculated as: 277 a 0 A = TOA (4.3) In order to minimize the retarding effect of the side walls, and to hopefully better capture the behavior that would occur in the infinite width layer, the central average wavelength, Ac" , was computed by taking the distance between the second roll boundary located and the location that completed a number of complete roll-pair cells, but remaining at least one (typically two) rolls away from the right side wall, and dividing by the number of roll-pairs 95 fi- counted over this distance. The central average wavenumber is then computed 27: actr : A (4.4) ctr 4.1.4 Wavenumber Results To evaluate differences between the two wavenumber measurement methods, and to look for convergence in the wavenumber with increasing aspect ratios, both methods of wavenumber measurement were applied to each simulation run. Values of the central wavenumber and overall wavenumber are plotted in figures 4.23 and 4.24 with increasing Rayleigh numbers for the range of aspect ratios studied from 2 to 25. These plots clearly demonstrate the decrease of the wavenumber with increasing Rayleigh number, which is consistent with the increasing cell widths (wavelengths) seen in the streamline plots. We can see that the a vs. Ra"n characteristics tend, in general, to converge as gamma increases. 96 . gamma=5 1 -~- gamma=6 + gamma=7 ...._ gamma=8 --— gamma=9 —- gamma=10 -o— gamma=12 : . gamma=15 l + gamma=18 -— gamma=20 1 + gamma=25 Central Waveumber.a ctr 30 40 50 60 7O 80 90 100 110 120 130 140 150 160 170 180 Porous Media Rayleig'l Nunber. Ram Figure 4.23 Central Average Wavenumber vs. Rayleigh Number 5 _ ,, 7 . ., a... . . , 1 W +gamma=2 —e— gamma=3 4.5 —.— gamma=4 3 gamma=5 4 -—- gamma=6 . + gamma=7 _.._ gamma=8 5 3.5 ——gamma=9 C. \ --—- gamma=10 5 “4. ‘ c s Ar 5 A e e c 4v c e c s c e gamma=12 a 3 ‘ _\I I I I - gamma=15 g &_h- q + 90mm“; 2'5 \e . H.— X gamma= > ° 8 ‘ ‘\ " 2%.. x \ I gamma=25 g \ I U ‘2 O I g I \g I X 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 Porous Media Rayleigh Number. Ram Figure 4.24 Overall Average Wavenumber vs. Rayleigh Number 97 1m? -- This trend is better observed in Figure 4.25, where the wavenumber vs. gamma is plotted for several values of Ram. Here we can see that the overall average wavenumber measure is consistently slightly higher than the central measure. This is presumably caused by the smaller rolls near the walls which are included in the overall average. Both measures are seen to be reasonably close, within 10% at most, and averaging 2.6% difference. The Ram=50 values show a clear trend toward convergence to a ”wide box" value of 3.0. 4., 4 - 2 , , 3.5 3 3 c 05 B 2.5 c If; 2 "‘A‘" Ra=50. central .2 +Ra=50.0A 3 1.5 00" Re=75. central 3 § +Re=75.0A '3' 1 --e~Re-100.centrel +Ra-100,0A 0.5 «3.- Ra=120. central —e—Re=120.0A 0. a 35554511+1111111111 4. 0 5 10 15 20 25 aspect ratio. 3! Figure 4.25 Comparison of Overall and Central Average Wavenumbers For the Ram=75, 100, and 120 values, it is interesting to note that they all converge toward a common value at aspect ratios of 15 and wider, though all continue to vary together in the same way with increasing aspect ratio. It is 98 believed these jumps are caused by the sudden changes in the roll planform, as the number or shape of the rolls in the box takes sudden changes. 4-2 Non-Linear Results Next the non-linear parameters were varied to determine their effect on the wavenumber and Nusselt number of the system. The non-linear cases were run with a fixed gamma=15. Choices of the parameters for the Brinkman and Forchheimer momentum terms, 8 and M, were based on the ranges identified by Garcia (1991) for a wide range of physical cases of fluids, layer widths, porous media solid materials, and particle sizes Those ranges, and the Corresponding possible variations in M and B, are summarized in Table 4.1. The range of B is assumed to be similar to the range of the Darcy number, assuming the porous media viscosity ratio, cm, is of order one. Table 4.1 Ranges for Non-Linear Parameters Jiarameter Minimum Maximum Efrcy number, Da 3:110'12 3X1045 Forchheimer Number, Fom 0.09 0.3 \ Porous Media Prandtl Number, Prm 10'3 10‘ \ Brinkman Coefficient, s = o,,, Da 3x1012 3x106 \ F JD Forchheimer Coefficient, M = 3%7—0 2.7><10-17 0.52 m 99 4.2.1 Forchheimer effect Solutions were obtained for a range of Ram from onset of convection toward 180, as, the Forchheimer parameter, M, was successively increased, to determine the effect, if any, of the Forchheimer term on the heat transfer rate and wavenumber of convection. The resulting central average wavenumber and Nusselt number variations with Ram are plotted in Figure 4.26 for several values of M. 4 -_ ’1 r 9. i 3.5 ~_ 5' l 3‘ l . l " \‘. . ! 6 2-5 ‘: as l c _ v.‘:“~ 1 'o i ' 1 g 2 ‘: —.— Nu @ M=0 ; 3 t —0— Nu @ M=2E-4 Z 1 5 -1 +Nu @M=2E-3 ; ° C ..._ Nu @ 14:55-3 '1 C —I-— NU @ M=2E-2 1 1 4: mm a_ctr@M=0 I : m... a_ctr @ M=2E~4 0 5 .4: "'5'" 8_CU@M=2E-3 ; ' ~ ----- a_ctr @ M=5E-3 1 1 ...... a_ctr@M=2E-2 O L L L a TA L L 1 1 : 1 L L l + 1 1 1 A ? i 0 50 100 150 200 Ram Figure 4.26 Effect of the Forchheimer Parameter, M It can be seen in the plot that departure from the linear (M=0) solution is not observed for most of the range of Ra examined, until M approaches its extreme upper limit. The Nusselt number appears relatively unaffected for most of the range of parameters, exhibiting little change until the extreme 100 case of M=0.02, where it there is at most a 9% reduction from the linear case at Ram=80. There is some variation observed in the roll-shift wavenumber transitions seen around Rayleigh numbers of 70, 150, and 190. 4.2.2 Brinkman Effect Solutions were obtained for a range of Ram from onset of convection toward a value of 180, as the Brinkman term parameter, 8, was successively increased, to determine the effect, if any, of the Brinkman term on the heat transfer rate and wavenumber of convection. The resulting central average wavenumber and Nusselt number variations with Ram are plotted in Figure 4.27 for several values of B. 4 "z . 1 3.5 4E 1 I l 3 «3 “K. : \ 3 2.5 «E 1 . l 2 = a 2 21 . g t 1.5 «g .4— Nu@8=1E-6 1 ; --o»-NU@B=1E-S 1 1 .: --t-NU@B=1E-4 1 : ~~0---a_ctr@B=1E-6 : .....-a_ctr@s=1E-s 1‘ 0.5 1 ---4—--a_ctr@8=1E-4 1 l- l 0 L 1 A L # A 1 n n 11L 1 1 1 1 Ar 1 Ar - 0 50 100 150 200 Ram Figure 4.27 Effect of the Brinkman Parameter, B 101 l- Again, we see little influence of the Brinkman parameter on the wavenumber or the Nusselt number for most of the range of B values, until the extreme case of 8:104, which actually exceeds the predicted range of physically realizable values. We can also notice some influence on the roll shift changes in wavenumber around Rayleigh numbers of 70 and 190. 4.3 Comparisons with Wavenumber Prediction Theories This section will compare results of the porous box simulation of this work with data from simulation of the infinite porous layer. The porous layer simulation is based on the same mixed Galerkin and finite difference method that is employed in the porous box simulations. It is based on the dissertation and subsequent work of Somerton (1982) and Somerton and Jimenez (1999). In that code, four different methods are used to attempt to predict the wavenumber needed to provide closure to the problem. 1. Constant wavenumber, 1:, from onset of convection 2. Maximization of the Nusselt number 3. Minimization of the generalized excess entropy production stability functional, ‘1’ 4. Maximization of the global rate of entropy generation, G The porous layer simulation code was updated and revised to compute the entropy generation rate and generalized excess entropy production stability 102 functional in a manner consistent with the work presented in this dissertation. Work from the porous layer simulation will be published in a separate paper. The results, together with those from the porous box simulations, are summarized in Table 4.2. Table 4.2 Results from Porous Layer and Porous Box Simulations Fixed Maximizing Wavenumber Maxlmizing Nu Minimizing ‘I’ Entropy Pbox (M=0.Y=15) Mia—“0" Ram Nu a Nu a Nu a Nu a Nu actr am 48 1.369 3.142 1.370 3.154 1.295 2.585 1.354 2.882 1.311 2.859 2.925 60 1.780 3.142 1.783 3.264 1.516 2.163 1.709 2.638 1.631 2.667 2.746 75 2.193 3.142 2.201 3.269 1.701 1.888 2.010 2.388 1.859 2.206 2.283 90 2.547 3.142 2.577 3.540 1.966 1.847 2.268 2.275 2.145 2.207 2.263 110 2.968 3.142 3.034 3.845 2.269 1.816 2.505 2.092 2.473 2.225 2.254 130 3.354 3.142 3.464 4.025 2.495 1.750 2.740 1.996 2.720 2.223 2.240 150 3.718 3.142 3.872 4.227 2.709 1.697 2.968 1.924 2.945 2.226 2.215 170 4.075 3.142 4.265 4.417 2.994 1.708 3.213 1.881 3.144 2.218 2.209 190 4.435 3.142 4.645 4.594 3.262 1.708 3.434 1.827 3.322 2.949 2.757 2E 4.799 3.142 5.017 4.598 3.646 1 .782 3.688 1.809 The Nusselt number-Rayleigh number relationship from each method is plotted in Figure 4.28. All of the methods employed correctly predict a trend of increasing Nusselt number with increasing Rayleigh number. The results follow generally over the prediction of Bejan (1984). All of the simulation results lie below the correlation given by Somerton (1983) for the high porous media Prandtl number case corresponding to M=O. This may be attributable to the two-dimensional roll restriction, which corresponds to a box with depth 103 less than its height, while the enclosure of the Combarnous experiments correlated by Somerton was deeper, resulting in three-dimensional rolls. 5.0 -- ’ _ + Fixed Wavenumber ——o—MaximizingNu —-»-— Minimizing PSI 4.0 ‘1 —+— Maximizing G : —-— PBox (gamma=15) A E. C - - Bejan 0 . Somerton g 3.0 «- 3 Z T: in g 2.0 T Z I . f 1")...J.[....;.....r....i 0 50 100 150 200 Rayleigh Number, Ra m Figure 4.28 Nusselt Number Comparison for Wavenumber Prediction Schemes Let us now restrict discussion to comparisons between the porous layer and porous box simulations. Clearly, the constant wavenumber and Nusselt number maximization schemes over-predict the slope of the relationship, over- predicting the heat transfer rate by as much as 34% and 40%, respectively, at Ram: 190. Certainly a more accurate method is wanted for engineering purposes. Both entropy-based measures yield much closer predictions of the Nusselt number. Compared with the result from the porous box, the entropy 104 generation rate method consistently over-predicts the Nusselt number by as much as 8% at Ram=75 and by only 3% at Ram: 190. The stability functional method consistently under-predicts the Nusselt number by as much as 9% at Ram=75 and by only 2% at Ram=190. The actual wavenumber predictions from each method are plotted against the Rayleigh number in Figure 4.29. While the few published experimental results and the porous box simulations tend to predict a decrease in the wavenumber with increasing Rayleigh number, the constant wavenumber method, of course, ignores this trend; and the Nusselt number maximization scheme actually predicts an increase in the wavenumber. Again, both entropy-based methods appear to more accurately predict the convective wavenumber. The entropy generation rate method and porous box predictions are within 10% up to Ram of 130, where the porous box simulations begin trending in the opposite direction. The stability functional method demonstrates a trend very similar to the entropy generation rate prediction, but predicts wavenumbers consistently lower than any other method, as much as 21% lower than the porous box simulation at Ram of 130, and as much as 21% lower than the entropy generation prediction over the range of Rayleigh numbers examined. Though the agreement in values produced by these methods do not exhibit conclusive evidence, the agreement in the trends does appear to suggest the 105 extremum of the entropy generation rate and of the generalized excess entropy production functional as valid predictors, and the accuracy of their resulting heat transfer predictions is very encouraging. 5.0 -;—~— ”W- +~ w m w mu w 4.5 4.0 -g <2 3.5 '- E g 3.0 «5 1 § 2.5 -E 5 : > 2.0 -E w ' - "Ree— .. u Mi ; 15 J +leed Wavenumber .3 a n - E -—o—-MaximizingNu 1 O ‘E --—Minimizing PSI ' E +MaximizingG Q5 -5 +PboxCentral Avg. : +Pbox Overall Avg. 0.0 ‘ ‘ J ‘ i ‘ L ‘ i 4 v ‘ i O 50 100 150 200 Porous Media Rayleigh Number, Ram Figure 4.29 Wavenumber vs. Rayleigh Number Under Various Prediction Schemes 106 4.4 Numerical Validation The code for the simulation, PBox, was written in Fortran 95, using the Compaq Visual Fortran 6.6c compiler and development environment, operating under Microsoft Windows XP-Professional Edition. Modern Fortran 90/95 array operations were employed wherever feasible. All real variables were stored and computed in double-precision (8 bytes). The Compaq Array Visualizer v1.6 was used to produce the streamline and isotherm plots. The code was executed on 3-GHz Intel Pentium-4 based Dell and MPC desktop computers with 1.0 GB of RAM. Solution times were reasonably rapid. Nonlinear solutions were achieved with 3000 iterations taking only a few minutes or less, when using 40 Galerkin modes, 41 finite difference nodes for each mode, and 501 horizontal samples for field calculations and the non- linear term expansions. Excessive virtual memory swapping to disk was not observed until the number of modes approached 60 during validation testing. The criteria for convergence was a maximum change of 1.0E-5 in any of the Wk, and Th, Galerkin mode coefficients after each iteration. Anticipating that, at least at the onset of convection, the rolls would be somewhat ”square" (Beck, 1972), it was desired that the maximum number of Galerkin modes, K, be chosen at least twice the aspect ratio of the box, to exceed twice the expected number of rolls. Thus, K =40 was used for the simulations reported in the results. 107 The number of finite difference nodes used in the solution, N=41, proved to be workable within the computational memory requirements, allowed for reasonably rapid solution times, and provided a reasonable resolution for post- processing and human interpretation of the output. It is subsequently shown here to be provide a satisfactorily accurate numerical solution, as well. The number of sample points in the y-direction, NV, was set at 501. This parameter controls the number of points at which the flow field is computed from the Galerkin solution, and was chosen to achieve a fine resolution in the wavenumber search of increments of 0.03 for a box with aspect ratio of 15. Parabolic interpolation used to locate the roll boundaries between the computed location further improved the precision of the wavenumber computation. This parameter, NV, also affected the sampling input to the discrete Fourier expansion used for the Forchheimer term, and well exceeded the necessary sample rate of the K =40 modes. To conclude that the numerical solution parameters used in the previous results are sufficient, a validation study of the PBox code results was performed to confirm convergence in the Nusselt number and wavenumber results as the number of Galerkin modes, K, and the number of vertical finite difference nodes, N, are increased. The maximum value of the vertical velocity field component was also examined in order to evaluate the convergence of the flow field solution. The runs were performed for the linear model with box 108 of aspect ratio, y, of 15. The results are shown in Figures 4.30 through 4.35. For the number of finite difference nodes, the results are plotted against the finite difference step size parameter, h. 1 =—— 4.6 h N_1 < ) Note that the N of 41 nodes used in the simulations that were reported in the results corresponds to h=0.025. From the plots we can see that, for Ram=50 and Ram=80, the values converge very well with the chosen numerical parameters. The Nusselt number, wavenumber, and maximum velocity curves are all relatively flat for K greater than 20, and for node spacings smaller than 0.05 (N=21). Note that at least 12 modes were required to allow the solution to model any convective motion, and at least 16 modes were needed to properly capture the motion of all of the rolls at this aspect ratio. For Ram=120, it is seen that the values are convergent with the number of Galerkin modes up to about 15, then become somewhat skittish. This suggests that some caution should be used in the results at higher Rayleigh numbers. A ”change-over" in the number of and configuration of the rolls found in the simulated flow field was observed in the simulation results in this neighborhood. In th the validation study, the number of rolls (those exceeding the 10% threshold) that were recorded for runs exceeding 15 modes ranged 109 from 7 to 13. It is suspected that the flattening out of the wavenumber results with increasing number is affected by this behavior, with the number of rolls appearing in the numerically stable simulated flow field is affected by the number of Galerkin-Fourier modes available in the expansion. The Ram=120 solutions, like the Ram=50 and Ram=80 cases, were well converged for all of the finite difference mesh sizes examined. 110 1 -+- Ra = 50 1 -e-- Ra=80 1 + 128le0 2.50 "r 1 2.00 "r 1 1 .50 ‘- 1 .00 ‘ ’- 0 3.00 ., ._... . o number of Galerkin modes. K Figure 4.30 Nusselt Number Convergence with Number of Galerkin Modes :3 Z 3.00 mm b -0— Re = 50 1' W +R8=120 2.50 1» 2.00 4 1.50 4 fl} 1_go W? ............................ ; ,,,,,,,,, 0'00 0-02 0.03 0.10 Figure 4.31 Nusselt Number Convergence with Finite Difference Mesh Size g 3.00 ‘, -—o—Ra=50 } E _ ,- 0 2.50 +Ra-80 ‘ -° +Ra=120 , E j! E 2.00 -- 9 / / g 1.50 d. 3 1.00 E 73' 0.50 O 0.001““i““+““i““i““ 0 10 20 30 40 50 number of Galerkin modes, K Figure 4.32 Wavenumber Convergence with Number of Galerkin Modes a 3.00 -; +Ra=50 g E W +Ra=120 a 2.50 ‘: E I g 2.00—E > t $150~§ - 5 ¢ 3,, 61 t ' 5 1.004 E *5 0.501 8 . 000 ’ ......... % ......... § ......... ; ........ % ........ 0.00 0.02 0.04 0.06 0.08 0.1 0 h Figure 4.33 Wavenumber Convergence with Finite Difference Mesh Size 112 25.00 -_-~—-——-——~~ -- 2, We WW 2 - _ 1 : + Ra = 50 t R8280 20.00 -~ + : + Raz120 \N‘ 15.00 ~ TTj’ I Wmax 10.00 -E 5.00 f A ' 0.00 " 0 10 20 30 40 50 number of Galerkin modes, K l J Y I l I Figure 4.34 Maximum Vertical Velocity Component Convergence with Number of Galerkin Modes 25.00 . -_._, “Ww-m_--2__~__ -_._ —o— Ra = 50 20.00 3: —-—-—Ra=120 i 6 15.00 ‘1 E I 3 1‘ 10.00 -* 1 500 -E c: e c AW v ¢ v 1- I 0.00 ........I ......... ;..........r.........? ........ 21 0.00 0.02 0.04 0.06 0.08 0.10 Figure 4.35 Maximum Vertical Velocity Component Convergence with Finite Difference Mesh Size 113 Chapter 5: Conclusions and Recommendations 5.1 Summary of Contributions In the work of this dissertation, a mixed Galerkin and finite difference numerical solution was produced to provide numerical-experiment data for the flow and heat transfer occurring in a wide porous box filled with fluid and heated from below, for the purpose of measuring the wavenumber of convection. The implementation in the program code, PBox, accurately computes the flow and temperature fields, and the resulting Nusselt number and convective wavenumber, and provides colored-graphical output to visualize the streamlines and isotherms. Two entropy-based measures, the global entropy generation rate, G, and the generalized excess entropy production stability functional, ‘1’, were derived for the case of the porous box, and formulated to be computed using the mixed Galerkin and finite difference numerical solutions. Extrema of these measures were used to determine the wavenumber of convection for the solution of the infinite-width porous layer. 114 5.2 Conclusions Several of the conclusions that were resulting from this research are summarized here. 1. The mixed Galerkin and finite difference numerical solution is an effective method for computing the flow and temperature fields in a system such as a porous box heated from below. It is reasonable efficient in computation time and memory consumption, and can be effectively implemented on an ordinary desktop computer with accurate results. For the numerical solution of the porous box, an aspect ratio of 15 is wide enough for its central region to model the behavior of an infinite porous layer for porous media Rayleigh numbers less than about 100. For higher Rayleigh numbers, a wider box model may be necessary to minimize the constraints of the side walls. The Forchheimer and Brinkman terms of the momentum equation have little effect on the Nusselt number and convective wavenumber of the porous box solution at porous media Rayleigh numbers up to 200, for all but the most extreme parameters describing physical porous flow systems. The extremum of global entropy generation rate, G, appears to be the best available measure to predict the wavenumber of convection for an infinite porous layer. It is computationally straight-forward, reasonably 115 intuitive from a thermodynamic perspective, and its application enables very accurate predictions of the heat transfer rate across the layer. The generalized entropy production rate stability functional provides a reasonable measure to predict the wavenumber of convection for an infinite porous layer, and results based on the predicted wavenumbers tend to under-predict the rate of heat transfer. However, the arguments for computing the fluctuation components of the flow, pressure, and temperature variables based on the higher-order terms of the Fourier- Galerkin expansion are somewhat empirical. 5.3 Recommendations for Future Study Although the results of this research provide significant contributions to the study of natural convection in porous media and in enclosures, there exist several areas for continued study, and for new directions based on these conclusions. 1. The PBox routine should be used to study the behavior of the box with wider aspect ratios to reduce the effects of the walls at Rayleigh numbers between 100 and 200. This would likely require more Galerkin modes to accurately model the additional rolls in the wider box. Because of memory limitations in the PC, implementing this may require employing one or more of the following options: a. Porting the code to a higher-powered computing platform with more system RAM, and possibly parallel-processors 116 b. Adding additional physical RAM to the PC c. Making more efficient use of memory in the PC, such as sparse system solvers and allocatable arrays. Refine the analysis of the onset of convection to more precise values for the critical Rayleigh number, and the corresponding critical wavenumber. Extend the study to the dependence not only of the side walls, but of the Brinkman and Forchheimer momentum term parameters, as well. Visualization of the code results could be further enhanced by employing quiver plots of the flow, and a ”heat-lines” plot of the heat flux throughout the box (Bejan, 1984). A study of the roll planform shifts that occur in the box across certain Rayleigh numbers could be conducted, to see how the addition or removal of the additional rolls correlates with the extrema in the entropy generation rate and generalized excess entropy production functional. Additional physical situations can be identified that would require additional information for closure in the model. Maximizing the rate of entropy generation should prove to be a useful predictor of many other natural phenomena. 117 Appendix A: Galerkin Inner Products of Sine and Cosine Bases The inner products of the half-wave sine and cosine basis functions used in the Galerkin expansions factor into many of the calculations and derivations for the porous box solution. The inner product for each combination pair is derived in this section. We begin with the inner product of sine 8 sine, which is orthogonal, as expected 7 7n Isinaky sin amydy= —Isin kflsmmél d6 0 7’ 0 Z—j(—) if k=mandk¢ 0 (A1) = 1 k 0 if k¢mork=Oorm=O by recalling the definition 7dr ak s — ; k = O,1,2,3... (A2) 7 and defining the substitution = a_kX 1. 19’. (9.. k 7 (A3) to achieve the typical form on the right-hand side of equation (A1). 118 The same procedure for the product of the cosines yields a similar result 7r Icosakycosamyd = gjcoskficosméldfl O O l (1M5) if k=mandk¢0 7r 2 (A4) 0 if k¢mork=00rm=0 Now, for the product of sin and cosine, 7 7r Icosakysinamy dy = Z- Icoskfisinmfl d6 (A5) 0 O N consider special cases individually. If m=0, then the product and integral are obviously zero. If k=0 and m¢0, then the integral of the sine term is cos(m7t) _ 1: [1— (— 1w] (A.6) m m 71' IsiandB = O For the non-zero coefficient terms, the result depends on whether there is a fill-wave of half-wave difference. Again referring to the integral tables yields [(ZX—fl—J if (k—m) isodd 7’ 7, m2 _ k2 [coskesinmeda = < (A7) 0 0 if k=mor(k-m) iseven 119 Expanding the odd result in partial fractions, and combining the results: 7 Isinakysinamydy = < 0 :(kl-m k- 1m) --( I)“ 1 72' l m 0 l if m=0 ork=m or (k—m) iseven if (k-m) isoddandm¢ 0 if k=Oandm¢0 120 (AB) Appendix B: Triple Inner Products The Galerkin formulations for the thermal energy equation and nonlinear momentum equation produce triple inner products of sine and cosine. Their values are derived in this section. Beginning with the triple-sine product: Aklmm) a Isinakysinaly sinamy dy O (8.1) 71' l sin kasinlesinmada 0 1. § where the definitions of a, and the substitution for 6 are again as used in Appendix A, eqns (AZ) and (A3). Now the sine product relation is applied with the k and [terms sin asinfl= %[cos(a- ,6) — cos(a+ ,8” (8.2) to get 71' AW“) = i [[cos(k - [)08inm6— cos(k + 1)9sinm9]d9 (3.3) 0 Next the cosine-sine product relation is employed cosasinfl= %[sin(a+ ,6) - sin(a- ,6” (3-4) 121 to reduce to simple integrable sine terms (0) = 1, fl 8in(k - 1+ m)6- sin(k -1— m)6 A d6 “"1 47! O - sin(k + 1+ m)9+ sin(k +1- mm] (3.5) Recalling the elementary definite integral 7. 1 % if k odd (sin/cc dx = - Zooskxlg = (8.6) 0 0 if k even and noting that the factors (k-l+m), (k-l-m), (k+l+m), and (k+I-m) will all be either odd or even at the same time, evaluating the integral for each term yields the result — ‘ - fk+l+ dd (0) )2” k-l+m k-l—m k+1+m+k+l—m 1( m)o AM" = (3.7) 0 if(k+l+ m) eveng For the product of two sines and one cosine, we again begin by substituting 8 from equation (A3) and applying the sine product relation (8.2), 7 Aklmm s Isinakysinalycosamyajz 0 (8.8) 71' = —”— ][cos(k - l)6cosm6- cos(k + l)6003m6]d6 27: 0 Now apply the definite integral 7’ 7r - _ (cosaxcosbx dx = /2 .lf a — b (8.9) 0 0 If a at b 122 to yield the result 7/4 ifm=i(k-l) -%' ifm=(k+l) (1) =( 7 AW A ifm=0andk¢0and1¢0 (310’ o flk=0ml=0 0 otherwise Also required the product of two cosines and a sine. Applying the substitution of 8 from equation (A3) and the sine-cosine product relation, (8.4), yields 7 Aklm(2) E (cosakysinalycosamy dy 0 (3.11) 7! = 3’— )[sinuc + [)6cosm6- sin(k - l)6cosm6]d6 272' 0 Again applying the sine-cosine product relation, (8.4), simplifies to A (2) y T sin(k+l+m)6+ sin(k+l—m)6 “"1 ‘ 47:0 — sin(k— 1+ m)9+ sin(k— 1- m)6 (8.12) And the definite sine integral, (8.6), again can be applied to evaluate the expression as _ - - 'fkl dd 2n k+l+m+k+l-m k—[+m k-l-m 1(+ +m)o Ame) = l (3.13) ) 0 if(k+l+m)even 123 Appendix C: Galerkin Integrals for Entropy Generation Rate This section provides the details for the calculation of the total global entropy generation rate, G, which represents the integral of the rate of entropy generation, Sm, over the entire domain {ye(0,y), zE(0,1)}. Integrating SNa , given by equation (2.66) 17 Ga Hsgendydz 00 l lllHBT1.)2wn2+.‘,,,:(. luv-Bl mlldydz 00 , H3 072 072 . (0.1) l)’ (T+1)2 [ —) (I???) —OO +_1_ 1 2+ 2—B [221+fl]+ w(§2w+§2w] ( R(T+T)v w véfl 01.2 @2 01,2 8y assuming that the temperature fluctuations are small in the y-direction, the (T+t) coefficients are passed outside the y-integrals, with only the mode zero term, (Tow) - which depends only on 2, remaining. The integral is then broken down into six component terms, Gl G6 ; the Galerkin ansatz are substituted; and the component integrals are simplified using the inner product terms. 124 The first component term, GI , contains the conduction heat flux terms. (($213)) K _;_2 (aszkz + T152) ((3.2) k I l 0 61(2) K 2 K 2 [Z-aka sinaky] +[Z Tk' cosaky] dy k=l k=l where the final step utilizes the orthogonality of the sine and cosine in the squared terms, as given by equations (3.9) and (3.14). The G2 term contains the velocity magnitude. 7 7 K 7- K 2 62(2) 5 [(122 + W2)dy= I [2 VIC sinaky] + [ 2 Wk sinaky] dy o 0 k2] k=l (C.3) K =1 The Brinkman terms with second y-derivatives are collected in G3. 7 v .32.. 03(2); (J; vgi-Z-+w-é)—2— dy 7 K K K K = ([ Z Vk Shawl: - 012V13inaly] + [ 2 Wk Sinakyflz ' 012W! Si"(WNW 0 k=l [=1 [(7-1 [=1 K .—. -22: Z ak2(Vk2 + M?) (CO4) 125 Now, for the Brinkman terms with second z-derivatives, we shall first define the integral, G, (which does not appear in the final computation), as the double-integral form 17 05”“) 1+T)[v£21+wi—:fijdydz ((3.5) 00( Next we reverse the order of integration and integrate by parts, to reduce the order of z-differentiation. r ' -To (vfl+wflj \ 7 1 [a my 1(7’0”) d7 0} G E _ _ _ )dz d 4 ({(To+r) v&+w& 2:0 (J): 1 (fl)2+(fljz y(C6) \ ,(%+fi & & , ( ° Note that the first terms in the square brackets vanish where evaluated at the top and bottom wall boundaries. The remaining integral can now be expressed with the order of integration restored as 126 Next, the two y-integral terms within G4 can be addressed as the separate components, GS and G6, as .5...§(._.._). 7 K K K K ' =J ZVksinaky 2V1 sina1y+ ZWksinaky 2W, sinaly dy O k=1 [=1 kzl [:1 7 K ' :2 (C8) = _ Vka +Wka 2k=l and 7 2 (a) ( 2:) l (295i — dy Gé(O[ & 7 K 2 K ' 2 = I (XI/k SinakY] +[2Wk sinaky] dy (C.9) 0 16:! k=l K 2 (Vk'z + (”V/(’2) k=l I NIV Each of the component integrals, G,, G, G3. Gs. and G, can then be directly computed from the Galerkin coefficients of the velocity and temperature. Re-assembling these components enables G to be computed with Simpson- rule integration across the z-dimension as 4(02 _ 30315.2 J27. 1 + T) R (To + T) 127 Appendix D: Galerkin Integrals for Entropy Generation Rate This section provides the details for the calculation of the generalized excess entropy production functional, ‘1’, which calculated by integrating over the entire domain, as indicated in equation (2.90). (H Ra a{vfl+wdéfl_é’2(fi)_§2(fl)] W _,. B m é’ & @2 @2 (.22)-? ...,,,...~ ‘11 5 HJ 0y >dy 0'22 0 oo (7%») 32(5)) 520») 22W) — (To + r) + B 6v @2 + @2 + @2 + 62-2 - (1- Mliil(”‘ 1))(6322 + 6w2) (DJ) L -. To compute ‘1’ from the Galerkin solution, the appropriate Galerkin expansions are substituted for the velocity and temperature factors, and the expression is simplified using the inner product coefficients. 128 We begin simplifying by first considering the conduction terms, and integrating by parts to reduce the order of the derivatives ‘7 52w) flan] (M of 22 M 1 7 2 7 l I 2 (MTV 6T) Mr) 0W) =l(”7.=.'l(7) dy 216%-.“(75 Ml '7 5T) 2 01”) 2 (D 2) -— —-—— + ddz ' (5)) .J .J 5 where the y-derivatives evaluated at the left and right boundaries vanish because of the adiabatic boundary conditions, and the the z-derivatives evaluated at the top and bottom walls cancel each other out because of the steady-state thermal energy balance. In a similar fashion, we address the z-component of the pressure term fifio + r)[6w%sfl]dydz = HUD + r)5wéPlg=l - (ii-[(76 + r)5w]6sz] dy (D.3) I7 = ”[TO'aw (TO + r)é(:v)]67’dydz 00 129 l and the z-derivatives in the Brinkman terms ((vowlavKi?) (3:2))de ’....(.(232)..(gfl)(‘0 \ l an, 015v) on») ‘ (Ema...) a} ] & dy (M 0L+[-—O§w+(T0+r) )] )1 ) ll [y [%&+ (To + r) 01:0] 6(3) .__ -Hl +[@&+(TO+ r)5(gv)]5(:)( afvdz 62' where the velocity variations evaluated at the boundaries all vanish because of the boundary conditions. Now ‘I’ may be re-expressed as H.R.5[.T [0(1) Law (a§_))2.(5§)2 T0 (SW 570— (TO + t)Ram&véT .4} -..,..).[..é<_522.:z<_&v_>-(.@)i(®_v>f) 0 6’ W“ > 0(3). .)- o 93, 4,2 52 & L5”) (MM J + 3(a) Mb a2 02 J (0.5) + (TO + r)(1— M|a|("‘ 1))(6vz + fiwz) K 1 130 To evaluate ‘P from the Galerkin solutions, the temperature variation, 6T, is assumed to be the higher modes (1 through K) of the expansion, which represent the fluctuations about the ”mean field” profile given by the zeroth mode. K éT—> Z Tk cosaky (D.6) k=I The velocity variations, 6v and (SW, are assumed to be the modes 1 through K of the corresponding expansions, representing the fluctuations about the quiescent state. K 512—) Z Vk sinak y (0.7) k=l K 6w—) 2 Wk sinak y (0.8) k=l Now the appropriate expansions are substituted for the temperature and velocity terms and fluctuations, and simplified. For computational 131 manageability, the integral in equation (D.5) is broken into component terms, S, S8. The first two are the heat transfer terms, defined as 7 5(a) (1.57)] S (2) a OT[v—+ w— d . l 5, 5.5 5 ' K K “‘ 7 K 2 Vk sinaky - Z 7} sina1y k=l [=1 = I){ZTmcosamy] K K 0 +[ZW15inalyMZTk'cosaky] [=1 k=l _) (0.9) K K K ' z Z Z le- Aklm(1)VleTm + Aklm(2)T/c WITm] m: and 52(2) (“fly + [fljz dy o2 k=1 7 K 2 K I 2 = I [Z-aka sinaky] +[Z Tk cosaky] dy (0.10) O k=l K ,2 Z (aszkz + Tk ) using the orthogonal property of the sines (3.14) and cosines (3.9). 132 Next we simplify the pressure terms 7 53(2) 5 (((TO + r)[6v 0 (76+r) ll K k=l -( K - TO'[ 2 Wm sinamy][ m=l 0151)) + 51566 03/ 32' Z Vk sinaky][ K m=l The buoyancy term is simplified as 7 K K S4(z) a ((TO + r)6w5fdy = “To + t)( 2 Wm sinamyJ (k2 Tk cosak =1 0 7 0 K K = (To + T)z Z Akaka k: 1m: 1 K 2 Wm, sinamy]( 2 Pk sinaky] J -—-6w61’} dy _5T__g & «W K 2 “ 011’1 Sinazy 1:1 k=l é K 2 Pk sinaky k=l m=l A [fig (To-i- T)akaPk + Z]: lAkak[( TO-I- r)",W + Ton ]( (D.11) (D.12) The Brinkman terms are simplified in components S5 and 8., defined as 32(fv)+ 6w 7 55(2) 5 I(T0 + T 0 V (T0+ 1' II OC—sV l (E: ')2-/(T0+T .-.)( (12)-(4:211 K 2— am2 Vm sin amy 133 212922 02 2 y] K K +[Zle sinaky]l(22=- am 2mV sinamy] 032 k: m- l K K I ~21Vk sinaky 'ZWk sinak k: k=l K .2 .2 )M[ (sz +sz )+V,, +Wk ( k=l (D.13) and 6v 6w S6(z) = (To [&é&—l+ fivflaz—qdy o r K K K K . = To I [2:1V Vk sinaky]( Z Vm sinaky]+ (2 Wk sinakyJ[ 2 Wm sinaky] dy O k m= I k=| m=l K (13.14) = £71121 Vka + Wka l The Darcy term is included in the S7 component, defined as 7 S7(z) a “To + r)(§v2 + 61v2)dy 0 7 K 2 K 2 = (76+ 7” [ng sinaky] +[Z Wk sinaky] dy (ID-15) o k- k=l K = :(TO+T)I;I(VI‘2 +Wk 2) And, finally, the Forchhiemer term appears in S8, given as 7 58(2) 2 ((TO + r)|ii|"' '(5122 + 61v2)dy 0 7 K K 2 K 2 = (To + t)(( Z Hk sinaky] (2 V, sinaly] + [Z W, sinaly] dy (D15) 0 k=l 1:1 1:1 K =(To+7)zlzl mil/41cm“) )Hk(VIVm +W1Wm) by again employing the nonlinear term expansion from equation (3.20), K 775 Ifiln-l -) 2 Hit sinaky k=l 134 From these simplified components, the generalized excess entropy production rate functional can be computed as l q! = HHBR(Sl + Sz)+ S3 - S4 + 3(- S5 + S6)+ S7 - M58]dz (0.18) 0 with the integration across the z-direction evaluated numerically using a Simpson's rule routine. 135 Bibliography Anderson, D. A., Tannehill, J. C., and Pletcher, R. H., 1984, Computational Fluid Mechanics and Heat Transfer, Taylor 8 Francis. Arnold, J. N., Catton, |., and Edwards, D. K., 1976, "Experimental Investigation Of Natural Convection In Inclined Rectangular Regions Of Differing Aspect Ratios,” J. Heat Transfer 98, 67-71. Bau, H. H., and Torrance, K. E., 1982, ”Low Rayleigh Number Thermal Convection in a Vertical Cylinder Filled with Porous Materials and Heated from Below,” J. Heat Transfer, 104, 166-172. Baytas, A. C., 2000, ”Entropy Generation for Natural Convection in an Inclined Porous Cavity," Int J. Heat and Mass Transfer, 43, 2089-2099. Beavers, G. S. and Joseph, D. D., 1967, ”Boundary Conditions at a Naturally Permeable Wall," J. Fluid Mechanics, part 1, 30, 197-207. Beck, J. L., 1972, ”Convection in a Box of Porous Material Saturated with Fluid,” Physics of Fluids, 15, No. 8, 1377-1383. Bejan, A., 1979, “A Study of Entropy Generation in Fundamental Convective Heat Transfer,” J. Heat Transfer 101, 718. Bejan, A., 1979, in S. Kakac, R. K. Shah, Aung (eds) Handbook of Single Phase Convective Heat Transfer, Ch 16. Wiley, New York. Bejan, A., 1982, Entropy Generation Through Heat and Fluid Flow, Wiley, New York. Bejan, A., 1984, Convection Heat Transfer, Wiley, New York. Bejan, A., 1987, ”Convective Heat Transfer in Porous Media," Handbook of Single-Phase Convective Heat Transfer, ed. Kakag, S., Shah, R., and Aung, W., Ch. 16, 16-1-34, Wiley, New York. Bejan, A., 1996, ”Entropy-generation minimization: The new thermodynamics of finite-size devices and finite-time processes," J. Applied Physics, 79, No. 3, 1191-1218. Bejan, A. and Almogbel, M., 2000, ”Constructal T-Shaped Fins," lntJ. Heat and Mass Transfer, 43, 2101-2115. 136 Bénard, H., 1900, ”Les Toubillons Cellulaires dans une Nappe Liquide,” Revue General des Science et Appliques, Vol. 12, 1261-1271, 1309-1328. Bories, S. A., and Combarnous, M. A., 1973, ”Natural convection in a sloping porous layer," J. Fluid Mechanics, 57, part 1, 63-79. Borkowska-Pawlak, B. and Korylewski, W., 1985, "Cell-Pattern Sensitivity to Box Configuration in a Saturated Porous Medium,” J. Fluid Mechanics, Vol 150, 169-181. Busse, F. H., and Joseph, D. D., 1972, ”Bounds for Heat Transport in a Porous Layer," J. Fluid Mechanics, vol. 54, 521-543. Caldwell, D. R., 1970, ”Non-Linear Effects in a Rayleigh-Bénard Experiment,” J. Fluid Mechanics, Vol 42, 161-175. Callen, H. B., 1960, Thermodynamics, Wiley, New York. Catton, I. and Edwards, D. K., 1967, ”Effect Of Side Walls On Natural Convection Between Horizontal Plates Heated From Below,” J. Heat Transfer 89, 295-299. Carton, |., 1970, ”Convection in a Closed Rectangular Region: The Onset Of Motion," J. Heat Transfer 92, 186-188. Catton, |.,1988, "Wavenumber Selection in Bénard Convection,” J. Heat Transfer, 110, 1154-1165. Chan, B. K. C., lvey, C. M., and Barry, J. M., 1970, ”Natural Convection in Enclosed Porous Media With Rectangular Boundaries,” J. Heat Transfer, 92, 21-27. Chan, Y. T. and Banerjee, S. 1981, ”Analysis of Transient Three-Dimensional Natural Convection in Porous Media,” J. Heat Transfer, 103, 242-248. Chandrasekhar S., 1961 , Hydrodynamic and Hydromagnetic Stability, Oxford, London. Cheng, P., 1978, ”Heat Transfer in Geothermal Systems,” Advances in Heat Transfer, Vol. 14, 1-105. Close, D. J., Symons, J. G., and White, R. F., 1985, ”Convective Heat Transfer in Shallow, Gas-filled Porous Media: Experimental Investigation," Int. J. Heat Mass Transfer, Vol. 28, No. 12, 2371-2378. 137 Combarnous, M. A., and Bories, S. A., 1975, ”Hydrothermal Convection in Saturated Porous Media," Advances in Hydroscience, Vol. 10, . 231-307. Darcy, H., 1856, Les Fontaines Publiques de la Villa de Dijon, Victor Dalmont, Paris. Davis, S. H., 1967, ”Convection in a Box: Linear Theory," J. Fluid Mechanics, 30, 465-478. Demirel, Y., Al-Ali, H. H., and Abu-AI-Saud, 1997, ”Entropy Generation of Convection Heat Transfer in an Asymmetrically Heated Packed Duct," Int. Comm. Heat and Mass Transfer, Vol. 24, No. 3, . 381-390. de Vahl Davis, G., 1968, ”Laminar natural convection in an enclosed rectangular cavity," Int. J. Heat Mass Transfer 11, 1675-1693. Elder, J. W., 1966, ”Numerical Experiments with Free Convection in 3 Vertical Slot," J. Fluid Mechanics 24, part 4, 823-843. Elder, J. W., 1967, ”Steady free convection in a porous medium heated from below," J. Fluid Mechanics, 27, part 1, 29-48. Elder, J. W., 1967, "Transient convection in a porous medium," J. Fluid Mechanics, 27, part 3, 609-623. Ergun, S., 1952, ”Fluid Flow through Packed Columns,” Chemical Engineering Progress, 48, No. 2, 89-96. Finlayson, B. A., 1968, ”T he Galerkin Method Applied to Convective Instability Problems," J. Fluid Mechanics, Vol 33, 201-208. Garcia, L., 1991, "A Theoretical and Computational Study of the Rayleigh- Benard Convection Problem in a Porous Medium Considering the Forchheimer Term and the Absence of Local Thermal Equilibrium," proposal for Ph. D. Dissertation, Department of Mechanical Engineering, Michigan State University, East Lansing MI. Genik, L. J., 1998, ”A Computational Approach to Simultaneous Two-dimensional Heat And Mass Transfer in 3 Heat Generating Porous Media," Ph. D. Dissertation, Department of Mechanical Engineering, Michigan State University, East Lansing MI. Georgiadis, J. G. and Catton, |., 1985, ”Free Convective Motion in an infinite Vertical Porous Slot: the Non-Darcian Regime,” Int. J. Heat Mass Transfer, Vol. 28, No. 12, 2389-2392. 138 Georgiadis, J. G. and Catton, l., 1986, ”Prandtl Number Effect on Benard Convection in Porous Media,” J. Heat Transfer, Vol.108, . 284-290. Georgiadis, J. G. and Catton, |., 1988, ”An Effective Equation Governing Conductive Transport in Porous Media," J. Heat Transfer, 110, 635-641. Georgiadis, J. G. and Catton, |., 1988, ”Dispersion in cellular thermal convection in porous layers," lntJ. Heat and Mass Transfer, 31, No. 5, 1081-1091. Georgiadis, J. G., 1988, ”Natural Convection in Porous Media with Anisotropic Dispersive Thermal Conductivity," Natural Convection in Enclosures - 1988, presented at the Winter Ann. Mtg. of the A. S. M. E. in Chicago. [OC 330.A43 1988] Glansdorff, P., and Progogine, |., 1971, Thermodynamic Theory of Structure, Stability, and Fluctuations, Wiley-lnterscience, London. Givler, R. C. and Altobelli, S. A., 1994, ”A Determination of the Effective Viscosity for the Brinkman-Forchheimer Flow Model,” J. Fluid Mechanics 258, 355-370. Gray, D. D. and Giorgini, A. 1976, "T he Validity of the Boussinesq Approximation for Liquids and Gases," Int. J. Heat Mass Transfer, 19, 545- 551. Holst, P. H. and Aziz, K., 1972, ”Transient Three-Dimensional Natural Convection in Confined Porous Media," lntJ. Heat and Mass Transfer, 15, 73-90. Home, R. N. and O'Sullivan, M. J., 1974, ”Oscillatory convection in a porous medium heated from below," J. Fluid Mechanics 66, part 2, 339-352. Home, R. N., 1978, ”Three-Dimensional Natural Convection in a Confined Porous Medium Heated from Below," ASME Paper number 78—HT-56, presented at the AlAA-ASME Thermophysics and Heat Transfer Conference in Palo Alto, CA, May 1978. Hsiao, K. and Advani, S. G., 1999, "A theory to describe heat transfer during laminar incompressible flow of a fluid in a periodic porous media," Physics of Fluids, 11, No. 7, 1738-1748. Hsu, C. T., Cheng, P., and Wong, K. W., 1994, "Modified Zehner-Schlunder Models for Stagnant Thermal Conductivity of Porous Media," Int. J. Heat Mass Transfer, 37, No. 17, 2751-2759. 139 Hsu, C. T., Cheng, P., and Wong, K. W., 1995, "A Lumped Parameter Model for Stagnant Thermal Conductivity of Spatially Periodic Porous Media," J. Heat Transfer 117, 264-269. Jiménez-lslas, H., L6pez-lsunza, F., and J. A. Ochoa-Tapia, 1999, "Natural convection in a cylindrical porous cavity with internal heat source: a numerical study with Brinkman-extended Darcy model," Int J. Heat and Mass Transfer, 42, 4186-4195. Jonsson, T. and Catton, |., 1987, "Prandtl Number Dependence of Natural Convection in Porous Media," J Heat Transfer, 109, 371-377. Joseph, D. D., 1976, Stability of Fluid Motions, Part II, Springer Tracts in Natural Philosophy, Vol 28, Springer-Verlag, Berlin. Katto, Y. and Masouka, T., 1967, ”Criterion for the Onset of Convective Flow in a Fluid in a Porous Medium," Int. J. Heat Mass Transfer 10, 297-309. Kaviany, M., 1984, ”Thermal Convective Instabilities in a Porous Medium," J. Heat Transfer 106, 137-142. Kaviany, M., 1984, "Onset of Thermal Convection in a saturated Porous Medium: Experiment and Analysis," Int. J. Heat and Mass Transfer, 27, No. 11, 2101-2110. Kaviany, M., 1991, Principles of Heat Transfer in Porous Media, Springer- Verlag, New York. Kladias, N. and Prasad, V., 1989, "Natural Convection in Horizontal Porous Layers: Effects of Darcy and Prandtl Numbers," J. Heat Transfer 111, No. 4, 926-935. Kladias, N. and Prasad, V., 1990, "Flow Transitions in Buoyancy-Induced Non-Darcy Convection in a Porous Medium Heated from Below," J. Heat Transfer 112, No. 3, 675-684. Kladias, N. and Prasad, V., 1991, "Experimental Verification of Darcy- Brinkman-Forchheimer Flow Model for Natural Convection in Porous Media,” J. Thermophysics and Heat Transfer, 5, No. 4, . 560-576, 1991. Kondepudi, D., and Prigogine, 1., Modern Thermodynamics: From Heat Engines to Dissipative Structures, Wiley, Chicester UK. Krishnamurti, R., 1968, "Finite Amplitude Convection with Changing Mean Temperature. Part 1. Theory," J. Fluid Mechanics, vol. 33, 445-455. 140 — ". . . Part 2. An Experimental Test of the Theory,” J. Fluid Mechanics, vol. 33, 457-463. Krishnamurti, R., 1970, "On the Transition to Turbulent Convection. Part 1. The Transition from Two- to Three-Dimensional Flow," J. Fluid Mechanics, Vol 42, 295-. — ". . . Part 2. The Transition to Time-Dependent Flow," J. Fluid Mechanics, Vol42, 309-. Lage, J. L., Bejan, A., and Georgiadis, J., 1991, "On the effect of the Prandtl Number on the Onset of Bénard Convection,” Int J. Heat and Fluid Flow 12, No. 2, 184-188. Lage, J. L., 1992, "Comparison Between the Forchheimer and the Convective Inertia Terms for Bénard Convection within a Fluid Saturated Porous Medium," Fundamentals of Heat Transfer in Porous Media, A. S. M. E. HTD-193, 49-55. Lai, F. C. and Kulacki, F. A., 1991, "Experimental study of Free and Mixed Convection in Horizontal Porous Layers Locally Heated from Below," Int. J. Heat Mass Transfer, 34, No. 2, 525-541. Lauriat, G. and Prasad, V., 1987, "Natural Convection in a Vertical Porous Cavity: a Numerical Study for Brinkman-Extended Darcy Formulation," J. Heat Transfer, 109, 688-696. Lapwood, E. R., 1948, "Convection of Fluid in a Porous Medium," Proc. Cambridge Philosophical Society, vol 44, . 508- . Luijkx, J. M. and Platten, J. K., 1982, "Precise Measurements of the Wavelength at the Onset of Rayleigh-Bénard Convection in a Long Rectangular Duct," Int. J. Heat Mass Transfer, Vol 25, No 8, 1252-1254. Malkus, W. V. R., and Veronis, G., 1958, "Finite Amplitude Cellular Convection," J. Fluid Mechanics, vol. 4, 225-. Masuoka, T., 1972, "Heat Transfer by Free Convection in a Porous Layer Heated from Below," Heat Transfer — Japanese Research, 1, No. 1, 39-45. Matheron, G. and de Marsily, G., 1980, "Is Transport in Porous Media Always Diffusive? A Counterexample," Water Resources Research, Vol. 16, No. 5, 901-917. 141 McDonough, J. M., 1980, The Rayleigh-Bénard Problem on a Horizontally Unbounded Domain: Determination of the Wavenumber of Convection, Ph. D. Dissertation, University of California, Los Angeles. McDonough, J. M. and Catton, l., 1982, "A Mixed Finite Difference-Galerkin Procedure for Two-Dimensional Convection in 3 Square Box," Int. J. Heat Mass Transfer, Vol 25, No 8, 1137-1146. Neale G. and Nader, W., 1974, "Practical Significance of Brinkman’s Extension of Darcy's Law: Coupled Parallel Flows within a Channel and a Bounding Porous Medium," Canadian J. Chemical Engineering, 52, 475-478. Nield, D. A., 1964, "Surface Tension and Buoyancy Effects in Cellular Convection,” J. Fluid Mechanics, Vol 19, 341-352. Nield, D. A., 1982, "Onset of Convection in a Porous Layer Saturated by an Ideal Gas" Int. J. Heat Mass Transfer, Vol. 25, No. 10, 1605-1606. Nield, D. A., 1983, "T he Boundary Correction for the Rayleigh-Darcy Problem: Limitations of the Brinkman Equation," J. Fluid Mechanics, Vol 128, 37-46. Nield, D. A., 1985, Corrigendium to "The Boundary Correction for the Rayleigh-Darcy Problem: Limitations of the Brinkman Equation," J. Fluid Mechanics, Vol 150, 503. Nield, D. A., 1991, "The Limitations of the Brinkman-Forchheimer Equation in Modeling Flow in 3 Saturated Porous Medium and at an Interface,” Int. J. Heat and Fluid Flow 12, No. 3, 269-272. Nield, D. A., 1995, Discussion of papers by K. Vafai and P. C. Huang, J. Heat Transfer, Vol. 117, 554-555. Nield, D. A. and Bejan, A., 1998, Convection in Porous Media, Springer- Verlag, New York. Palm, E., Weber, J. E., and Kvernvold, O., 1972, "On steady convection in a porous medium," J. Fluid Mechanics 54, part 1, 153-161. Pearson, J. R. A., 1958, "On Convection Cells Induced by Surface Tension," J. Fluid Mechanics, Vol 4, 489-500. Pedlosky, J. 1987, Geophysical Fluid Dynamics, 2nd ed., Springer-Verlag, New York. 142 Prasad, V. and Kulacki, F. A., 1987, "Natural Convection in Horizontal Porous Layers With Localized Heating From Below," J. Heat Transfer, 109, 795-798. Prasad, V. and Kulacki, F. A., 1984, "Convective Heat Transfer in a Rectangular Porous Cavity --- Effect of Aspect Ratio on Flow Structure and Heat Transfer," J. Heat Transfer, 106, 158-165. Prasad, V., Kulacki, F. A., and Keyhani, M., 1985, "Natural Convection in Porous Media," J. Fluid Mechanics, Vol 150, 89-119. Prasad, V., Kladias, N., Bandyopadhaya, A., and Tian, O.,1989, "Evaluation of Correlations for Stagnant Thermal Conductivity of Liquid-Saturated Porous Beds of Spheres," Int. J. Heat Mass Transfer, 32, No. 9, 1793-1796. Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P., 1992, Numerical Recipies in FORTRAN: The Art of Scientific Computing, 2"d ed, Cambridge, New York. Rayleigh, Lord, 1916, "On Convection Currents in 3 Horizontal Layer of Fluid when the Higher Temperature is on the Under Side," Philosophical Magazine, Vol. 6, No.32, 529-546. Rees, D. A., 1997, "The Effect of Inertia on the Onset of Mixed Convection in a Porous Layer Heated from Below," Int. Comm. Heat and Mass Transfer, Vol. 24, No. 2, 277-283. Roberts, P. H., 1966, "On Non-Linear Bénard Convection," in Non-Equilibrium Thermodynamics, Variational Techniques, and Stability: Proceedings of 3 Symposium held at the University of Chicago, May 77-19, 1965, Ed. by Donnelly, R., et. al., University of Chicago Press, Chicago. Roberts, P. H., 1967, "Convection in horizontal layers with internal heat generation. Theory," J. Fluid Mechanics, 30, part 1, 33-49. Roebers, J., 1986, "Computer Simulation of a Natural Convection Driven, Porous Media Chemical Reactor," Diplomarbeit, RWTH Aachen, Germany, work performed in Department of Mechanical Enginereing, Michigan State University, East Lansing MI. Rubin, H. 1975, "On the Analysis of Cellular Convection in Porous Media," Int J. Heat and Mass Transfer, 18, 1483-1486. Rudraiah, N., Veerappa, B., and Balachandra Rao, S. 1982, "Convection in a Fluid-Saturated Porous Layer with Non-Uniform Temperature Gradient," Int J. Heat and Mass Transfer 25, No. 8, 1147-1156. 143 Saatjian, E., ”Natural Convection in a Porous Layer Saturated with Compressible Ideal Gas," Int J. Heat and Mass Transfer, Vol. 23, 1681-1683. [also see response by Nield, 1982] San, J. Y., Worek, W. M., and Lavan, Z., 1987, "Entropy generation in combined heat and mass transfer," Int J. Heat and Mass Transfer, 30, No. 7, 1359-1369. San, J. Y., Worek, W. M., and Lavan, Z., 1987, "Entropy Generation in Convective Heat Transfer and Isothermal Convective Mass Transfer," J. Heat Transfer, 109, 647-652. Scanlon, J. W. and Segel, L. A., 1967, "Finite Amplitude Cellular Convection Induced by Surface Tension," J. Fluid Mechanics, 30, part 1, 149-162. Scriven, L. E. and Sternling, C. V., 1964, "On Cellular Convection Driven by Surface-Tension Gradients: Effects of Mean Surface Tension and Surface Viscosity," J. Fluid Mechanics, Vol 19, 321-340. Seki, N., Fukusako, S., and Inaba, H., 1978, "Heat Transfer in a Confined Rectangular Cavity Packed with Porous Media," Int. J. Heat Mass Transfer, 21, No. 7, 985-989. Singh, P., 1976, "T he Application of the Governing Principle of Dissipative Processes to Bénard Convection," Int J. Heat and Mass Transfer, 19, 581-588. Slattery, J. C., 1967, "Flow of Viscoelastic Fluids Through Porous Media,” AlChE Journal, 13, No. 6, 1066-1071. Slattery, J. C., 1969, "Single-Phase Flow Through Porous Media,” AlChE Journal, 15, No. 6, 866-872. Somerton, C. W., and Catton, l. 1982, "On the Thermal Instability of Superposed Porous and Fluid Layers," J. Heat Transfer, 104, 160-165. Somerton, C. W., 1982, Natural Convection and Boiling in Porous Media, Ph. D. Dissertation, University of California, Los Angeles. Somerton, C. W., McDonough, J. M., and Catton, l., 1982, "Natural Convection in Porous Media: A Mixed Finite Difference-Galerkin Solution with Wavenumber Predictions," Proc. Vll Int Heat Transfer Conference, paper NC40. A.S.M.E. Somerton, C. W., 1983, "T he Prandtl Number Effect in Porous Layer Convection," Applied Scientific Research, Vol 40, 333—344. 144 Somerton, C. W., and Jimenez, J., 1999, Independent Study Project, Dept. of Mechanical Engineering, Michigan State University, East Lansing. Stork, K., and Muller, U., 1972, "Convection in Boxes: Experiments," J. Fluid Mechanics, vol. 54, 599-611. Tian, Y. S., and Karayiannis, T. G., 2000, "Low Turbulence Natural Convection in an Air Filled Square Cavity, Part I: the Thermal and Fluid Flow Fields," Int J. Heat and Mass Transfer, 43, 849-866. Turner, J. S., 1968, "T he Behaviour of a Stable Salinity Gradient Heated from Below," J. Fluid Mechanics, vol. 33, 183-200. Vaifai, K., and Kim, S., 1995, "On the Limitations of the Brinkman-Forchheimer Extended Darcy Equation," Int. J. Heat and Fluid Flow 16, 11-15. Wang, C. Y., 1994, "Thermal Convective Instability of 3 Horizontal Saturated Porous Layer with a Segment of lnhomogeneity," Applied Scientific Research, 52, 147-160. Wang, C. Y., 1999 "Onset of Convection in a Fluid-saturated Rectangular Box, Bottom Heated by Constant Flux,” Physics of Fluids, Vol 11, No 6, 1673-1675. Weber, J. E., 1974, "Convection in a Porous Medium with Horizontal and Vertical Temperature Gradients," lntJ. Heat and Mass Transfer, 17, 241-248. Willis, G., Deardorff, J. W., and Somerville, R. C. J., 1972, "Roll-Diameter Dependence in Rayleigh Convection and Its Effect Upon the Heat Flux," J. Fluid Mechanics, Vol 54, 351-367. Wooding, R. A., 1957, "Steady state free convection of liquid in a saturated permeable medium," J. Fluid Mechanics, 2, 273-285. Wooding, R. A., 1963, "Convection in a saturated porous medium at large Rayleigh number or Péclet number," J. Fluid Mechanics, 15, part 4, 527-544. 145 111111111111111111