CONDENSATION IN A POROUS MEDIUM by LLOYD JAMES BRIDGE M.Math. (Mathematics) University of East Anglia, 1998 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA November 2001 © Lloyd James Bridge, 2001 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Abstract One-dimensional, steady-state heat and mass transfer with phase change are studied, for a two-phase zone in a water-saturated porous medium. A model problem for the saturation and fluid pressures is formulated, and an explicit temperature dependence for the saturation vapour pressure, together with an explicit saturation dependence for the capillary pressure allows us to solve the model problem for temperature and saturation profiles within the two-phase zone. A boundary layer analysis identifies existing models as approximations to the full model in the large vapour-pressure gradient limit. This approximation yields an uncoupled problem for saturation and temperature, and asymptotic analysis shows that a singularity is introduced if we accept the approximate model over the full model. An iterative method is described which allows both the full and outer models of the two-phase zone to be coupled to the two singlephase zones, and computations are performed with realistic control parameters for the entire three-zone system. Numerical solutions for both the two-phase zone and the three-zone system show excellent agreement between the full and outer formulations. ii Table of Contents Abstract ii Table of Contents iii . List of Figures v Acknowledgements vi Chapter 1. Introduction 1 Chapter 2. Modeling the 2-phase zone 2.1 Conservation of mass 2.2 Daxcy's Law 2.3 Energy Equation 2.4 Some empirical results 2.5 Boundary conditions 2.6 Nondimensionalisation 2.7 Typical data 5 6 7 8 8 11 12 13 Chapter 3. Asymptotic Analysis 3.1 Local analysis of full model at upper boundary 3.2 Local analysis of full model at lower boundary 3.3 A boundary layer problem 3.3.1 Outer solutions 3.3.2 Comparison with Udell's formulation, and physical interpretation 3.3.3 Inner solutions 15 17 19 20 21 23 24 Chapter 4. Numerical Solutions 4.1 The full coupled problem 4.1.1 Validation 4.2 Boundary layer problem 4.2.1 Outer problem in the limit 0 —> oo 4.2.2 Solving the outer problem alone 4.3 Regularisation 4.3.1 Error analysis 25 25 27 30 32 34 35 36 • Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem 5.1 Formulation of the problem 5.2 Solving the "outer" inverse problem 5.2.1 Results and discussion 5.3 Solving the inverse problem for the full model 5.3.1 Results and discussion, comparison with outer formulation iii 39 39 41 43 45 46 Table of Contents Chapter 6. Conclusions and Future Work Bibliography iv List of Figures 1.1 1.2 The PEM fuel cell Udell's two-phase zone problem 2 3 2.1 2.2 The three distinct zones Saturation pressure 5 10 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 Numerical solution to full problem Numerical and asymptotic solutions to full problem near the lower boundary. Numerical and asymptotic solutions to full problem near the upper boundary. Numerical solutions to full problem, varying heat flux Length of two-phase zone as a function of applied heat flux Pull and boundary layer solutions. Pull and boundary layer solutions near lower boundary. Pull and boundary layer solutions near upper boundary. Full and outer saturation as <f> —• co Convergence to constant temperature Outer solution and Udell's problem Error in computed length L using regularised boundary conditions for outer model Error in computed length L using regularised boundary conditions for full model Error in computed length L using regularised boundary conditions for outer model 26 28 28 29 29 30 31 32 5.1 5.2 5.3 5.4 5.5 5.6 5.7 Porous layer and three distinct zones Profiles for entire pack, with variation in To Profiles for entire pack, with variation in W Profiles for entire pack, large variation in heat Profiles for entire pack, for outer and full formulations Profiles for entire pack, for outer and full formulations Profiles for entire pack, for outer and full formulations v flux. 33 33 34 37 38 38 40 43 44 45 47 48 • 48 Acknowledgements I am grateful for the supervision and support of Dr. Michael Ward. This work was funded in part by MITACS, and I thank the MITACS Mathematical Modeling and Scientific Computation group, in particular Dr. Brian Wetton and Dr. Radu Bradean, for many helpful comments. vi Chapter 1 Introduction An understanding of heat transfer and fluid flow through porous media is central to research in many environmental and technological applications. Soil is one of many geological materials which is both porous, and permeable for the liquids and gases with which it naturally interacts. Waste disposal often requires fluid transport through a porous structure, as does oil recovery, where thermal effects are also significant. The design of porous insulation clearly raises issues of heat transfer and the migration of moisture. The modeling of these and other such processes must take into account the geometry of the porous solid, which impedes the flow of the fluid through the medium. The simplest physical problem of modeling heat and mass transfer for a single-phase fluid flowing through a porous medium is itself very complicated. However, many physical processes are more complicated, such as those involving a porous structure that contains a fluid undergoing a phase change. The drying of wood is a classic example considered by Whitaker [10], while fluid phase change in porous media is also at the centre of many technological problems. In configurations involving condensation and evaporation, regions of single-phase fluid and two-phase fluid often coexist within the porous medium. One technological advance inspired by environmental concerns is the development of the proton exchange membrane (PEM) fuel cell. In recent years, as automotive manufacturers have recognized the need for a low-emission alternative to the internal combustion engine, the design of fuel cells has attracted much interest, and demanded a mathematical model for the condensation of water vapour inside a porous electrode. Consider the simplified fuel cell configuration 1 Chapter 1. Introduction illustrated in Figure 1.1. The two electrodes indicated consist of a porous material. Hydrogen Fuel Fuel - Membrane Graphite Plates Gas Channels Figure 1.1: The PEM fuel cell. and oxygen diffuse through the electrodes to the membrane, where they react, producing water vapour, which is then transported back through the cathode. Condensation of this water vapour is observed to occur at certain locations within the cell, resulting in liquid water within the cathode, or emerging into the oxygen channel. The accumulation of liquid water in the gas flow channels impairs the efficient delivery of reactant gas to the electrode, whilefloodingof a porous electrode inhibits the diffusion of gas through the medium, thereby reducing the supply of reactant to the membrane. However, the membrane must remain hydrated for the reaction to occur, and the topic of "water management" draws on the effects of water in many inter-related processes inside the cell. The effects of water condensation can clearly be detrimental to the efficiency of operation of the fuel cell, and as such, a mathematical investigation of phase change in a porous medium is of interest to fuel cell manufacturers. In particular, a model problem amenable to computational study is desired. Bradean et al. [4] identify regions of water vapour oversaturation within a fuel cell electrode, where condensation is likely to occur, without including phase change effects in their model. In this work, we present a one-dimensional, steady-state model of heat and mass transfer in a porous layer, including phase change effects. Our aim is to form an understanding of this phase-change problem, by developing solution techniques which may subsequently be extended to a two-dimensional model of a porous electrode. 2 Chapter 1. Introduction A one-dimensional study by Udell [9] investigates phase change effects in a porous layer containing water. This layer is heated at the top and cooled from below. Experimental results indicate that, at steady state, there are three distinct zones within the porous pack: a vapour zone at the top, a liquid zone at the bottom, and a two-phase zone in between, through which there is a counter flow of liquid, driven upwards by capillary forces, and vapour, driven downwards by a pressure gradient. A model of the two-phase zone is presented, which assumes a constant temperature throughout this zone, and condensation and evaporation occurring at the lower and upper boundaries respectively, as illustrated in Figure 1.2. The model problem is solved to heated , \ . vapour liquid / cooled Figure 1.2: Udell's two-phase zone problem. give a saturation profile, which indicates the length of the two phase zone. A similar study is performed by Torrance [8], with heating from the bottom, and similar results are obtained. In this thesis, we present an extension to existing models of the two-phase zone. By adding an 3 Chapter 1. Introduction energy equation to the system considered by Udell [9], and assuming an explicit temperature dependence for the vapour pressure, we are able to solve not only for the saturation, but also for the temperature. In identifying a boundary layer near the lower boundary of the two-phase zone, we are able to show that Udell's formulation is an approximation to the full problem in the large vapour pressure gradient limit. Furthermore, the inclusion of temperature in our model allows us to compute the position of the two-phase zone, as well as the length of it. By coupling the boundary conditions for the two-phase zone to the simple saturation and temperature profiles in the two single-phase zones, we formulate a problem based on the continuity of temperature throughout all three zones, and obtain continuous saturation and temperature over the entire domain. The solution of the approximate model, which is valid only outside the boundary layer, is significantly simpler than solution of the full problem. However, results for both the full and approximate models are shown to be in good agreement, both for the two-phase zone problem alone, and for the entire three-zone system. 4 Chapter 2 Modeling the 2-phase zone Here we develop a one dimensional model for steady state heat and mass transfer within the two-phase zone of a porous layer, heated from the top, cooled from the bottom, and containing a specified mass of water. Consider the configuration shown in Figure 2.1. Suppose that the heated at top vapour only z=L 2-phase (liquid & vapour) z=0 liquid only bottom held at constant temperature Figure 2.1: The three distinct zones. 5 Chapter 2. Modeling the 2-phase zone heat flux through the porous pack maintains three distinct zones as shown. At the top, there is a zone containing only water vapour. In this top zone, phase change does not occur, and heat transfer is through conduction. Similarly, at the bottom, there is a conductive zone of only liquid water. In between, there is a zone where both liquid and vapour occur. In this middle zone, the liquid is driven upwards through the pores by capillary forces, and the vapour flows countercurrent to the liquid, driven by a slight pressure gradient. Phase change in our model is assumed to be possible anywhere in this two-phase zone, and one of the aims of this work is to find the length, L , of this two-phase zone. The model we present for the heat and mass transfer in the two-phase zone is based on those given by Bradean [3] and Whitaker [10]. This model incorporates the governing equations of mass conservation, energy conservation, and Darcy's Law for flow in a porous medium. 2.1 Conservation of mass We write mass conservation for each of the two phases separately. Firstly, consider the liquid, having material density pi, and moving with average velocity [7j through a medium of porosity (void volume fraction) e. The saturation s is defined to be the volume of liquid per volume of void space, and so the mass (per unit volume) in any control volume is espi. Hence we write conservation of mass for the liquid as | [cam] + ^ [espM] = F, (2.1) where the source term V is the condensation rate. Similarly, if the vapour has material density p v and moves with average velocity U over a control volume, then conservation of mass for the v vapour phase is ^Hl-8)p ] v + ^[e(l-s)p U ] v v = -T. (2.2) In (2.2), the term (1 — s) is the vapour volume fraction. Now, denote the discharge rate (volume flux per unit area) for the liquid and vapour as u\ and 6 Chapter 2. Modeling the 2-phase zone u respectively. Then we have v §i [espl] ih ' [plUl] = r + — [e(l - s)pv] + ^ [p u ] v = -T. v Assuming a steady-state condition, these equations reduce to d {pm dz + p Uv) 0, = v (2.3) and -T = j- {p u ). z v (2.4) v Further, working in one space dimension allows us to integrate equation (2.3), by considering the mass flux balance at either boundary of the two-phase zone. Since the fluid in each of the two conductive zones must be stationary, there must be zero mass flux across the boundaries of the two-phase zone. Thus, at both boundaries, the upwards liquid mass flux must be equal in magnitude to the downwards vapour mass flux. That is, piui = -p u . v v Thus, equation (2.3) integrates to give the overall mass balance piui + p u v v = 0. (2.5) 2.2 Darcy's Law Darcy's law gives relations between the discharge rates and the pressure gradient for each of the two phases as follows Ul = u = v « ( dpi 1 ([dz ~7~ + PW ) ) ( dp + + Pv9 Pv9 V\-rdz v Pv • (2-6) I-) 2 7 Here, K is the permeability of the porous material, n denotes the relative permeability of a r particular phase, p, is the dynamic viscosity, and p is the pressure. The body force term includes the acceleration due to gravity, g. Chapter 2. Modeling the 2-phase zone 2.3 Energy Equation The energy equation describing the heat transfer and phase change in the 2-phase zone is (piciui + p c u )^ v v = ^ v + hvapT- Here, c denotes specific heat capacity of a particular phase, h ( -8) 2 is the specific latent heat of vap water, K is the averaged thermal conductivity and T the temperature, which is the same for both phases at any given z. 2.4 Some empirical results Substituting for ui and u in equation (2.3), using (2.6) and (2.7), we obtain v (dpi d_ Pi — Krl ( -j Pi \dz dz \ r- Pl9 1 J H P p v K v ' dp r v I ~ v + V dz pg v = 0. (2.9) The mass balance (2.5) gives the liquid velocity in terms of the vapour velocity, and hence we may write the energy equation (2.8) as . ^ - C l ) p dT d ( dT\ T z T z \ T z ) , T> v U v = K + K a p _ V - Using (2.4) and (2.7) to write u and T in terms of the vapour pressure, p , we can write the v v equation above as KPv( l C ~ Cy)Kf V ( dp dT ^ v d ydT K— + az KPvKap fx -K v (dp — v r v \ dz + pg v (2-10) Equations (2.10) and (2.9) form a coupled system of two equations for the variables pi, p and v T, where the coefficients K i and n r rv will depend on the state of the system. We now make use of some empirical and analytical results to reduce the number of variables. Firstly, from Baggio et al. [2], we note that the vapour pressure in a porous medium may be written as -PcM Pv = Psat exp 8 I PiRT (2-11) Chapter 2. Modeling the 2-phase zone where p sat is the saturation vapour pressure, M the molar mass of water, and R the universal gas constant. From the data given in [2], the exponent in (2.11) has a typical order of magnitude of 10~ , and we shall hence take p = p 5 v for water vapour, p , sat sat in the analysis that follows. The saturation pressure is tabulated against temperature" in [1]. We note from this data the exponential-like appearance of the relationship. Assuming (approximately) a relationship of the type Pv = ae , bT we fit the data to an exponential curve to find the constants a and b. Fitting data for temperatures close to 370 K that are relevant to our model, we obtain the dimensional quantities a = 0.19743 Pa and b = 0.03525 K - 1 . (2.12) In Figure 2.2, we compare our exponential fit with the data given in [1], and note the closeness of the fit over the range of temperatures indicated near 370K. Secondly, the capillary pressure, p , is defined as the difference between the vapour and liquid c pressures, (2-13) Pc=Pv~Pi- As shown by Leverett [7], the capillary pressure is found to be a function of the saturation. The functional form for the capillary pressure, p = p (s) is known as the Leverett function. c c Udell [9] correlates this Leverett function to write the capillary pressure as p c = 5J(s) = 6 [1.417(1 - s) - 2.120(1 - s) + 1.263(1 - sf] . 2 (2.14) Here, ' = -(;)'• where o is the vapour-liquid interfacial tension, and e is the porosity of the porous medium. Also, we should note that the function J[s) given in 2.14 is an empirical relationship which 9 Chapter 2. Modeling the 2-phase zone saturation pressure as a function of temperature x10 280 300 320 340 Temperature T (K) 360 380 400 Figure 2.2: Saturation pressure. will depend on the particular porous medium being used. The correlation given in [9] is for a particular type of sand. However, in the absence of another model, and to allow comparison with our results, we will continue to use this model. Finally, the relative permeabilities are functions of the saturation. The relative permeabilities account for the decrease in mobility of one phase due to the presence of another. Thus, we expect n i to be zero when only liquid is present, and unity for only vapour. Similarly, n r rv should be zero when s = 1 and unity when s — 0. We will use the cubic forms suggested by Udell [9] to represent the results reported by Fatt and Klickoff [5]. 10 Chapter 2. Modeling the 2-phase zone In summary, we have the following empirical relationships: (2.15) r^rv p — (1 (2.16) ^) i = ae ] (2.17) bT v = P l ae -6J(s). (2.18) bT Substituting (2.15)-(2.18) into equations (2.9) and (2.10), we obtain two coupled equations for 8 and T: o, (2.19) and np (ci-c ) _ v d_ ,dT K~ + dz dz (i_ v 3 s) f d (_ Kp h v ] b T [ae vav 3 / d . b p— ^(l-sy[—[ae \dz } l 01 v \ dT g) + P v T z T + g Pv (2.20) 2.5 Boundary conditions The system (2.19), (2.20) is second order in both s and T, and thus requires four boundary conditions. At the top of the two-phase zone, that is, the boundary with the vapour only zone, the saturation (liquid volume fraction) is zero. The heat flux across the boundary, q, which we will initially specify, accounts for the temperature gradient and the evaporation of the upwards flowing liquid water at the boundary. Thus, two conditions at the upper boundary are = q- h piui, s= 0 vap at z = L. (2.21) Since we are considering a one-dimensional problem with no radial heat losses, the heat flux across the lower boundary must also be q. Here, the associated phase change is condensation of the downward flowing vapour, and the saturation is unity. Then the conditions at the lower boundary are frdT _ s= 1 11 at z = 0. (2.22) Chapter 2. Modeling the 2-phase zone The final coupled system for s = s(z) and T — T(z) is given by the second order equations (2.19) and (2.20), subject to the boundary conditions (2.22) and (2.21). 2.6 Nondimensionalisation We now nondimensionalise the system using the following scalings: z = —z*' and Pl9 L = —L*, Pl9 Pc = $P*c, T = T r*, ref k = kik\ KhyqppyPtg „ = Pi where * denotes a dimensionless quantity, and Equation _d_ (2.19) _d_becomes dz* dz* a bT T* +a 6 ref = 0, (2.23) and equation (2.20) becomes bT T* e ref ref Pv d_ dz* + oi T k^*^^— + dT* ~dz*~ ° P ( l -gffUJL \ t>TrefT (2.24) +a 6 dz* p \ 5 dz* I where a = p /pi and B — pi/p . We see from the right-hand side of equation (2.24) that one K p v h v e v v v suitable choice of T f for this problem is re Tref — KPvh pS va (2.25) PiKi With this choice, (2.24) becomes Ra(l - sf a d bTrefT* ~b~dz* + CX e _d_ dz* ,bT T* ref 6dz* dT* dz* +a (2.26) where Ra = KPv(ci -c )5 v PiKi defines the Rayleigh number for this problem. 12 (2.27) Chapter 2. Modeling the 2-phase zone 2.7 Typical data It is useful now to include typical data from Udell's experiment, so that we can identify the orders of magnitude associated with certain parameters and variables in our model. Table 2.1 lists typical values for dimensional and dimensionless quantities of interest, using data given by Udell [9] and Bradean [3]. One particular point of interest is the magnitude of the Rayleigh number. From the data given in Table 2.1, we see that Ra = O(10 ), while terms on the right hand side of equation (2.26) _1 are on the order of 10 . We have a small Rayleigh number problem. That is, diffusive effects 2 are dominant over convective effects. Thus, the convective left hand side of (2.26) is negligible, or equivalently, the energy equation (2.24) may be written as z> £r*Tref dT* K l K ~r IF K h Pv + „\3 I , vap ~^r {l ~ s ) a d + CX C, (2.28) [jd? where C is a constant to be determined. Though the Rayleigh number is a useful measure of the relative effects of diffusion and convection in the heat transfer problem, we will continue with equation (2.28) in preference to (2.26), since we can then choose the temperature scale T f conveniently such that T* re 1. Simply setting T f = 380K, a typical temperature inside re the two-phase zone, we guarantee T* w 1, which will prove to be convenient subsequently. 13 Chapter 2. Modeling the 2-phase zone Table 2.1: List of parameters and variables used. Symbol Interpretation Typical value Dimensional parameters K permeability 6.4 x 10" vapor density 1 Pv ' liquid density 10 Pi c specific heat of vapor 10 specific heat of liquid water 4.1 - 4.3 x 10 Q viscosity of water vapour 2.2 x 10~ Pv viscosity of liquid water 2.5 x 10" Pe k thermal conductivity of vapor saturated medium 1.0 Ki thermal conductivity of liquid saturated medium 2.5 heat flux 10 Q latent heat (water liquid-vapor) 2.5 x.10 hvap o surface tension (water liquid-vapor) 72.4 x 10~ R universal gas constant 8.31 M molar mass of water 18 x 10" acceleration due to gravity 9.8 9 a characteristic vapour pressure - see (2.12) 0.19743 b characteristic temperature scaling - see (2.12) 0.03525 H total height of porous pack 0.254 Tf reference temperature 380 Dimensional variables T temperature 360 - 390 vapor pressure 0.7 - 1.4 x 10 Pv capillary pressure 0 - 1 x 10 Pc Non-dimensional parameters e porosity 0.38 a 10" density ratio (p /pe) viscosity ratio (pt/p ) 11 (3 see (3.6) 13.4 <t> see (3.1) 13.5 V see (3.10) 7.3 D total height of porous pack 0.14 12 3 3 v 3 5 4 v 3 6 3 3 re 5 4 3 v v Non-dimensional variables s liquid volume fraction T temperature vapor pressure Pv capillary pressure Pc 0-1 0.95 - 1.05 3.98 - 8.13 0-0.6 Sources: Bradean [3], Udell [9]. 14 Units (SI) •y mr kg/m kg/m J/kgK J/kgK kg/ms kg/ms W/mK W/mK W/m? J/kg kg/s J/ mole K kg/mole m/s Pa Km K 3 3 2 1 K Pa Pa — — — — — — — — — — — Chapter 3 Asymptotic Analysis In this chapter, we derive asymptotic forms of the saturation and temperature profiles near the boundaries of the two-phase zone. Any singularities identified as a result of this analysis will highlight potential difficulties that the numerical solution of the problem will encounter. Further, the convenient temperature scaling described in the previous chapter will aid us in the comparison of magnitudes of terms in our equations, and thus facilitate the identification of regions within the domain where particular effects are negligible. First, let us formulate the problem as neatly as possible. Multiplying the energy equation (2.28) through by 6/KiT f, re we obtain where fi KiT f v re Again, the constant C is to be determined. We notice that (3.1) is simply a nondimensional version of - dT K——h h pu vap v v = constant. (3.2) The boundary condition (2.22) gives us that the constant in (3.2) is equal to the heat flux, q. In view of the overall mass balance, we note that the condition (2.21) is then automatically satisfied by equation (3.2). That is, neglecting the convective term in the energy equation means that one of the heat flux boundary conditions is trivially satisfied once the other one has 15 Chapter 3. Asymptotic Analysis been applied. We also see that the constant C is a scaled heat flux, and is given by • q. (3.3) KlTrefPW The problem now is reduced to solving the first order system " [? 3 * - - Tre/T dz* 8 1 , v \6 ) 3 {ah [? Tre/T 1 +°)=°' - (3 4) dz* subject to the conditions s(0) = 1, s(L*) — 0 for the unknown nondimensional two-phase zone length L*, and T* specified at one of the boundaries. We need to specify this extra condition on T* since the trivially satisfied heat flux condition leaves us with insufficient data to solve the system. In practice, we will later specify T* at a boundary of the two-phase zone by coupling the temperature profiles in all three distinct zones together such that the temperature is continuous throughout the entire porous pack. The analysis of the system is assisted by the identification of any large parameter associated with the problem. Defining the parameter 4> = bT , (3.6) ref we use the value of b from (2.12) and the value T f = 380K to calculate that <p « 13.4. Since re our choice of T f ensures that T* = 0(1), we have that e f * = e^ * is large. Physically, we bTre re T T have an exponentially large vapour pressure, and also an exponentially large vapour pressure gradient throughout the two-phase zone. Nonlinear exponential functions that can range over many orders of magnitude are common in combustion problems, and are associated with high activation energies in Arrhenius kinetics. A n asymptotic theory, referred to as high activation energy asymptotics, has been developed to treat such problems in combustion. Here, we take a similar approach to Kapha and Matkowsky [6]. We wish to have the variable exponential term in equation (3.5) being of order one. In order to ensure this, we factorise the exponential term as e 6T r e / T* = 4>T* e 16 = e * *(T--l)_ e ( 3 . 7 ) Chapter 3. Asymptotic Analysis Now the system (3.4), (3.5) may be written as (3.9) where (3.10) No terms have been neglected here, so we refer to this as the full problem. This nonlinear, coupled problem requires numerical solution, but we first investigate the asymptotic forms of the saturation and temperature profiles near the upper and lower boundaries. 3.1 Local analysis of full model at upper boundary In the proceeding work, we will drop the * for convenience, but we still work in nondimensional quantities. Near the upper boundary, we know that 5 —> 0 as z —» L. Also, the temperature T —> T , for top some constant T . Assume that, to leading order, s ~ By and T ~ T r top — Ey as z —> L, q top where B, D, q and r are positive constants, and y — L — z. Then, to leading order, dJ^ dz dJds ^ as dz = ^ Substituting for s and T in the system (3.8), (3.9), we obtain + a/3{l - 3By + W y r 2 2r - By ) {^e^ ^- ^{Eqy - ) - By ) {^e^ ^^<j>(Eqy - ) 3 3r Tt l) q 1 + a) = 0 , (3.11) + a) = C. (3.12) and KiEqyi- ) 1 + r?(l - 3By + W y r 2 2r 3 3r 17 Tto q x Chapter 3. Asymptotic Analysis In (3.11), we identify the powers of y as (0), (q - 1), (4r - 1), (3r), (r + q- 1), (2r + q - 1), and (3r + q - 1), (3.13) while (3.12) has the following exponents: (0), (q - 1), (r + q - 1), (2r + q - 1), and (3r + q - 1). (3.14) From (3.14), we see that the dominant balance in equation (3.12) as y —> 0 is given by q = 1. Then a consistent balance from (3.13) is clearly given by 0 = <?-l = 4 r - l . Thus, as z —* L, we have s ~ B(L - z ) / , 1 (3.15) 4 and T~T -E(L-z). (3.16) top To find the constants JB and £ we simply substitute these asymptotic forms back into equations (3.11) and (3.12). Retaining only the dominant terms, we have + + ) = o, a / 3 a and KE + r) (te'K »- ><l>E + a) = C. Tu r Thus, the constants are given by E= . d - V , a (3.17) and f ^ { C - K E ) V ' \ ( 3 . 1 8 ) It is clear from the local form (3.15) that s'(z) is singular at z — L, and as such, a finite difference solution to this problem will not be able to use the condition s(L) = 0. A regularised model will be discussed in the next chapter. 18 Chapter 3. Asymptotic Analysis 3.2 Local analysis of full model at lower boundary Near the lower boundary, we know that s —> 1 as z —> 0. Also, the temperature T —> Tb , for ot some constant Tb . Assume that, to leading order, s ~ 1 — Fz and T ~ Tj, + Gz as z —> 0, m n ot ot where £?, G, m and n are positive constants. Then, to leading order, dm^dJds dz ^ _ / ( 1 ) j F W n - i as dz Again, we substitute the asymptotic forms into the system (3.8), (3.9), and neglecting the (1 — s) term in (3.8), we obtain 3 & ?*x-V<p{Gnz - ) + J'(l)Fmz <t>{ n x m +1 =0 x (3.19) and K{Gnz - ) n 1 (^e^ - <t>{Gnz ~ ) + r,F z z Zm Tbot l) n l + oj=C. (3.20) In (3.19), we identify the powers of z as ( 0 ) , ( n - l ) , and (m - 1), (3.21) while (3.20) has the following exponents: (0),(n-l),(3m), and (3m + n - l ) . (3.22) Here, a consistent dominant balance is given by m — n — 1. That is, both the saturation and temperature profiles are locally linear near the lower boundary, and as z —* 0 we have the asymptotic forms s~l-Fz, (3.23) and T~T bot + Gz. (3.24) The constants F and G are found by substituting the dominant terms back into equations (3.19) and (3.20) to give ^ - <t>G Tbot l) + J'(1)F + 1 = 0, 19 Chapter 3. Asymptotic Analysis and KG = C. So finally, the constants F and G are given by C_ G (3.25) and -l-^e^^-^G (3.26) 3.3 A boundary layer problem The local solutions we have found near the upper and lower boundaries will be used to verify the numerical solution to the governing system. The full coupled system will be solved numerically, but in this section, we first present a simpler problem suggested by the identification of a boundary layer, in the asymptotic limit (j> —• oo. Consider again the coupled system (3.8), (3.9). Recalling that we have found 4> « 13.4 to be a "large" parameter, we see in the integrated energy equation / HT* dz* ~1 term 1 \ ,<Mr-i) C (3.27) 0(i) v_ term 2 that term 1 is small compared with term 2, provided that that is, (3.28) Now, using the data from Table 2.1, we compute f- w 7.3, 77 « 13.5, and K* = 0(1). We may write this conveniently as £ = 0(1) and rj = 0(0), 20 Chapter 3. Asymptotic Analysis and hence we see that (3.28) may be written as 1 -s > (f>- . (3.29) 2/3 The inequality (3.29) holds unless s ~ 1 - S i 0 / , where si ~ 0(1). That is, term 1 in -2 3 equation (3.27) is negligible except in a thin layer near z = 0, since s(0) = 1. The width of this layer then clearly tends to zero in the limit 0 —> oo. In the work which follows, we shall refer to this layer as the boundary layer. However, note that we are not presenting a classical boundary layer argument. The neglected term in a classical argument would typically be a small parameter multiplying the highest derivative, and then the width of the layer is accordingly a function of the small parameter. Here, we have in effect a small parameter multiplying a derivative which is of the same order as another term in the equation. While the result is a layer at the lower boundary in which the solution has a different form to the "outer solution" far from the boundary, the techniques of matching and width finding cannot be so easily employed in this case. 3.3.1 Outer solutions In the "outer" region, away from the lower boundary at z = 0, we have seen that term 1 in equation (3.27) is negligible, and thus dz* e 0(T*-l) a. r?(l - sf q (3.30) The neglect of the small term uncouples the system (3.8), (3.9). Substituting (3.30) in equation (3.8) we get '(*rV-s'H f + ( 3 ' 3 1 ) Equation (3.31) is a first order equation in s only, valid in the outer region. We can write this equation in the form 1 (!-<*) + -outer V~7 J / ( g , ) The solution to this outer problem for the saturation has a singularity at s = 0, in common with the full problem, but it also has a singularity at s = 1. Since J'(s) is nonzero at s = 0 21 Chapter 3. Asymptotic Analysis and s — 1, we see from equation (3.32) that the asymptotic form of the saturation profile at the boundaries is s~ 1-Biz / 1 s~ 4 B (L-z) l l A 2 as z^O, (3.33) asz^L. (3.34) Substituting these forms into the outer equation (3.31), and retaining only the dominant terms, we obtain 6 t / l , = 1 ] and J J'(0)B$ 4 | aBC ^ Q T] So the constants B\ and JE?2 are given by and Notice that the expression given for is consistent with neglecting K in (3.18), as we would expect. The saturation in the outer region will be found by the numerical solution of (3.32) subject to the "outer" boundary condition s(L) = 0. Once the saturation profile has been computed, the outer temperature profile follows from (3.30). Integrating, we find *(TM-D _ ^ - 1 ) = \ f ( _ ) Z e a UL \v[l's(t)f dt, J so that the outer temperature is given by Tauter (z) = 1 + T <P m I f Z 2 ( _ - )dt UL \v[l-s{t)f a / + e^^-D where the temperature at the upper boundary, T , must be specified top 22 (3.37) Chapter 3. Asymptotic Analysis 3.3.2 Comparison with Udell's formulation, and physical interpretation Here, it is interesting to look at another formulation. The saturation problem considered by Udell [9] reduces to solving the first order equation + S { Z } - a(3 (l-a)J'(s) ' ( 3 > 3 8 ) subject to s(0) = 1. Now, a = 0.001, and so we have almost perfect agreement with the outer problem (3.32). So Udell's solution corresponds to our outer solution, and so is not expected to be valid near the lower boundary. Consider the physical meaning of the outer solution. The neglect of term 1 in equation (3.27) corresponds to having no heat conduction through the two-phase zone, a common assumption in the engineering literature (eg. [9], [8]). Furthermore, in neglecting this term, we force, in dimensional variables, that puh v v vap = constant. Equivalently, —j-^ [pv^vhyapl ~ = 0, o, by equation (2.4). This says that the condensation rate is zero everywhere inside the domain, and thus phase change can only occur at the boundaries. So Udell's model clearly describes an evaporation front at the upper boundary and a condensation front at the lower boundary. Our boundary layer formulation similarly suggests that evaporation takes place in a front at the upper boundary, with no phase change in the outer region, but that condensation may occur in a layer near the lower boundary. Further insight is gained by writing equation (3.27) in the form I F < » - > * ( £ + *)= * t e r m 1 term 2 23 <-> 3 39 Chapter 3. Asymptotic Analysis We see that an 0(1) variation in vapour pressure corresponds to an 0(^) variation in temperature. The popular isothermal assumption in the literature (eg. [9], [8]) suggests that diffusive effects are negligible throughout the entire two-phase zone. However, we see that this is onlytrue in the limit <f> —• oo, and we have identified a layer near the lower boundary for which diffusion must be considered. 3.3.3 Inner solutions We have identified a layer near the lower boundary, where, to leading order, the saturation is s ~ 1. Thus, inside this layer, equation (3.8) reduces to d d (3.40) + The boundary layer system for the inner solution is thus given by equation (3.40) together with (3.27), where term 1 is now retained. The system is then easily decoupled. Firstly, we notice that equation (3.40) integrates to give f *(r-i) _ j ( ) e s + Z = A, (3.41) where A is a constant. If we denote the temperature at the bottom of the two phase zone by Tj, , then A is given by ot A = ^( ^*- )-J(l) T 1 = teW**- ), 1 (3.42) since J(l) = 0. So now we can write T in terms of s, T = +\ 1 l n A + J(s) - z (3.43) which is valid in the boundary layer. Now we can substitute for T and ^ in the second boundary layer equation (3.27), resulting in the following equation for s only: </> \A + + (l-sf(j'(s)^-l J(s)-z V +ay=C. (3.44) / Here, we have an equation of first order, which is to be solved numerically, subject to the initial condition s(0) = 1, and a specified temperature at z = 0. As in the case of the outer solution, the numerical solution for the saturation s is then used to compute the temperature in (3.43). 24 Chapter 4 Numerical Solutions In this chapter, we present numerical solutions to the full coupled model, as well as the boundary layer inner and outer problems described in Chapter 3. We have seen that for a <C 1, the outer solution is equivalent to the solution to Udell's problem [9], and so it will be useful to compare the features of the solution to the full model with those of the outer solution. The singular nature of solutions at the boundaries of the two-phase zone prompts a regularisation which is required to proceed with the finite difference computation, and we shall verify that the solutions agree with the corresponding asymptotic forms near the boundaries. 4.1 The full coupled problem The full problem (3.8), (3.9) is a first order system in two variables, and will be solved as an initial value problem subject to s and T specified at one boundary. We will employ an explicit 4th/5th order Runge-Kutta scheme, and first write the full problem in the form (4.1) From equation (3.9), we find that C - an(l - sf /2(*,T) = -5 K + ^0(1 - sfe^ - ) T 1 (4.2) and then substituting in equation (3.8) yields (6-kf )+s ( ^ e ^ - ^ + l) 3 h(s,T) = 2 s J'(s) 3 25 (4.3) Chapter 4. Numerical Solutions We then solve the problem (4.1)-(4.3) subject to the initial condition at either end. Suppose we choose the initial condition at the lower boundary, that is 1 •ot 2=0 where the temperature T bot (4.4) must be specified. Then, once the solution is computed, the unTbot=1,q=1000 1.004 0.14 Figure 4.1: Numerical solution to full problem. known length L is found from the second boundary condition s(L) = 0. Notice, though that since s'(z) —* co as z —• L, the condition s(L) = 0 is difficult to implement numerically. Rather, we must consider the condition at the upper boundary to be (4.5) '-top z=L where £ < 1 , with the temperature T top to be found. In Figure 4.1, we show the saturation and temperature solutions for the values of nondimensional temperature, T bot = 1, and dimensional heat flux, q = 1000, which determines C. Equation (3.3) and the data given in Table 2.1 give the value C = 1.89 x 10" q. 3 26 Chapter 4. Numerical Solutions Firstly, let us note the qualitative features of the two profiles. The gradients of both profiles appear steep but finite at z = 0, while the singularity in s(z) at z = L is apparent. The value L — 0.1248 is found from the computation. We report here that the saturation profile is similar to that shown by Udell [9], whose model is equivalent to our outer model. We will later present a comparison between results from the full model and results from the outer model. The temperature profile shown indicates that the isothermal assumption is not valid, and in particular, there is a sharp increase in temperature near the lower boundary. However, the change in temperature over the two-phase zone is small, on the order of 0.1%. In view of the scaling T f = 380K, the variation over the entire zone will only be on the order of one Kelvin, re whereas the difference in temperature over the entire porous pack ranged from 20K to 120K in Udell's experiments. In the absence of highly accurate measurements along the entire porous pack, the two-phase zone would clearly have appeared approximately isothermal. 4.1.1 Validation Firstly, let us check that the numerical solutions obtained approach the asymptotic forms at the boundaries, as described in Chapter 3. In Figures 4.2 and 4.3, we plot the numerical solutions near the lower and upper boundaries respectively, together with the asymptotic solutions given in Chapter 3, namely 8~1-Fz, T~T bot + Gz, . as2-» 0; as z L, (4-6) 8~B(L-zfl , A T~T -E(L-z), top where the constants B, E, F and G are calculated accordingly. Certainly, the full numerical solution is in good agreement with the local solutions near the boundaries of the domain. A further validation investigates the effects of changing the value of the heat flux q. Udell [9] reports that the length of the zone increases with decreasing heat flux. In Figure 4.4, we plot the saturation and temperature profiles for various values of q, keeping T t fixed. Indeed, the twobo phase zone length increases from L ss 0.06 when q = 3200 to L pa 0.14 when q = 800. We plot in Figure 4.5 the length L as a function of q, and note the smooth relationship. Interestingly, if we consider Udell's porous pack, which has a nondimensional length of 0.1411, we see that a heat 27 Chapter 4. Numerical Solutions comparison of asymptotic and numerical solutions near z=0 full numerical solution local solution 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 full numerical solution local solution 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01 Figure 4.2: Numerical and asymptotic solutions to full problem near the lower boundary. comparison of asymptotic and numerical solutions near z=L 1.0034 1.0034 1.0034 1.0034 full numerical solution local solution 0.122 0.12 0.123 0.1235 0.124 0.1245 0.125 0.1255 0.126 Figure 4.3: Numerical and asymptotic solutions to full problem near the upper boundary. 28 Chapter 4. Numerical Solutions Tbot=1, q varying 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Figure 4 . 4 : Numerical solutions to full problem, varying heat flux. dependence on two-phase zone length on heat flux 0.14 0.12r- 0.1 0.08 0.06 0.04 0.02 1000 2000 3000 _i 4000 i 5000 i 6000 i_ 7000 8000 9000 q (W/m ) 2 Figure 4 . 5 : Length of two-phase zone as a function of applied heat flux. 29 Chapter 4. Numerical Solutions flux of less than a critical value of approximately 600 W m - 2 gives L > 0.1411, and thus will not support the three distinct zones in a pack of this size together with the given temperature Tbot = 1- 4.2 Boundary layer problem comparison of numerical solutions for full, inner and outer problems 0.14 1.005 1.004 1.003 • 1.002 1.001 full problem outer solution inner solution 1 0.999 0.02 0.04 0.06 z 0.08 0.1 0.12 0.14 Figure 4.6: Full and boundary layer solutions. In Chapter 3, we derived a first order equation for the outer solution, and noted that this was a restatement of Udell's problem. A comparison between the inner and outer solutions and the full solutions near the respective boundaries should show close agreement. Further, the outer solution is expected to be valid throughout the domain, except only near the lower boundary. By considering the difference between the outer solution and the full solution on the whole domain, we can examine the effect of neglecting diffusion, and thus the error made by accepting Udell's formulation. To find the inner solution, we solve numerically the initial value problem for the saturation given by equations (3.42) and (3.44) with the initial conditions s(0) = 1, and T bot specified, again using a 4th/5th order Runge-Kutta method. The temperature T(z) is then trivially 30 Chapter 4. Numerical Solutions comparison of numerical solutions for full, inner and outer problems near z=0 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.016 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 Figure 4.7: Full and boundary layer solutions near lower boundary. computed using the algebraic equation (3.43). The outer solution is similarly obtained by solving the saturation problem (3.32) subject to the condition at z — 0, and then applying a numerical quadrature method to compute the temperature given by (3.37). We choose a midpoint method to deal with the improper integral as s —> 1. Notice that the length L is unknown in the "outer" formulation, and for comparison with the full solution, we simply set L to be that which is found from the full solution, and then compute the outer solution. In the Figures 4.6-4.8, we plot the full solution together with the inner and outer solutions for both the saturation and the temperature. Over the length of the two-phase zone, we can hardly see any difference between all three solutions, except near the lower boundary, where the outer solution looks poor. Zooming in on the two boundaries in Figures 4.7 and 4.8, we note the appropriate agreement of the full solution with only the outer and inner solutions and the upper and lower boundaries respectively, as desired. The close agreement between the solution to the full problem and the outer solution or inner solution at the upper or lower boundary respectively suggests that a matching method should result in a uniformly valid solution in good agreement with the full solution over the whole 31 I Chapter 4. Numerical Solutions comparison of numerical solutions for full, inner and outer problems near z=L 0.4 full problem outer solution inner solution 0.3 "> 0.2 0.1 0.12 0.121 0.122 0.123 0.125 0.126 1.0034 full problem outer solution inner solution 1.0034 1.0034 l - 1.0034 1.0034 1.0034 0.12 0.121 0.122 0.123 0.124 0.125 0.126 Figure 4.8: Pull and boundary layer solutions near upper boundary. domain. However, as previously mentioned, such matching would not follow the classical approach of van Dyke, and will not be explored here. Rather, we accept the outer problem as a good approximation to the full problem, and below, we will examine the behaviour of the outer problem in the limit —> oo. 4.2.1 Outer problem in the limit <>/ —• CXD In Section 3.3, we showed that the width of the boundary layer tends to zero in the limit <>/ —> oo, and thus the outer solution tends to the full solution. In Figure 4.9, we plot saturation profiles for <f> = 5, 10, 13.4 and 50, computed with q = 1000 and Tb — 1. We see that the outer solution ot quickly converges to the full solution as (j> becomes large. Also, the full solution converges, and the profile for our computed value (j) — 13.4 is in close agreement with the profile for <p = 50. We have also shown that as (j) —> oo, the temperature in the two-phase zone becomes isothermal. In Figure 4.10, we plot the temperature profiles obtained as <j> increases through large values. Our computations clearly show the convergence to an isothermal two-phase zone. 32 Chapter 4. Numerical Solutions saturation profiles in the large phi limit z Figure 4.9: Full and outer saturation as < > / —» oo. Temperature profiles for increasing phi 1.015 1.01 1.005 1 0 0.02 0.04 0.06 0.08 0.1 0.12 Figure 4.10: Convergence to constant temperature. 33 0.14 Chapter 4. Numerical Solutions 4.2.2 Solving the outer problem alone Here, since we have chosen to accept the outer model as a good approximation to the full model, we describe a method for computing the outer solution without prior knowledge of the length L , which we have so far taken from the full solution. Consider the problem of solving the outer equation 1 ~, /3 (1-s) *~ a l^" 3 (z) = outer 3 (4.7) J'(S) subject to s(L) = e, where £ < C l and L is unknown. Suppose we make an initial guess for the length, Li u, and solve backwards in z until s(z\) = 1 — e, as the solution is singular at s = 1. n Since equation (4.7) is autonomous, the solution is invariant under the translation z i—> z — z\, which would hence give a solution with s(L) = e and s(0) = 1 — e, as required by our outer model, where L = Li u — z\. We illustrate the idea in Figure 4.11, although equivalently, Udell n solves the problem forward in z. Solving the outer saturation problem 1.4 with initial guess translated Udell forward solution translate 0.16 Figure .4.11: Outer solution and Udell's problem. 34 Chapter 4. Numerical Solutions 4.3 Regularisation We have noted throughout this chapter that the singularity in the saturation gradient at the upper boundary of the two-phase zone in both the full and outer formulations, as well as at the lower boundary for the outer problem, is difficult to implement numerically. The regularisation method we have used so far for the outer problem has been simply to replace the boundary conditions s(0) = 1 and s(L) = 0, with s(0) = 1 — e and s(L) = e, where e C 1. For the full problem, we need only replace the condition for z = L. For later work in two dimensions, these regularised conditions on .s will be implemented. As such, it is interesting to calculate the error made in calculating the length L as e —> 0 . For this, we + seek an "exact" value for L, and so require a different solution method. Also, in the following chapter, we present an iterative method for computing saturation and temperature over the entire porous pack, for which we require an accurate value for L. Consider first the outer equation (4.7). For s € [0,1], the numerator is positive. Also, from the definition of J(s) in equation (2.14), we have that J'(s) = -1.417 + 4.24(1 - s) - 3.789(1 - s) . 2 (4.8) For s G [0,1], we see that J'(s) < 0. Hence, we deduce from equation (4.7) that for the outer model, s'(z) < 0 everywhere in the two-phase zone. Thus we have proven that s(z) is monotonic decreasing. Then it is valid to replace the outer equation (4.7), where s'(z) has singularities, with J' 35 (4.9) Chapter 4. Numerical Solutions If we now solve (4.9) for s e [0,1], with inital condition z(l) = 0, then z'(s) is zero at the boundaries and regular throughout the zone, and we compute for s = 0 and s = 1 exactly. Then the length L is given "exactly" from the condition z(0) = L. We can prove similarly that for the full model (4.1)-(4.3), the saturation, s(z) and the temperature, T(z), are monotonic decreasing and increasing respectively. So we may write that dz ds (4.10) h(s,Ty (£) and dT ds dTdz dz ds f2(s,T) (4.11) h(s,TY So we solve the system (4.12) hi*?) for s € [0,1], subject to the condition (4.13) 4.3.1 Error, analysis For the outer problem, if we choose s = 1 - e with e <C 1 as the lower boundary condition, then the asymptotic form (3.33) suggests that the error we make in z is on the order of (e/Bi) . 4 Similarly, from (3.34), the error we make in z at the upper boundary by choosing s = e as the boundary condition would be on the order of ( e / i ^ ) . We conjecture that if we compute the 4 solution for s G [e, 1 — e], then the error in L is approximately given by Lpxnrt — L £ « ( + ) £ For q = 1000, the crude approximation in (4.14) gives L 4 (4.14) 0. as e — L « 159e , and in Figure 4.12, 4 exact £ we plot the computed error and this conjectured error, showing reasonable agreement as e —> 0. Similarly, for the full model, the error we make in using the one regularised boundary condi36 Chapter 4. Numerical Solutions x 1 cf I e r r o r 1 i 1 Lexact-L when setting s(b)=1 -eps, s(L)=eps 1 1 1 1— 1 I Z 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 i — error 159* eps 0.04 4 0.045 eps Figure 4.12: Error in computed length L using regularised boundary conditions for outer model. tion s(L) = e can be approximated by assuming the asymptotic form (3.15). Then a crude approximation for the error in the computed value of L would be given by Lexact -L tt £ -^e as £ -> 0, 4 where B is computed from (3.17) and (3.18). Computing with q = 1000 and % approximation in (4.15) gives L (4.15) o t = 1, the — L sa 157e , and in Figure 4.13, we plot the computed 4 exact e error and this conjectured error, showing reasonable agreement as e —> 0 for the full model. For both the outer and full models, an error of less than 10~ is made in the computed length 4 is made for e chosen to be less than approximately 0.03. Finally, in Figure 4.14, we plot the computed error on a logarithmic scale, again illustrating the agreement with the conjectured error of order e . 4 37 Chapter 4. Numerical Solutions x10 error Lexact-L when setting s(L)=eps Figure 4.13: Error in computed length L using regularised boundary conditions for full model. error Lexact-L when setting s(0)=1-eps, s(L)=eps eps Figure 4.14: Error in computed length L using regularised boundary conditions for outer model. 38 Chapter 5 Coupling to the Single-Phase Zones: An Iterative Inverse Problem In this chapter, we formulate and solve a problem which has so far not been considered in the literature. Reports on the one-dimensional study of fluid phase change in a heated porous pack have noted the existence of distinct zones containing liquid only, vapour only and both phases. Typically, a model problem is described in the two-phase zone, from which the length of that zone may be found (see [9], [8]). However, the problem of finding the position of the two-phase zone in a three-zone problem has so far been neglected. The heat transfer in the two singlephase zones is by conduction only. Thus, the temperature in these zones is harmonic, and has a linear profile. The problem of finding the length and also position of a two-phase zone lying between two single-phase zones rests on the fact that the temperature is continuous throughout the entire pack. Thus, while the heat transfer problem in a single-phase zone with specified boundary temperatures is mathematically trivial, we must consider this problem for both the liquid and vapour zones, coupled with the more complicated problem for the two-phase zone which lies between them. 5.1 Formulation of the problem Consider the one-dimensional porous pack configuration shown in Figure 5.1. The height D is the height of the entire pack, which is a measurable constant for any particular experiment. For the one-dimensional, steady-state problem under consideration, we have so fax been able to compute the length, L, of the two-phase zone given the heat flux, q, across the boundaries 39 Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem 1 1 vapour zone D-L-X heat flux = q 'top 2-phase zone D 'bot liquid zone heat flux = q b Figure 5.1: Porous layer and three distinct zones. of that zone, and a boundary temperature, either T * or T^ . In an experiment such as that b ot op performed by Udell [9], however, the parameters which we are able to explicitly control are the temperatures applied to the top and bottom of the entire pack. Thus, we consider our temperature control parameters to be T\ and To, the temperatures at the top of the vapour zone and at the bottom of the liquid zone respectively. In the case of a porous electrode inside a fuel cell, the temperatures at the boundaries of the plate would be measurable, and we would wish to use these temperatures to predict the length and location of the three distinct zones. Conceivably, we could formulate an inverse problem such that, given the control parameters T\ and To, we can compute the values of the heat flux, q, through the vapour-only zone, and the temperature T * . Then the length, (D — L — X), of the vapour zone would be easily computed, t p since the temperature profile in that zone is linear. Since we would now have q and T * specified, t p we could compute the values of L and T * using our solution technique from Chapter 4, and b ot 40 Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem then X is known. However, with T * , To and X known, and a linear temperature profile in h ot the conductive liquid-only zone, we can calculate the heat flux qb t through this bottom zone. 0 Since we are considering a one-dimensional problem, there is no radial heat loss from the pack. Experimentally, this is achieved by insulating the walls of the container which holds the porous medium. Thus the downward heat flux across the vapour zone must equal the downward heat flux across the liquid zone, and so we have the requirement that Q = Qbot- (5-1) A further constraint must be incorporated into our model in order to meet this requirement, and so we seek a third parameter to introduce into the inverse problem. A third control parameter for the experiment is the total mass of water inside the porous pack. This parameter may obviously be varied for different experiments. The mass in the two conductive zones is trivially computed if their lengths are known, while the mass contained within the two-phase zone will depend on the saturation profile. Suppose the entire pack contains a measured mass (per unit cross sectional area) of water. Then this mass, W, must satisfy the integral constraint Xpi + f Jo L (spi + (1 - s)p ) dz + {D-Lv X)p v = W. (5.2) We now pose an inverse problem which, given inputs T\, TQ and W, will return all the unknowns shown in Figure 5.1, and be consistent with the heat flux condition (5.1). In solving this problem, we compute the saturation and temperature profiles within the two-phase zone, and we present the problem for both the "outer" equations (3.32), (3.37) and the full equations (4.1)(4.4). Further, given the inputs T\, TQ and W, the saturation and temperature profiles over the entire porous pack may then be compared for the outer and full models. 5.2 Solving the "outer" inverse problem We firstly consider the outer formulation, which gives a computationally simpler problem than the full model. Recall the uncoupled equations for saturation and temperature, derived in Chapter 3, namely SouteA ) = j7^) z 41 ' (.5-3) Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem and outer(z) = 1 + -jln (5.4) T 77(1 - s(t))3 where c=-^r Tn - - q (5 r 5) If we give the data T\, To and W, then the unknown parameters in Figure 5.1 which we require are the dimensional heat flux q and the nondimensional temperature T* . We first need to op initialize the value of q, to calculate C, in order to solve the saturation problem (5.3). Once we have s(z) and L, the integral constraint (5.2) gives v- _ W - (D - L)p - JQ (s v Pl + (1 - )p ) dz 8 v Pi- Pv Further, the dimensionless temperature at the top of the two-phase zone, T * , is given by t •'•top ii (D — L — X) * lengthscale p /T , (5-7) ref K v where lengthscale = 8/pig, as given in section 2.6. Once the temperature T * is known, we t p compute T(z) according to (5.4), hence finding the temperature at the bottom of the zone, T * . b ot Then the dimensional bottom temperature is simply T = T * * T f, and the computed heat bot b ot re flux through the liquid-only zone is given by £ 0 Tbot ~ TQ X * lengthscale' .J, g. Having specified the three control parameters T\, To and W, together with an initial value for the parameter q, we have used equation (5.7) to compute the unknown parameter T^, and (5.8) to find a second value for the heat flux. Now, in order to meet the heat flux requirement (5.1), we define the function /(?) = Q ~ Qbot, (5-9) and employ a quasi-Newton method, using centered differencing to approximate f'(q), to approximate the root of the equation f(q) = 0. 42 Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem The stopping criterion we use for the Newton iteration is that ||/(<z)|| < tol , where tol is a q q user specified tolerance on the heat flux. This criterion ensures that the requirement (5.1) is met to within some tolerance, and setting tol = 1 is sufficient to give only a 0.1% difference q between q and q t if q is on the order of 1000 W m " . To solve the outer equations in the 2 bo two-phase zone we use the numerical techniques and regularisation methods described in the previous chapter. The integral in equation (5.2) is computed using Simpson quadrature. It is worth noting that for this outer formulation, we are able to find the two parameters q and by implementing a single variable root finder. This is due to the fact that the outer problem for s(z) and T(z) may be uncoupled. The inverse problem for the full model, which we will describe in the next section, does not have this convenient property. saturation and temperature over entire pack for outer model, T1 =385, W=60 0 J L 0.05 0.1 0.15 0.2 0.25 z Figure 5.2: Profiles for entire pack, with variation in To. 5.2.1 R e s u l t s a n d discussion Our quasi-Newton code is seen to converge to within the specified tolerance, usually in less than 10 iterations. We take our initial value of heat flux q from the data from Udell's experiments [9]. 43 Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem saturation and temperature over entire pack for outer model, T1=385.1, T0=359.7 0 0.05 0.1 0.15 0.2 0.25 ! 1 1 z I C- , / • — • 0 0.05 0.1 0.15 W=55 W-60 W-62 q-1055.1886 q-1113.6908 q-1138.1157 0.2 0.25 z Figure 5.3: Profiles for entire pack, with variation in W. In Figure 5.2, we plot results in dimensional variables for three runs of our code. In each run, we keep the top temperature T\ and the water content W the same. The three distinct zones within the porous pack are clear from both the saturation and temperature profiles. Also, we can see clearly that the temperature variation over the two-phase zone is very small in comparison with the variations over the two single-phase zones, as expected. By way of validation, as we lower the bottom temperature To, we note that the heat flux increases, as expected. With this increase in heatflux,the decrease in the length of the two-phase zone is also apparent. The results of a further test are shown in Figure 5.3. Here we hold the temperatures T\ and To constant. Then the effect of an increase in.total water mass W is clearly evident, as the length of the liquid-only zone is seen to increase. Finally, in Figure 5.4, we choose the three control parameters in such a way that we can demonstrate the effects of large variations in heat flux. Here, we see again the decrease in length of two-phase zone with increasing heatflux,while the increase in liquid-zone length with increasing water mass is again clear. 44 Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem saturation and temperature over entire pack for outer model, varying heat flux 0.25 420 400 380 • 360 340 320 300 0 0.05 0.1 T1=385.1,TO=359.7,W=55,q=1055 T1 =385,TO=340,W=60,q=1283 T1=420,TO=300,W=70,q=2027 0.15 0.2 0.25 Figure 5.4: Profiles for entire pack, large variation in heat flux. 5.3 Solving the inverse problem for the full model For the full problem, we have usually specified the boundary conditions at the lower boundary of the two-phase zone, so now consider the problem shown in Figure 5.1, except with the heat flux through the liquid zone being equal to q. Now we will follow a similar approach to that described in section 5.2, and compute a new heat flux in the vapour-only zone, which we will label q . So we now have an inverse problem for the two parameters q and T * . Again, we top b ot initialize the value of q, and also initialize T * . We simply take the values we obtain from b ot solving the outer inverse problem. Now, to obtain s(z), T(z), L and T^, we solve the full coupled problem (4.1)-(4.4). The length of the liquid zone, X, is given by equation (5.6), and we have a second value for the temperature at the lower boundary of the two-phase zone, bot — q(X * lengthscale) Ki 45 +T 0 (5.10) Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem and also the heat flux through the top, vapour-only zone, given by kv(Ti-T *T ) ~ (D-L-X)* lengthscale' top °P qt ref Now we must meet the heat flux requirement q = q , together with the additional requirement top that T * = T * . We define the functions gi, g and Q by b ot b ot 2 (q, T * ) = (q- q )/1000, (5.12) g (q,T *ot) = T *ot-T * , (5.13) gi b 2 ot top b b b ot and and then employ a quasi-Newton method, using centered differencing for the Jacobian, to approximate roots of the vector equation G(q,T ) = 0. Here, we have defined the function gi bot so that it is of order one. The stopping criterion we use for the Newton iteration is that < tol, where tol is a user specified tolerance on the L2 norm of the error vector. The regularisation in the system (4.10)-(4.13) is used to ensure we compute to s = 0 exactly within the two-phase zone. 5.3.1 Results and discussion, comparison with outer formulation For good initial values of q and T , our code converges, again usually in less than 10 iterations. bot Further, the heat flux q and the temperature T bot converge to values close to those obtained from the outer formulation. In each of the Figures 5.5-5.7, we present the saturation and temperature profiles over the entire pack for both the outer and full models, given the same control parameters. The agreement is excellent. The computed heat flux for the two different formulations agrees to within 2% in all cases we have tried. The values of mean temperature across the two-phase zone for the full and outer models agree to within 1% in all cases computed. 46 Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem In conclusion, we note that the simplified "outer" formulation derived from a boundary layer argument not only approximates the two-phase zone results for the full model with close agreement, but also the results over the entire porous layer for the full model. Despite the extra singularity which we introduce into the problem by accepting the outer model, we have a much simpler problem to solve. The computational expense in solving the outer problem is significantly less than that for the full problem, in both the solution of the differential equations, and in the calculations for the inverse problem. saturationandtemperatureoverentirepack,T1=385.1, T0=359.7,W=62 outer model, q=1138.1157 full model, q=1152.91 0.8 - \ 0.6 0.4 0.2 - —— l I L_ 0 0.05 1 0.1 1 0.15 0.2 0.25 1 0 0.05 i 0.1 i 0.15 _i 0.2 ul 0.25 390 355 z Figure 5.5: Profiles for entire pack, for outer and full formulations. 47 Chapter 5. Coupling to the Single-Phase Zones: An Iterative Inverse Problem saturation and temperature over entire pack, T1=420, T0=300, W=70 0.25 outer model, q=2027.17 full model, q=2050.1105 0.1 0.15 0.25 Figure 5.6: Profiles for entire pack, for outer and full formulations. saturation and temperature over entire pack, T1=385, T0=340, W=60 1 1 1 1 — outer model, q=1283.2265 • full model, q=1307.6357 0.8 - 0.6 0.4 0.2 0 1 0.05 1 1 0.1 0.15 0.2 0.25 390 380 370 ' 360 outer model, q=1283.2265 full model, q=1307.6357 350 340 0.05 0.1 0.15 0.2 0.25 Figure 5.7: Profiles for entire pack, for outer and full formulations. 48 Chapter 6 Conclusions and Future Work A new formulation for one-dimensional, steady-state heat and mass transfer in a porous medium with phase change has been presented. We have extended the work of Udell [9] to allow for variations in temperature throughout the two-phase zone of a three zone system. The asymptotic analysis in Chapter 3 reveals a singularity in the saturation profile at the upper boundary of the two-phase zone, consistent with an evaporation front at this boundary. However, the saturation is found to be analytic at the lower boundary, and so the condensation does not necessarily take place at a front as previously assumed [9], but rather we allow for phase-change within the two-phase zone. We have identified a boundary layer near the lower boundary of the two-phase zone, in which phase change effects and heat diffusion are important, and an outer region, where these effects are negligible. The outer equations are consistent with a previous model [9], and approximate the full model in the limit of a large vapour pressure gradient. In Chapter 5, we have presented and solved a problem which has so far not been considered in the literature. We consider realistic control parameters for the entire three zone system, and have developed an iterative method whichfindsboth the length and the position of the twophase zone at steady-state. Heat transfer is considered separately for the liquid, two-phase, and vapour zones, and then the solutions for the three distinct zones are coupled together through continuity of temperature. The iterative inverse problem has been solved for both the full and outer formulations, with close agreement. 49 Chapter 6. Conclusions and Future Work Since the outer model is identified as a good approximation to the full model for both the two-phase zone calculations and those for the entire system, it is proposed that we accept the outer model in preference to the full model, as it is seen to be significantly less expensive computationally. Further work should include an asymptotic analysis of the matching between the inner and outer solutions, with the aim of replacing the full model with the outer model subject to an interface condition at the lower boundary. A new model of the entire three zone system will be presented in later work, using computational capturing techniques to find the boundaries of the three distinct zones. The new work in Chapter 5 of this thesis will serve as a reference for the evaluation of these capturing methods, which should be further extended to two-dimensional and three-dimensional models of a porous fuel cell electrode. Finally, the work in this thesis has considered the steady-state condensation problem. Later work will extend the models and methods described here to include time-dependence. Capturing the moving interfaces and boundary layers in a three zone system will be of particular interest in the design of a fuel cell, to aid the modeling of the start-up and shut-down of a cell. 50 Bibliography [1] M . M . Abbott and H. C. Van Ness. Theory and Problems of Thermodynamics (2nd edition), 1989, McGraw-Hill. [2] P. Baggio, C. Bonacina and B. A . Schrefler. Some considerations on modeling heat and mass transfer in porous media. Transport in Porous Media, 28, 1997, pp. 233-251. [3] R. Bradean. Heat transfer in porous media with phase change. Unpublished notes, 2000. [4] R. Bradean, K . Promislow and B. Wetton. Heat and mass transfer in porous fuel cell electrodes. Submission to International Symposium on Advances in Computational Heat Transfer, (Palm Cove, Queensland, Australia), 2001. [5] I. Fatt and W. A . Klikoff. Efffect of fractional wettability on multiphase flow through porous media. AIME technical note #2043. AIME Transactions, 216, 1959, pp. 246. [6] A.K. Kapila and B.J. Matkowsky. Reactive-diffuse systems with Arrhenius kinetics: multiple solutions, ignition and extinction. SIAM Journal of Applied Mathematics, 36, 1979, no. 2, pp. 373-389. [7] M . C. Leverett. Capillary behaviour in porous solids. AIME Transactions, 142, 1941, pp. 152-169. [8] K. E . Torrance. Phase-change heat transfer in porous media. Heat Transfer 1986, 1, 1986, pp. 181-188. [9] K.S. Udell. Heat transfer in porous media heated from above with evaporation, condensation, and capillary effects. Journal of Heat Transfer, 105, 1983, pp. 485-492. [10] S. Whitaker. Simultaneous heat, mass and momentum transfer in porous media: a theory of drying. Advances in Heat Transfer, 13, 1977, pp. 119-203. 51
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Condensation in a porous medium
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Condensation in a porous medium Bridge, Lloyd James 2001
pdf
Page Metadata
Item Metadata
Title | Condensation in a porous medium |
Creator |
Bridge, Lloyd James |
Date Issued | 2001 |
Description | One-dimensional, steady-state heat and mass transfer with phase change are studied, for a two-phase zone in a water-saturated porous medium. A model problem for the saturation and fluid pressures is formulated, and an explicit temperature dependence for the saturation vapour pressure, together with an explicit saturation dependence for the capillary pressure allows us to solve the model problem for temperature and saturation profiles within the two-phase zone. A boundary layer analysis identifies existing models as approximations to the full model in the large vapour-pressure gradient limit. This approximation yields an uncoupled problem for saturation and temperature, and asymptotic analysis shows that a singularity is introduced if we accept the approximate model over the full model. An iterative method is described which allows both the full and outer models of the two-phase zone to be coupled to the two single-phase zones, and computations are performed with realistic control parameters for the entire three-zone system. Numerical solutions for both the two-phase zone and the three-zone system show excellent agreement between the full and outer formulations. |
Extent | 3258078 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-08-12 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080047 |
URI | http://hdl.handle.net/2429/12035 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2002-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2002-0028.pdf [ 3.11MB ]
- Metadata
- JSON: 831-1.0080047.json
- JSON-LD: 831-1.0080047-ld.json
- RDF/XML (Pretty): 831-1.0080047-rdf.xml
- RDF/JSON: 831-1.0080047-rdf.json
- Turtle: 831-1.0080047-turtle.txt
- N-Triples: 831-1.0080047-rdf-ntriples.txt
- Original Record: 831-1.0080047-source.json
- Full Text
- 831-1.0080047-fulltext.txt
- Citation
- 831-1.0080047.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080047/manifest