UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Computational interface capturing methods for phase change in porous media Bridge, Lloyd James 2006

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2007-266878.pdf [ 9.32MB ]
Metadata
JSON: 831-1.0080046.json
JSON-LD: 831-1.0080046-ld.json
RDF/XML (Pretty): 831-1.0080046-rdf.xml
RDF/JSON: 831-1.0080046-rdf.json
Turtle: 831-1.0080046-turtle.txt
N-Triples: 831-1.0080046-rdf-ntriples.txt
Original Record: 831-1.0080046-source.json
Full Text
831-1.0080046-fulltext.txt
Citation
831-1.0080046.ris

Full Text

Computational Interface Capturing Methods for Phase Change in Porous Media by Lloyd James Bridge  M . M a t h . (Mathematics) University of East Anglia, 1998 M.Sc. (Applied Mathematics) University of British Columbia, 2002  A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF Doctor of Philosophy in The Faculty of Graduate Studies (Mathematics)  The University Of British Columbia November 7, 2006 © Lloyd James Bridge 2006  11  Abstract In this thesis, computational interface capturing methods for mathematical models related to fluid phase change processes in porous media are studied. The mathematical models are often singular and degenerate, which contributes to the computational difficulty. A n analysis of a smoothing method applied to a one dimensional free interface problem is presented. A n asymptotic analysis shows the dependence of the error in the computed interface location on the chosen small smoothing radius. Numerical convergence studies are performed for existing capturing methods applied to simple, scalar, moving interface problems, for later comparison with convergence rates for a new capturing method applied to a coupled, vector model problem. A model problem for two-phase fluid flow and heat transfer with phase change in a porous medium is described. The model is based on a steamwater mixture in sand. Under certain conditions, a two-phase zone, in which liquid and vapour coexist, is separated from a region of only vapour by an interface. T w o numerical methods are described for locating the interface in the one-dimensional, steady-state problem; one of these is based on an existing method, while the other uses the method of Residual Velocities. Agreement between solutions from these two methods is demonstrated, and the results from the steady-state computations are used as benchmarks for the numerical results for the transient problem. It is shown that methods such as front-tracking and the level-set method are not practical for the solution of the transient problem, due to the indeterminate nature of the interface velocity, in common with similar degenerate diffusion problems. A n interface-capturing method, based on a two-phase mixture formulation, is presented. A finite volume method is developed, and numerical results show evolution to the correct steady-state. Furthermore, similarity solutions are found, and the interface is shown to propagate at the correct velocity, by way of a numerical convergence study. Numerical results  Abstract for the two-dimensional problem are also presented.  iii  iv  Contents Abstract  ii  Contents  iv  List of Tables  vi  List of Figures  vii  Acknowledgements  x  1 Introduction and Background '.  1  1.1 1.2 1.3 1.4  Free and moving interface problems . Numerical methods for interface problems Two-phase flow with phase change in porous media . . . . . . Thesis overview  2 Solutions to prototype interface problems 2.1  2.2  2.3  The steady-state, two-phase Stefan problem 2.1.1 The one-dimensional problem 2.1.2 The two-dimensional problem 2.1.3 A smoothing method and asymptotic results 2.1.4 The method of residual velocities The time-dependent, two-phase Stefan problem 2.2.1 Mathematical formulation . 2.2.2 Analytical solutions 2.2.3 Formulation of the enthalpy method 2.2.4 Mathematical justification for the enthalpy method 2.2.5 Numerical results and convergence study. The Porous Medium Equation 2.3.1 Examples and applications 2.3.2 Analytical solutions  2 6 8 14  16 16 16 18 22 27 29 30 32 36 . . 42 45 49 49 51  Contents  2.3.3 2.3.4  v  Numerical methods Numerical results and convergence study  55 58  3 Nonisothermal, steady-state phase change in porous media 60 3.1 3.2 3.3  3.4  Mathematical formulation of the model problem A n iterative disjoint-domain method 3.2.1 Numerical results and discussion The method of residual velocities 3.3.1 Numerical results and discussion 3.3.2 A model problem in higher dimensions Benchmark solutions  61 68 71 75 76 76 • • • 79  4 The M 2 mixture method for the transient Udell problem . 80 4.1  4.2 4.3 4.4 4.5  Mathematical formulation of the model problem 4.1.1 Modified Stefan conditions at the interface z = L(t) 4.1.2 Model summary Fixed domain, mixture formulation Computational method Numerical results and discussion . Similarity solutions and numerical convergence study 4.5.1 A reduced model problem 4.5.2 Travelling wave solution for case p t(T) — aT 4.5.3 Numerical results and convergence study 4.5.4 Other choices for p t(T) Computations in higher dimensions 4.6.1 Mathematical model 4.6.2 Numerical results and discussion sa  sa  4.6  5 Conclusions and future work -5.1 5.2  Summary of results Future work 5.2.1 The "Big-H" regularisation  Bibliography  80 . . 82 86 87 89 93 95 96 99 104 107 107 108 109  117 117 118 ..119  121  vi  List of Tables 2.1 2.2  Constants for the Stefan problem . Errors for the enthalpy method, Forward Euler time-stepping, withe = 1 Errors for the enthalpy method, Forward Euler time-stepping, with c = 3 Errors for the enthalpy method, Forward Euler time-stepping, for a Neumann problem. Errors for Evje's method, Forward Euler time-stepping, for the Porous Medium Equation on —0.5 < x < 0.5, with \i — 0.1 and c — 3 • • Errors for Evje's method, Forward Euler time-stepping, for the Porous Medium Equation on —0.5 < x < 0.5, with JJL = 0.5 and c = 0.5 . • • • •'• •  49  3.1 3.2 3.3  Constants for Udell problem The effect of varying boundary temperature The effect of increasing mass  63 73 74  4.1  Errors for reduced Udell problem, implicit ( B E ) time-stepping with fi = 0.2, c = 4 Errors for reduced Udell problem, implicit ( B E ) time-stepping with fi = 0.2, c = 2 : Errors for reduced Udell problem, implicit ( B E ) time-stepping with // = 0.2, c = 1  2.3 2.4 2.5  2.6  4.2 4.3  45 47 48  59  59  105 106 106  VII  List of Figures 1.1  The temperature profile for the solution of the steady-state Stefan problem (1.3)-(1.4). Here, we have taken T = —1, T i = 1, K = 2.2, and K = 0.55 . Geothermal power generation . . . • A P E M fuel cell A three-zone system for Udell's experiment 0  ice  1.2 1.3 1.4 2.1 2.2 2.3 2.4 2.5 2.6  water  The steady-state, two-phase Stefan problem in two dimensions. Temperature profile and contours for steady-state two-phase Stefan problem. . . Temperature profile and contours for a steady-state two-phase Stefan problem with multiple interfaces. Smooth approximations converging to exact, nonsmooth steadystate temperature (dots) as e —* 0, for Example 2.1 '. Computed smooth approximation and asymptotic form, for Example 2.1. . A sketch of a monotonic increasing smoothed Heaviside function, H . Clearly, J ° H {£) d£ = 0(e), for X < 0. The steady-state, one-dimensional Stefan problem (2.1)-(2.2). Convergence to steady-state interface location, using Residual Velocity computations Heat balance at the freezing interface - the Stefan condition. . A typical Neumann solution to the freezing problem (2.32)(2.34) A typical travelling wave solution to the freezing problem (2.37). Enthalpy as a function of temperature Enthalpy and smoothed enthalpy as functions of temperature. The (z, £)-plane for the enthalpy formulation Evolution of numerical results using the enthalpy method. . . e  2.7 2.8 2.9 2.10 2.11 2.12 2.13 2.14 2.15  e  4 10 11 13 19 20 21 24 25 26 27 29 31 34 35 37 41 42 46  List of Figures  viii  2.16 Evolution of the interface using the enthalpy method, together with associated Neumann result. 47 2.17 Temperature history at the point z = 0.5 48 2.18 Computed interface location for travelling wave conditions, with c = 1 and c = 3 49 2.19 A thin droplet spreading under gravity over a solid substrate. 51 2.20 A solution of the steady-state Porous Medium Equation. . . . 52 2.21 A solution of the time-dependent Porous Medium Equation. . 53 2.22 Barenblatt-Pattle solutions of the Porous Medium Equation. . 55 2.23 Evolution of profiles given by conservative and non-conservative schemes applied to u — (u u ) 58 3  t  3.1 3.2 3.3 3.4  x  x  Udell's experiment A two-zone system for Udell's experiment. System for steady-state solution of the free interface problem. Solution profiles for free boundary problem using the iterative disjoint-domain method. Close-up of temperature profile for free boundary problem using the iterative disjoint-domain method Condensation rate in the two-phase zone. Convergence to correct interface location, using method of residual velocities, for T = 320, T i = 450, W = 30 Convergence to correct interface location, using method of residual velocities, for T = 375, T i = 500, W = 36 Steady-state system in higher dimensions  60 62 69  Mass conservation at the moving interface. Temperature profile and the moving interface G r i d and staggered grid Evolution to correct steady-state Interface location with grid refinement Temperature history at a point Timestep and interface location. Travelling wave profiles for mixture density, temperature and saturation 4.9 Computed and exact interface location for the reduced travelling wave problem 4.10 Steady-state., one-dimensional profiles. .  83 85 91 94 95 96 97  3.5 3.6 3.7  0  3.8  0  3.9 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8  72 73 74 76 77 78  103 104 109  List of Figures  ix  4.11 Saturation and temperature evolution. Saturation contours are i n the top row, temperature contours are in the bottom row. 110 4.12 Saturation and temperature profiles - initial condition and long time Ill 4.13 Saturation contours with liquid and vapour fluxes . 112 4.14 Saturation contours for an initial "blob" of two-phase fluid in the centre of the domain 113 4.15 Long time saturation and temperature profiles at fixed x for the spreading, migrating blob, together with one-dimensional, steady-state solution 114 4.16 Saturation contours for an initial condition with multiple twophase regions in the centre of the domain 115 4.17 Long time saturation and temperature profiles at fixed x for initial multiple two-phase zone problem, together with onedimensional, steady-state solution. 116  X  Acknowledgements This thesis was completed under the supervision of Brian Wetton. I wish to thank B r i a n for his support, guidance, kindness, and most of all, his encouragement. I am also grateful for the helpful comments of my supervisory committee, Roger Beckie, Chen Greif, Brian Seymour and Michael Ward. I have been lucky to work with members of the M I T A C S Mathematical Modeling and Scientific Computation group. In particular, I would like to thank Radu Bradean, Atife Caglar, Clarina Chan, W a n Chen, Roger Donaldson, A r i a n Novruzi, Keith Promislow, Akeel Shah, John Stockie, J F Williams and Jianying Zhang, for all their help and support. I am very grateful for the guidance and supervision given by Huaxiong Huang, during two summers spent at York University. The code for the twodimensional computations presented in this thesis was written during a very productive two weeks spent at the University of Saskatchewan. I wish to thank Ray Spiteri and his group there for their work and for making me feel so welcome. The support of friends and family has been so important, and I am forever grateful to them all. Special thanks go to Jason Mitchell and Johnson G o for their friendship, support, and rescuing skills during my time in Vancouver. This work was in part funded by M I T A C S and Ballard Power Systems. M y time in Saskatchewan was funded through the M I T A C S Mobility Fund.  /  1  Chapter 1 Introduction and Background Accurate numerical prediction tools and simulations are required for many physical processes which arise in industrial and environmental applications. The starting point for developing such numerical methods is the mathematical statement of the physical laws which must be obeyed. This typically involves one or more partial differential equations, which must be satisfied on the domain of interest. The mathematical model is completed by the specification of conditions on the boundaries of the domain and on any interfaces between regions inside the domain, and the initial state of the system. The solution of the mathematical model for its dependent variables forms the basis for a numerical simulation. It is very rare that a mathematical model for a complex physical process yields a closed form solution. A s such, a numerical approximation to the exact solution is usually sought. The development of an algorithm to find a good approximation relies in part on the careful implementation of the correct boundary and interface conditions. This can be a nontrivial task for a problem which is posed on a domain with fixed, known boundaries. Clearly, both the correct statement of a mathematical model, and the subsequent numerical procedure, become more complicated if there are moving boundaries or interfaces. The location of any free or moving boundary becomes an extra unknown in the model, and an extra condition is required. Free and moving boundary problems have attracted an enormous amount of interest in recent years. In this chapter, we give a brief introduction to free and moving boundary problems, and methods available for their numerical solution. We then discuss applications and model problems arising in the study of phase change, and in particular, two-phase flow with phase change in porous media, which is the focus of the work in this thesis. In this thesis, we formulate model problems for phase change in porous media, and develop numerical methods for the solution of such problems. These models and methods represent a contribution to the literature, in that they avoid many of the simplifying assumptions and computational regulari-  Chapter  1. Introduction  and  Background  2  sations that we will describe along the way, and we have demonstrated their validity by way of convergence studies using analytical solutions which we have constructed.  1.1  Free and moving interface problems  Many problems from applied science involve interfaces which separate domains in which different physical processes or flow regimes occur. The problem of a liquid droplet spreading under gravity concerns the interface between the liquid and the surrounding air, and the moving contact line which separates the wetted area beneath the droplet and the dry substrate ahead of it. Nonlinear hyperbolic equations represent gas dynamic quantities such as the gas velocity in a tube, where a region of moving gas and a region of stationary gas are separated by a moving interface, or shock. A melting block of ice contains both a region of liquid water, and a region of solid ice, and we conceive of a sharp interface between the two regions. We will refer to the problem of locating such an interface, which is stationary, as a free boundary problem or free interface problem. Correspondingly, we refer to the problem of finding a moving interface as a moving boundary problem, or moving interface problem. Many heat transfer and fluid dynamic processes result in free and moving boundary problems, and as such, these problems have generated much interest through industrial applications. A n alytical progress in the study of these complex industrial problems is often limited to asymptotic results and similarity solutions for idealised physical situations. A s such, the primary tool for solution of these nonlinear problems is usually numerical computation. However, the analysis of certain free and moving boundary problems has received attention, and is well understood. The prototype Stefan problem for phase change is well documented in the literature, in particular in [3, 20]. Further, the development of existing numerical methods for free and moving boundary problems has often relied on the mathematical insight into such problems. Classical solutions to free boundary problems are not usually available. In models where physical parameters are discontinuous across an interface, solutions will lack regularity at the interface. Consider, for example, the one-dimensional, steady-state, free interface problem for a partially melted block of ice. Suppose the ice-water system occupies a domain z € [0, D], and that an interface z = s lies between a region of ice (0 < z < s) and a region  Chapter 1. Introduction and Background  3  of liquid water (s < z < D). Assuming that the liquid is stationary, the temperature T throughout the system is given by the mathematical model: 0 0 < z < s, To ( < 0 ) ,  T(0)  T(s-)  \ (1.1)  °>  0  {.KwaterTz) z  T(s ) . T(D) +  s <z <D,  0, T (>0), x  together with the heat balance at the interface 8+  /  KT  .  K  ^waterT?  Z  T  -) =  0.  (1.2)  The problem of locating the phase-change interface, and determining the temperature throughout the domain, is commonly referred to as a Stefan problem, and so we shall call the system (1.3)-(1.4) the mathematical formulation of a S t e a d y - S t a t e S t e f a n P r o b l e m . Here, the known melting temperature is zero, and K is a thermal conductivity. The steady-state position, s, of the melting interface is given by s=  TK 0  T()Kice  i(  D.  T\K t wa  (1.3)  er  Then the temperature profile is given by 0<z<s, T(z)  (1.4) ^ (  Z  - D ) + T x  s<z<D.  In Figure 1.1, we show a temperature profile for the solution with T = —1, T i = 1, K = 2.2, K = 0.55, and D = 1. While the temperature is continuous throughout the entire domain 0 < z < D, it is clear that its derivative is not continuous across the interface z = s. For the corresponding time-dependent Stefan problem, the temperature profile may evolve with a discontinuous derivative at the interface. A s such, if we are able to reformulate such a problem as a P D E for temperature over the fixed domain 0 < z < D, then classical solutions will, i n general, be unavailable. That is, we will not be able to find solutions, with continuous derivatives up to 0  ice  water  Chapter  4  1. Introduction and Background Steady-state temperature profile for 10, two-phase Stefan problem  Figure 1.1: The temperature profile for the solution of the steady-state Stefan problem (1.3)-(1.4). Here, we have taken T = —1, T i = 1, K = 0  2.2, and  K er wat  ice  = 0.55.  the order of the P D E , which satisfy the P D E at all points in the domain. If the conditions under which a "solution" is defined are weakened, then we may admit solutions with discontinuous derivatives, as so-called weak, or generalised solutions. Such solutions need only satisfy an integral form of the original P D E . We will discuss this concept in detail i n Chapter 2, in our treatment of the so-called enthalpy method for numerical solution of phasechange problems. Weak formulations for free and moving boundary problems are presented in detail by Elliott and Ockendon [23], as well as others [3, 20]. In fact, it is the weak formulation of the Stefan problem which provides a basis for the mathematical justification of the enthalpy method. Another moving interface problem of physical interest may be described mathematically by the P o r o u s M e d i u m E q u a t i o n . This is a nonlinear diffusion equation whose diffusion coefficient is a power of the dependent variable. Mathematically, we have u = V . (u Vu), n  t  (1.5)  where n > 0, and t is time. Isothermal, ideal gas transport through porous media results in (1.5) with n = 1, where u represents the gas density (see, for example, [74]). The spreading of a thin liquid droplet on a solid substrate under the effect of gravity is described by (1.5) with n = 3 [23]. In this case,  Chapter 1. Introduction and Background  5  u is the height of the free surface of the droplet. The behaviour of solutions to equation (1.5) is significantly different to that of linear diffusion equations. A solution of the equation (1.5) is said to have compact support if u(x, t) vanishes outside a region of x , which varies with t. The region for which the solution is non-zero is called the support of u. In contrast with solutions of linear diffusion equations, equation (1.5) has the property that, for an initial condition it(x, 0) with compact support, the solution i/(x,t) also has compact support, for t > 0. The "interface", which is the boundary of the support of « ( x , t), will move with finite speed. Also, this speed may be zero until the solution at the interface has adjusted to a particular structure required for motion, giving so-called "waiting-time" solutions [32, 43]. Again, any interface-type solutions that we construct for (1.5) are understood to be weak, or generalised, solutions. T h e equation will not, in general, be satisfied by a constructed solution at the interface. Furthermore, singular gradients may appear at the interface, where the diffusion coefficient vanishes, and the equation degenerates. That is, equation (1.5) ceases to be parabolic at the boundary of its support. We will discuss some implications of the singularity and degeneracy of such problems further i n Chapter 2. There is an enormous number of other free and moving boundary problems which we will not describe here. The interface between liquid and the surrounding air i n Hele-Shaw flow [53], and the free liquid-air surface created by seepage of liquid through a porous dam are two well-studied examples, which bear some relation to our work. Crank's book [20] presents a variety of other free and moving boundary problems. However, our discussion will largely concentrate on the two prototype examples of Stefan problems and the Porous Medium Equation, which provide a good starting point for the mathematical detail of the interface problems discussed in this thesis. These problems may be extended, generalised and combined to model the more complex physical problems i n which we are particularly interested here. While the few analytical solutions which are available for free and moving boundary problems may only be applicable to a narrow range of physical problems, the insight gained from them is often valuable. Furthermore, the performance of the numerical methods which'are developed for the solution of such problems is typically evaluated by way of convergence study using test cases with analytical solutions.  Chapter 1. Introduction and Background  1.2  6  Numerical methods for interface problems  The formulation of an initial-boundary value problem requires the specification of boundary conditions. In a free or moving interface problem, one or more conditions are therefore required on the interface, whose position is an unknown in the problem. In particular, in a moving interface problem, at least one of these conditions will be on the interface velocity, and a successful numerical scheme must accurately capture this velocity. Front-tracking methods (see [20], C h . 4) are those which explicitly compute the interface velocity at each time step, and use this velocity to advance the interface. The use of fixed, uniform grids is not possible for such methods, since special difference formulae are required in the vicinity of the interface, as the interface will not, i n general, coincide conveniently with a grid point. Crank [20] describes various approaches to solving this problem, including the modified difference formulae, and adaptive time-stepping to ensure the coincidence of the interface with grid nodes. He also presents the related approach of frontfixing. Here, a coordinate transformation is performed after each update of the interface location, to ensure that the interface falls at a fixed coordinate value. Therefore, this method involves not only the implementation of the interface velocity at each time step, but also remeshing at each time step. Certain computational challenges become apparent when considering fronttracking and front-fixing methods. Both of these approaches require a numerical implementation of the interface velocity. In degenerate diffusion problems such as the porous medium equation, the interface velocity is often seen to be the limit of an indeterminate form (see [43], and [54] C h 7 ). Numerical methods are available for such scalar problems [25], which avoid the need for explicit implementation of the interface velocity, but in more complex problems to which these methods do not extend, front-tracking type methods would be impractical unless the indeterminate velocity could be accurately computed. Also, in greater than one dimension, the solution of boundary value problems on either side of the interface requires consideration of the interface geometry, and hence, unstructured grids and coordinate transformations are common. Donaldson [22] presents a front-tracking approach to a generalised Stefan problem using the finite element method in two dimensions. Representation of the interface as a spline contributes to a lengthy computation. For many Stefan-type problems, front-tracking methods have  Chapter  1. Introduction  and Background  7  been avoided. There has been much activity in recent years in developing methods for moving interface problems which avoid the need for explicit tracking of the interface. A method which computes the solution to a free or moving interface problem without explicitly computing the interface location, but rather recovers the interface location from the numerical solution to a reformulated problem, is generally referred to as a front-capturing method. The enthalpy method for the Stefan problem is described in detail by a number of authors [3,. 20, 23, 54], due to its wide applicability. The method takes a reformulation of the problem based on a transformation of the dependent variable from temperature to enthalpy, leaving a problem for the enthalpy and Kirchoff temperature over a fixed d o m a i n which contains the phasechange interface. The Stefan condition for the interface velocity is implicitly absorbed into the fixed domain formulation. Finite difference and finite volume schemes for the enthalpy method typically result in a lower order of accuracy than corresponding front-tracking schemes [24, 31], but the enthalpy method is widely used in practice, due to its ease of implementation (see, for example, [11, 24]). Since the method does not require any explicit handling of the interface, its extension from one to two or more dimensions is straightforward. In Chapter 2, we will describe the enthalpy method for the Stefan problem in detail. A n alternative method which has been developed more recently is the Level Set Method [55, 61]. This is also a capturing method which results from a reformulation over a fixed domain containing an interface. Rather than the interface velocity being implicitly absorbed into the new formulation, it appears explicitly in a Hamilton-Jacobi type equation. While the Level Set Method does not require a calculation of the interface location at each time step, as in front-tracking methods, it does require an accurate implementation of the interface speed function. In [18], Chen et al present a Level Set method for solving a Stefan problem. Their computation is more expensive than an enthalpy method for the same problem, but higher order of accuracy is achieved. The Immersed Interface Method is a modern numerical method which may be used for certain moving interface problems. Like the enthalpy method, it is a fixed grid method, but one which explicitly incorporates the interface conditions into the discrete scheme. Correction terms are added to the discretization in the neighbourhood of the interface, and interface conditions may be implemented by way of a smoothed delta function centered on the  Chapter 1. Introduction and Background  8  interface. L i [46] describes the method in some detail, and presents numerical convergence studies for several example problems. In particular, he demonstrates second order accuracy when the method is applied to a twophase Stefan problem. Careful choices for the correction terms and smoothed delta function can result in second order accuracy, in cases where the corresponding implementation of the enthalpy method can only achieve first order accuracy. Hou et al [36] present a hybrid method which combines the Immersed Interface and Level Set methods for interface problems, paying particular attention to achieving second order accuracy. However, for relative ease of implementation, the enthalpy method often appears to be the preferred method of solution for simple Stefan problems. Conservation of mass is an important concept in the development of numerical schemes for moving interface problems. It is well known that finite difference schemes for hyperbolic conservation laws with shock-type solutions must be derived in conservative form in order for the shock to propagate at the correct speed (see [44], for example). It is not enough for a scheme simply to be consistent with the governing P D E , when non-unique weak solutions exist. Critically, the conserved "mass" in the continuous physical problem must also be conserved in time in the discrete approximation. The necessity for discrete conservation is also a feature of numerical methods for equations of degenerate parabolic type [25, 26], which model moving interfaces which are not necessarily shocks.  1.3  Two-phase flow with phase change in porous media  A n understanding of heat transfer and fluid flow through porous media is central to the analysis of many environmental and technological processes. Soil is one of many geological materials that 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 may also be significant. The design of porous insulation clearly raises issues of heat transfer and the migration of moisture. The modelling 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.  Chapter 1. Introduction and Background  9  The problem of modelling heat and mass transfer of a single-phase fluid flowing through a porous medium is somewhat challenging. A significantly more complicated modelling problem concerns the heat and mass transfer, and phase change, of a fluid which flows through a porous structure. Twophase flow and transport with phase change in porous media has attracted much interest in recent years, from researchers in such diverse fields as microwave heating of foods and biomaterials [52], geothermal energy recovery [59, 74], and fuel cell technology [58, 21]. The process of wood drying is examined in detail by Whitaker [72]. In these examples, it is important to consider not only the flow and heat transfer, but also the effect of phase change, which further complicates the modelling and numerical effort. In configurations involving condensation and evaporation, regions of single-phase fluid and two-phase fluid often coexist within the porous medium. The location of interfaces between single-phase and two-phase zones are often of primary interest. Woods [73] analyzes the processes involved in liquid injection into a hot porous medium, as related to models of geothermal reservoirs. Typically, cool liquid water is pumped into a porous superheated reservoir below the Earth's surface. The injected water boils at a front, and the resulting water vapour is extracted by a well in order to drive turbines for the generation of electricity (see Figure 1.2). In [73], the flow of vapour ahead of the boiling front is analyzed, and also the flow of liquid behind the front. The rate of vapour recovery through the well depends on the rate of migration of the boiling front, which in turn depends on the liquid temperature and injection rate. While two-phase flow and phase change in geothermal reservoirs arise through the injection of cold liquid into a hot porous solid, similar flow phenomena arise in certain oil recovery processes. Typically, a displacing fluid is pumped into an underground reservoir, in order to extract oil from the porous rock. Allen et al [5] look in detail at the mathematical modelling of such flows, with particular emphasis on reservoir simulation. Peaceman [57] presents a comprehensive introduction to modelling and numerical methods for reservoir simulation, and Allen [4] reviews numerical methods for isothermal flows in natural porous media. A number of reductions to the mathematical model of multiphase flow in porous media are made under assumptions specific to reservoir flow. Flows are typically taken as isothermal, with no phase change throughout the reservoir, and capillary effects are assumed to be negligible. Under such assumptions, the fluids are incompressible, and  Chapter  1. Introduction  and Background  10  Turbine power station  Pump  cool liquid pumped down hot vapour extracted through well  hot porous reservoir  Figure 1.2: Geothermal power generation. the system of P D E s governing the flow decouples, and reduces to an elliptic equation for the total pressure, and a hyperbolic equation for the saturation of the wetting phase. These simplified problems may be solved numerically by a well established sequential Implicit Pressure, Explicit Saturation ( I M PES) time-stepping method (see [57]). Karlsen et al [38] take an alternative approach, by using a fast marching Level Set method for the saturation equation. The same assumptions of isothermal, incompressible flow without phase change are key to their formulation. Oil recovery processes can be thermally enhanced, whereby the injected fluid is a hot vapour such as steam. This acts to transfer heat to the oil, lowering its viscosity, which increases its mobility. Clearly, thermal and phase change effects become important in modelling such processes, as the injected steam will eventually condense at sites in the reservoir. Hanamura and K a viany [33] describe such a situation, where the injected steam condenses at a front, which then propagates through the porous medium. Bruining et al [16] formulate a model problem for steam injection into hot porous rocks, neglecting capillary effects, but including the effects of phase change, at a known phase change temperature. W i t h the thermal and phase change effects, the assumptions behind the I M P E S method are no longer valid, and alternative computational strategies must be devised. Furthermore, the modelling of two-phase flow and phase change in porous media on smaller scales is often  Chapter 1. Introduction and Background  11  d o m i n a t e d b y c a p i l l a r y effects, w h i c h alter the s t r u c t u r e of the m a t h e m a t i c a l m o d e l . I n p a r t i c u l a r , where c a p i l l a r y effects are i m p o r t a n t ,  saturation  equations are often of degenerate p a r a b o l i c type. T h e development of P r o t o n E x c h a n g e M e m b r a n e ( P E M ) fuel cells is a t e c h n o l o g i c a l advance p r o m p t e d b y e n v i r o n m e n t a l concerns. F u e l cell technology has recently received m u c h interest as the a u t o m o t i v e i n d u s t r y has recognized the need for l o w emission power supplies as a n alternative t o the i n t e r n a l c o m b u s t i o n engine. A s w i t h steam i n j e c t i o n processes i n o i l recovery, c o n d e n s a t i o n of s t e a m is also observed to o c c u r i n the porous electrodes of P r o t o n E x c h a n g e M e m b r a n e ( P E M ) fuel cells, a n d has received a t t e n t i o n from a n u m b e r of authors, for example, [12, 21, 56, 58, 70]. C o n s i d e r the  Graphite plates  Gas channels  Porous electrodes  F i g u r e 1.3: A P E M fuel cell. simplified fuel cell c o n f i g u r a t i o n i l l u s t r a t e d i n F i g u r e 1.3. T h e two electrodes i n d i c a t e d consist of a porous m a t e r i a l . H y d r o g e n a n d o x y g e n diffuse t h r o u g h the electrodes t o the m e m b r a n e , where they react, p r o d u c i n g water v a p o u r , w h i c h is t h e n t r a n s p o r t e d back t h r o u g h the cathode.  C o n d e n s a t i o n of this  water v a p o u r is observed t o o c c u r at c e r t a i n l o c a t i o n s w i t h i n the cell, resulti n g i n l i q u i d water w i t h i n the cathode, or emerging into the oxygen channel. T h e a c c u m u l a t i o n of l i q u i d water i n the gas flow channels i m p a i r s the efficient delivery of reactant gas t o the electrode, w h i l e flooding of a porous electrode i n h i b i t s the diffusion of gas t h r o u g h the m e d i u m , thereby r e d u c i n g the s u p p l y of reactant t o the m e m b r a n e . However, the m e m b r a n e must r e m a i n h y d r a t e d for the r e a c t i o n t o occur, a n d the t o p i c of "water management" draws o n the  Chapter  1. Introduction  and Background  12  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. Models of two-phase flow with phase change in porous fuel cell electrodes are seen to include thermal and capillary effects, resulting in degenerate parabolic equations as described above. A number of recent studies have used simplifications and regularisations of the model problems in order to deal with the numerical difficulties presented by their singular, degenerate nature. Bradean et al [12] identify regions of water vapour oversaturation within a dry fuel cell electrode, where condensation is likely to occur, but do not include phase change effects in their model. He et al [35] present a steady-state model problem, in which the liquid saturation dependence is removed from the capillary pressure and relative permeability terms. Therefore, the degeneracy and singularity in the problem as the liquid saturation vanishes is avoided. In [70], another steady-state model problem is described, for two-phase, multicomponent flow in porous electrodes, but neglecting thermal variations and phase change effects. Natarajan and Nguyen [51] solve numerically a time-dependent problem which includes phase-change effects. Saturation dependence in the relative permeability is included in their model, but regularised in order to simplify the computation. The same regularisation is used by Mazumder and Cole [50]. The effect of such regularisations is, in general, to smear out sharp interfaces. A time-dependent problem which makes no mention of any such regularisation, but which assumes isothermal conditions and does not include phase change, appears in [56]. A more recent study by Birgersson et al [10] considers the steady-state flow and phase-change with no apparent artificial regularisation added to the problem. In this thesis, we are primarily concerned with model problems for phase change in porous media which encompass all the mathematical difficulty of those appearing in the reservoir and fuel cell literature, but with a view to developing methods which may be applied in quite general, rather than specific, settings. W i t h this in mind, we note here that we shall only consider single component, two-phase flow, rather than the multicomponent flows described in the fuel cell literature, where oxygen, nitrogen and water may all coexist in the porous media. Multiphase, multicomponent computations will, in general, be more expensive than the computations which we describe here. However, the issues of singularity and degeneracy which we  Chapter 1. Introduction and Background  13  tackle here will not be further complicated by adding extra components to our model. More theoretical and general studies have been presented by other authors. A steady-state, one-dimensional study by Udell [64] investigates the effects on a sand pack, which contains water, of heating the layer at the top and cooling the layer from below. Experimental results indicate that, at steady state, there may be 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. The basic set up is shown in Figure 1.4. In the two-phase zone there is a counterflow of liquid, driven upwards by capillary forces, and vapour, driven downwards by a pressure gradient. Udell presents a model of the two-phase zone that assumes a constant temperature throughout this zone, with condensation and evaporation occurring at the lower and upper interfaces, respectively. The vapour in the two-phase zone is assumed to be fully saturated. The model problem is solved to give a saturation profile, which indicates the length of the two phase zone. A similar study is performed by Torrance [63], with heating from the bottom, and similar results are obtained.  tw>-phase  liquid  Figure 1.4: A three-zone system for Udell's experiment. The isothermal two-phase assumption made by Udell [64] appears to be popular throughout the literature. Temperature variation throughout a twophase zone may be critical in applications such as fuel cell design. A recent study by Wang and Wang [71] specifically examines the fuel cell setting  14  Chapter 1. Introduction and Background  under nonisothermal conditions, and shows numerical results for a two-phase zone. Also, a more complete steady-state model for Udell's problem [64] is presented i n [14]. Temperature effects are included i n this model, and phase change is allowed to occur throughout the two-phase zone. A numerical method is developed for locating interfaces between the single-phase and two-phase regions, and is described algorithmically in [15]. Now we consider the time-dependent Udell problem. Specifically, we consider the problem bf locating the interface between a single-phase region and a two-phase region in a single component, two-phase mixture, before the system has reached a steady-state. Such problems are of interest to fuel cell manufacturers when considering the effects of start-up and shut-down of cells. For the numerical solution of such a problem, fixed-domain, front capturing methods appeal. A method, based on mixture quantities, and using the same saturation assumption as Udell, is described by Wang and Beckermann [68], and implemented in [66]. Another model which has been used in the fuel cell literature penalizes any vapour not at saturation pressure into condensing at a large rate (see, for example, [21, 50]). Both of these methods require the solution of a fixed-domain problem, from which the interface location can be recovered. Computations using these methods have often been performed under the various simplifying assumptions which we have described above. Due to these assumptions, and the lack of analytical solutions to the coupled model problems, the challenge remains to show, either analytically or numerically, that such a computational capturing method for nonisothermal, two-phase flow with phase change in porous media yields an accurate solution in time. We aim to formulate model problems and develop capturing methods for their solution, for which we can show that the interface evolves with the correct location and velocity.  1.4  Thesis overview  ,  In this thesis, we formulate mathematical models for phase change i n porous media, and related model problems, and develop numerical methods for the solution of these problems, both steady-state and time-dependent. Throughout the thesis, we highlight points of particular mathematical interest. These are intended simply to aid the motivation for and understanding of our work, as well as to suggest directions of further interest. Most of our mathematical discussion will be illustrated using one-dimensional model problems. The  Chapter 1. Introduction and Background  15  emphasis throughout is on developing understandable, reproducible solution methods, rather than rigorous mathematical analysis. The remainder of this thesis is organised as follows. In Chapter 2, we discuss in some detail two prototype interface problems, which are related to the process of phase change in porous media. We describe formulations of the Stefan problem of phase change, and the Porous Medium Equation, and discuss numerical capturing and tracking methods to approximate their solutions. A new asymptotic analysis of commonly used smoothing strategies is presented for a smoothed, steady-state Stefan problem. Also, careful numerical convergence studies are presented, in preparation for a comparison with the results of Chapter 4. In Chapter 3, the steady-state, one-dimensional problem of phase change in porous media is described. This is an extension of an existing model, allowing for compressible vapour. Two numerical methods are described for the solution of this free interface problem, and numerical results are presented, showing good agreement between the two methods. These results are used as benchmark solutions for the time-dependent and two-dimensional computations in Chapter 4. In Chapter 4, the phase change in porous media problem is extended to include time-dependence, and reformulated as a fixed-domain problem. A numerical capturing method is developed for the solution of this problem, avoiding several simplifications which have commonly appeared in the literature. Computational solutions are shown to evolve to the correct steadystates predicted by the methods of Chapter 3. A n analytical solution is found for a reduced model problem, and numerical convergence studies using this exact solution show that solutions from our capturing method are indeed convergent. Furthermore, the convergence rates are comparable with the methods used for much simpler scalar problems, as described in Chapter 2. The implementation is extended to two dimensions, and computational results are shown. * In Chapter 5, we summarize our findings, describe ongoing work, and suggest directions for future work.  16  Chapter 2 Solutions to prototype interface problems In this chapter, we study two prototype free and moving interface problems of particular relevance, namely the two-phase Stefan problem and the Porous Medium Equation. We will formulate these problems mathematically, and discuss some analytical and numerical solution techniques for them. In particular, one-dimensional model problems are used for illustration. Much of this chapter is a review, intended as a mathematical background to accompany and guide the applied work in later chapters. We also present a new asymptotic analysis of smoothing strategies applied to interface problems, and suggest further work.  2.1 2.1.1  The steady-state, two-phase Stefan problem The one-dimensional problem  First of all, consider the steady-state, one-dimensional ice/water problem which we briefly discussed in Chapter 1 (see (1.1), ff.). We suppose that ice occupies a region 0 < z < s, and that liquid water occupies a region s < z < D. A t z = s, there is a melting/freezing interface which separates the two regions. In each of the two distinct regions, the temperature is governed by a steady-state heat equation. A t z = 0, suppose we have a given temperature, T , which is below the melting temperature. A t z = D, suppose we have a given temperature, T i , which is above the melting temperature. Also, take the melting temperature to be zero (This is for ease of illustration, and temperature is measured in degrees Celsius in the physical problems in this chapter. In later chapters, where phase change temperatures are not fixed, and the ideal gas law is important, we revert to temperature measured 0  Chapter 2. Solutions to prototype interface problems 17 in Kelvins). Then a mathematical model for the temperature throughout the system consists of two boundary value problems in each of the two regions, as follows: (K T ) = 0 0<z<s, T(0) = T ( < 0 ) , T(s-) = 0, (2.1) (K T ) = 0 s < z < D, T(s+) = 0, T{D) = T ( > 0 ) . ice  z  z  0  water  z  z  x  Now, since the interface position s is an unknown in the problem, we require one more condition. Energy should be conserved across the interface. In other words, the heat flux should be continuous across the interface, giving = 0.  KT,  (2.2)  It is a straightforward matter to solve the problem by calculating s using the three interface conditions, and then finding the temperature T(z) in the ice and the liquid water regions. Let us also consider an alternative formulation, to show how we may "capture" the interface from the solution of a transformed problem. We can reformulate the disjoint domain problem (2.1)-(2.2) as a linear problem for a transformed variable over the domain [0,D]. Let  v{T)= fK(0di,  (2.3)  JTo  where the conductivity K is given by K(T)  = |  ^K  T < 0  ice  Kwater  T > 0  or equivalently K(T)  = K  ice  + (K  -  water  K )H(T), ice  where H is the Heaviside function. In the context of the time-dependent problem, which we shall discuss in the next section, (2.3) is known as the K i r c h o f f t r a n s f o r m a t i o n . Then (2.1)-(2.2) becomes HT))  ZZ  = o,  v\ =o = 0, z  v\  (2.4) z=D  = v(Ti)  Chapter  18  2. Solutions to prototype interface problems  Trivially, this has solution (2.5) and the temperature T(z) can then be recovered using (2.3). The interface z = s is then "captured" using the condition that T(s) = 0. For the steadystate problem here, an exact solution is available. The steady-state position, s, of the freezing/melting interface is given by O-'Mee  (2.6)  D.  \ToKice ~ TiK t  J  wa er  Then the temperature profile is given by ( -^z T(z)  +T  0<z<s,  0  (2.7)  = {  [ ^ ( z - £ > ) + Ti  s<z<D.  It appears that the discontinuity in the conductivity function K at the interface does not hinder us analytically, or indeed, numerically.  2.1.2  The two-dimensional problem  Now let us consider how the problem is solved in two dimensions. Consider the problem shown in Figure 2.1. Liquid water and ice occupy regions Q i , respectively, and are separated by a freezing/melting interface T. We now simply extend the mathematical model (2.1)-(2.2) to two dimensions, giving V.(K VT) T(x,0) T V.(K VT) ice  water  T T{x,D ) T(0,y) 2  = = = =  0 (x,y) G f t T (x) (<0), 0 on r ~ , 0 (x,y)6ni, 2  0  = 0 onTV = T {x) (>0), = TuMv), T(D ,y)  \  (2-8)  x  2  = T (y), right  J  together with the heat balance across the interface: (KVT).n  = 0.  (2.9)  Chapter  2. Solutions to prototype interface problems  19  T = Tj  A  WATER \  r  ICE  0  0  A  T = To  Figure 2.1: The steady-state, two-phase Stefan problem in two dimensions. Here, n is the unit normal to the interface T. Once again, we may solve for the temperature and interface location by a change of variables which leaves a fixed-domain problem. Making the change of variables (2.3) v(T)=  f  K(Od£,  T  for a chosen reference temperature T f, the problem reduces to re  Av  = . 0,  v{x,0) .= v(x,D ) = v{0,y) = v(A.y) = 2  (x,y)e(0,A)x(0,A)  v(T (x)), v(Ti(x)), v{Tu t(y)), v{T (y)). 0  } (2.10)  f  right  That is, we may solve for the interface location and the temperature i n the disjoint regions by first solving Laplace's equation for v on the fixed domain. Given the solution v, the temperature T is easily recovered from the following: v+K T,  fsirsi ice  if  v  <  —K T f,  if  v  >  -K T .  ice  re  P- ice  ^fs^sl  * i water  ice  ref  (2-11)  Chapter  2. Solutions to prototype interface problems  20  Analytical solutions to the problem (2.10) may be found using Fourier series methods. For illustration here, consider a problem where we take To to be constant, using T = T , Di = D = 1, take Ti(x) to be 1-periodic, and give periodic conditions in the x-direction. Then the solution v of the fixed domain problem is given by ref  0  2  v(x, v) = —V + ^ I — — cos 2mrx + , , " — sin 2n7ra; ) sinh 2n7ry, \ ,»/ y Z^\ i h2n7r sinh2n7r J 2  n=l  s  x  n  (2.12)  where the Fourier coefficients are given by f  = ft viTtii))  a  = 2 f v(Ti(x)) cos2n7rx dx, n = 1,2,... = 2 JQ V(TI(X)) sin 2rwnr dx, n = 1, 2,...  n  b  dx,  (2.13)  Q  1  n  In Figure 2.2, we plot the steady-state solution to the problem (2.8), using Temperature profiles  Temperature contours  1  Figure 2.2: Temperature profile and contours for steady-state two-phase Stefan problem. periodic conditions in x, with K = 2.2, K = 0.55, T = —2 and T\(x) = 3 + 3cos27T2;. This two-dimensional problem is simply a sinusoidal ice  water  0  Chapter  21  2. Solutions to prototype interface problems  perturbation to the one-dimensional problem (2.1), with Xi = 3, which, by (2.6), has the interface at s = 8/11 ~ 0.73. This average interface position is clear in Figure 2.2. It is worth noting here that our solution method, namely solving the fixed-domain problem (2.10) for the transformed variable v, then recovering temperature using (2.11), is valid for more general problems than (2.8). Since v has continuous derivatives everywhere, the heat balance (2.9) will hold across any interface. The temperature recovery relies only on the value off, and does not require explicit knowledge of the location of any interfaces. Therefore, this method may be applied in general to problems where there are no interfaces, a single interface, or multiple interfaces. This idea extends to capturing methods for time-dependent problems, which gives these methods a particular appeal over front-tracking, as we shall see in the next section. In Figure 2.3, we plot the steady-state solution to a two-phase Stefan problem using periodic conditions in x, with K = 2.2, K = 0.55, X = —2 and T\{x) = 3 + 3.5cos2vra;. Clearly, there are two interfaces, which separate the liquid water from two regions containing ice - one near the cold interval along the upper boundary, and one near towards the lower boundary. ice  0  water  Temperature contours with two interfaces  0.6  »05M 0.41 0.31 0.21 0.11 o "  0  0.1  0.2  0.3  0.4  0.5  x  0.6  0.7  0.8  0.9  1  Figure 2.3: Temperature profile and contours for a steady-state two-phase Stefan problem with multiple interfaces.  Chapter  2.1.3  2. Solutions to prototype interface problems  22  A smoothing method and asymptotic results  The two-phase Stefan problem which we have described here has been illustrated by examples for which exact solutions are available. Such analytical solutions are not generally available for problems on irregular domains, problems with different boundary data, generalized Stefan problems with extra heat sources, for example, and for time dependent problems. In these cases, we inevitably resort to numerical solution methods. A n y code which involves the step of temperature recovery from the transformed variable will require "if" statements, and discretizations may be.hindered by the discontinuous temperature gradients. One approach to take is to smooth the discontinuous conductivity function K(T) over some radius in T , and solve the resulting smoothed problem. While some implementations of the enthalpy method for the transient problem, which we will describe in the next section, often make use of the exact conductivity, some degree of smoothing is often applied in order to aid computations (see, for example, [3, 20]). It appears that little analysis of such smoothing methods has been presented in the literature. Here, we present a brief analysis of a smoothing strategy applied to our one-dimensional, steady-state Stefan problem. We aim to give a regularization of the problem 2.1 by smoothing the conductivity function, leaving a smooth, fixed-domain problem. Let us write K (T) £  = K  + (K ^-k )H (T),  ice  wat  ice  (2.14)  e  where H is a C°° symmetric regularization of the Heaviside function, with smoothing radius e. Then the problem £  v(T) = (v(T))  zz  ff K (Od^ o  £  = 0,  ^U> = 0,  v\  (2-15) z=D  =  v(T ) 1  has a C°° solution T(z) which should tend to the exact solution of the Stefan problem as e —» 0. We do not prove this here, but observe this numerically. Of particular interest is the difference in the location of the interface where T = 0 between the exact and smooth solutions, and the dependence of this error on e. We see from (2.5) that the interface location z = s is given by * =^U  (2-16)  Chapter  2. Solutions to prototype interface  problems  23  Since the smooth solution satisfies v(T) = K (T  - To) + (K  ice  - K )  wateT  f  H {0  T  ice  d4,  e  (2.17)  JTQ  we see that the value of s for this smooth solution can be found approximately by considering the asymptotic evaluation of the integrals  and P H {£) d£,  f°H (Od(, £  e  JT  JTa  n  with T < 0, T i > 0 and e | T | , T i , for our choice of smoothed Heaviside H . We illustrate this for two common choices in the examples below. 0  0  e  Example 2.1 Let H (X)  = J + -tan- f-). 2 7T \ e '  (2.18)  1  e  Then  ' H (0dt T  =  e  To  hT-To)  + -I(T),  2  where  tan" 1 1  I(T) = T t a n " - - T 1  T  0  (2-19)  7T  - \ log  (J^) •  (2-20)  Recalling the identity tan  + tan- i  _i 1 -  A  v  f  X = |  - T^/ 2  X < 0 X  >  Q  ,  and the expansions tan- X = X  -  1  ^ + ...,  log(l + X ) = l - ^ + ...,  we find 1(0) ~ | T o + ( l + l o g | T | ) e - e l o g £ , 0  /(Ti) -  ^(Tt+To)-(loge. (2.21)  2. Solutions to prototype  Chapter  interface  problems  24  Hence, (0# J T  7T  [(l + l o g | r | ) £ - e l o g e ] , 0  0  and [ H (Odt  ~  Tl  £  ri-f^log^yje.  So we have -K T  v(0)  ice  0  +  E(Fe-eloge),  where E =  and  IT  F = l + log|T |, 0  and also ~  A-E\og-^-e,  where A  - KwoterEl  KiceTo-  Finally, substituting into (2.16), recalling the exact interface location (2.6), Smooth approximations as e -> 0, withTO=-3 T1=7 1^=0.8 k =0.25 N=80 wJBr  0  0.02  0.04  0.06  0.08  0.1  0.12  0.14  0.16  0.18  0.2  Figure 2.4: Smooth approximations converging to exact, nonsmooth steadystate temperature (dots) as e —• 0, for Example 2.1.  Chapter  2. Solutions to prototype  interface  problems  25  and noting that l/v(T\) ~ \ (l + f log j ^ e ) , we find the interface location to the smoothed problem is given by E S  ~  e — e log e D.  Sexact " t "  (2.22)  Clearly, s —* s act as e —> 0. In Figure 2.4, we show the convergence of the smooth approximations to the nonsmooth steady-state temperature. In Figure 2.5, we show a close-up of the computed smooth solution, showing good agreement with the asymptotic interface location (2.22). • ex  e=0.01 wlthT0=~3 T W 1^=0.8 * - „ = ° m  • O  0.11  25  N  =  4 0 0  exact solution computed smooth solution asymptotic interface position  0.112  0.114  0.116  0.118  Figure 2.5: Computed smooth approximation and asymptotic form, for E x ample 2.1.  Example 2.2 Let  tf (20 = - < | l + tanh e  §)}'  Here, we use the fact that log I cosh  el  ~ ^1*1 - log 2 e  for e < | X | ,  (2.23)  Chapter  2. Solutions to prototype  interface  26  problems  to obtain r-0  / fT«(Ode U  JIT To  ~  and  *  0  ci£ ~  Ti+o(e,eloge).  JT  JTo  0  Proceeding as in Example 1, we find that /  fe  _ r>  -water  \ice  Ii  1  KwaterTi — Ki To ce  , l?g2 I  (2.24)  T»£.  2  We notice the order of error in each of these examples, and conjecture that we cannot achieve a higher order error if our choice of smoothed Heaviside is monotonic increasing. The reason is due to the fact that the integral  JTo  for any monotonic increasing H , as suggested by Figure 2.6. £  He(X)  1  0.5  r •  0(1) X  >  Figure 2.6: A sketch of a monotonic increasing smoothed Heaviside function, H . Clearly, / ° H ( £ ) d £ = 0 ( e ) , for X < 0. e  e  The smoothing method presented here can, of course, be applied to higher dimensional problems, and time-dependent problems. We suggest further investigation of smoothing methods and a more formal analysis of our technique for the steady-state Stefan problem as open problems.  Chapter  2.1.4  2. Solutions to prototype interface  27  problems  The method of residual velocities  Recent work by Donaldson [22] presents a trial method for solving free interface problems, with a view to solving more complex problems for which fixed domain methods may not be straightforward to formulate or implement. For illustration of the method, particular attention is paid to the steady-state Stefan problem, i n two dimensions. Here, we describe the method for the one-dimensional problem. Consider again the problem (2.1)-(2.2), as shown in Figure 2.7. We can think of the problem as consisting of two boundary value problems, one on either side of the interface z ~ s. If the value of s is known, then we require only two conditions at the interface, but since s is an unknown i n the problem, we require three conditions there.  0 I '  T  (KiceT )z — 0 z  ICE  {KwaterT )z — 0 z  s \  '  Tbot  WATER  |  'KT ] _ +  Z  D  _  1  •  z  1  T — T  top  = O  (i)  T+ = 0 (ii) T~ = 0 (iii)  Figure 2.7: The steady-state, one-dimensional Stefan problem (2.1)-(2.2). The idea of the Residual Velocity method described i n [22] is to make an initial guess for s, then to solve the boundary value problems i n the ice and water regions, using two of the three interface conditions. T h e third condition, i n general, will not be satisfied, unless the interface is at the exact steady-state location. A n interface "velocity" is defined to be equal to the computed residual in this third condition, and the interface is subsequently moved according to this velocity and the chosen time-stepping scheme. This process is repeated until the interface approaches some limit for large "time". This is essentially a front-tracking method applied to steady-state problems on either side of the front. While the solution method treats the residuals as velocities, the evolution of the interface in "time" is completely artificial, and only the steady-state is meaningful.  Chapter  28  2. Solutions to prototype interface problems  For the one-dimensional problem in Figure 2.7, we now demonstrate the Residual Velocity method using each of the three conditions as velocities. Suppose the initial guess for the interface is s . 0  Example 2.3 Using the residual in condition (i) as the velocity (ie. solving the boundary value problems using conditions (ii) and (iii)), the interface velocity is given by . = - [ K dt \  w  a  t  e  r  - ^ - +K D —s  i  c  e  ^ \ , sJ  5(0)  =  a o  .  . (i)  Example 2.4 Using the residual in condition (ii) as the velocity (ie. solving the boundary value problems using conditions (i) and (iii)), the interface velocity is given by KiceTbot s — D  8(0) = s . 0  (h)  K, water  Using the residual in condition (iii) as the velocity (ie. solving the boundary value problems using conditions (i) and (ii)), the interface velocity is given by .. ds  I  m  KyjaterTfap  S  \  . .  .....  In Figure 2.8, we plot solutions of the O D E ' s from Examples 2.3-2.5. We take D = 1, T = - 1 , T = 1, K = 2.2, K = 0.55, and compute solutions given two different initial conditions, s = 0.5, and s = 0.97. In both cases, we see that all three velocity choices give an interface which evolves and converges to the correct steady-state solution s = 0.8. For the one-dimensional problem considered here, the analytical solution of the boundary value problems is trivial. In higher dimensions, where the interface geometry must be taken into account, solution of these problems requires numerical methods. Also, the choice of residual velocity is not bot  top  ice  water  0  0  Chapter  2. Solutions to prototype interface problems Steady-state interface location, using Residual Velocity method  11  1  1  1  1  1  1  1  1  29  r  "time" t  Figure 2.8: Convergence to steady-state interface location, using Residual Velocity computations. limited to the three shown here. The interface conditions used in the computations may be constructed using linear combinations of the three physical conditions (i)-(iii), and the interface should evolve to the same steady-state. In [22], a discretization of the interface for a two-dimensional problem is described, and the numerical properties of various choices of velocity applied with certain time-stepping schemes are analyzed.  2.2  The time-dependent, two-phase Stefan problem  In this section, we discuss the extension of the model two-phase Stefan problem to include time-dependence. W i t h a view to solving the phase change  Chapter  2. Solutions to prototype interface problems  30  problems in porous media described in later chapters, we look in particular at the formulation of the problem, exact solutions, and the capturing methods which have been developed for numerical solution. The contents of this section are largely drawn from well established work, which is particularly well described in detail by Crank [20] and Alexiades & Solomon [3].  2.2.1  Mathematical formulation  Let us consider the extension of the problem (2.1) to the time-dependent problem. Conceptually, we imagine a partially melted block of ice, with a moving freezing front, before a steady-state has been reached. Mathematically, we have pCi T T(0) T(s-) ce  t  = (K T ) = -T (< 0), = 0^ (K T ) = 0, = T (> 0). ice  +  0<z<s(t),  z  bot  water  T(s ) T{D)  z  z  s(t) < z < D, (  z  top  Here, p represents density, which we assume is the same in both liquid and solid phase, and c^water are the specific heat capacities of ice and liquid water respectively. Once again, we require an extra condition at the interface, which represents conservation of energy across the interface. Consider the diagram shown in Figure 2.9. Suppose that the interface is moving to the right, so that it is a freezing front. In small time St, the front moves a small distance Sz. During this time interval, the heat which flows out (per unit area) into the ice region is approximately K T St. This must equal the heat which flows in from the water region, plus the heat released upon freezing the mass (per unit area) pSz. Thus we have ice  z  where L is the specific latent heat of freezing. Taking the limit Sz, St —• 0, we find (2.26) pLs(t) = KT l water Z  J  ice  Chapter  2. Solutions to prototype interface problems  31  WATER  Figure 2.9: Heat balance at the freezing interface - the Stefan condition. where s is the interface velocity. This condition for the interface velocity is generally referred to as the S t e f a n c o n d i t i o n . It is worth noting that we will arrive at this same condition regardless of whether we consider a melting or freezing process. In order to complete the specification of the time-dependent problem, we also add an initial temperature distribution, and an initial interface location. That is, we also give T(z,0) = T (z),  0<z<D,  init  s(0) = s . 0  (2.27)  Now, defining Ot-ice. —  Kice  Kwater &water — PCwater  ce  our two-phase Stefan problem becomes T T(0) T(s-) T T(s+) T(D)  = a T 0<z.<s(t), = T (< 0), = 0, = a T s{t) < z <D, = 0, = T (> 0),  pLs(t)  = - KT  T(z,0) s(0)  = T (z), =s .  t  ice  zz  bot  t  water  zz  top  r ^ L  init  0  -j water Z  J ice  0<z<D,  (2.28)  Chapter  2.2.2  32  2. Solutions to prototype interface problems  Analytical solutions  Exact, analytical solutions to two-phase Stefan problems such as (2.28) are not generally available. The two most common types of solution, namely the Neumann similarity solution and the travelling wave solution, appear to have limited application in real world problems, but such solutions are valuable in giving mathematical insights into the Stefan problem, providing useful approximations, and evaluating the performance of numerical schemes. Here, with a view to constructing analytical solutions to the moving interface problems which we will encounter in later chapters, we describe both a Neumann and travelling wave solution.  The Neumann solution A Neumann similarity solution is available for a problem on a semi-infinite domain, with an initial condition of all liquid water or all ice, thus with a freezing or melting interface moving from one throughout the domain, starting at the fixed boundary. Here, we consider the problem (2.28), with D = oo and T = T . That is, we have a freezing problem, with the freezing interface moving through the domain from z = 0. We follow similar analyses for melting problems by Alexiades and Solomon [3], and Zwillinger [75]. Guided by the heat equations either side of the interface, we seek solutions for t > 0 of the form init  top  T(z,t) = f(ri)  in ice,  T(z,t) = g{n) in water,  where V = A=vt  ( -29) 2  For such solutions to exist, we require (2.30)  s{t) = f3>/i,  for some constant (3 to be determined as part of the solution. The Stefan problem (2.28) then becomes 0 /(0) = T  boU  Ciwaterg". + | W  9iP) = o, \pLl3  0<n<p,  f(P) = o, 0  f3<T) <00,  >  g(oo) = T , top  K f'(f3)-K g'(f3). ice  water  )  (2.31)  Chapter  2. Solutions to prototype interface problems  33  The solution of this system is given by (  f( )  = T  V  1-  bot  and g(y) = T  top  e r f(5( 5 ^ = ) \ erf J^A  0<v<P,  erfGr7==)  (2.32)  (l.-^HH)  /3<»7<oo.  er ' /  where the functions erf and erfc are the error function and complementary error function respectively. The value of /?, which gives the interface location, satisfies the equation V^r  2  :  P  _ K  T  water  e  top  Jb-^  i  a  ^ K T  ^ r  ice  bot  erfc(2^^)  r  e  4  °»ce  evi(^=Y  In order t o construct the Neumann solution given above, we must first solve (2.34) for (3. This may be achieved using a Newton iteration, or bisection method, for example. A typical solution is shown in Figure 2.10.  Travelling wave solutions When seeking an analytical solution to a two-phase Stefan problem, it is sometimes more straightforward to construct a travelling wave solution rather than a Neumann solution. A travelling wave solution is one whose profile remains unchanged in a reference frame moving with constant speed. For the following problem, T T(s-)  t  Tt  = a T = 0, ice  =  Oi t T wa er  T{a+)  = 0,  Ls(t)  =  P  z<s(t),  zz  z > s(t),  zz  (2.35)  water  -[KT^  ice  we seek solutions of the form T =  for f < 0 (ice),  T = F ( £ ) for £ > 0 (water), 2  £ = z - ct. • (2.36)  (2.33)  Chapter  2. Solutions to prototype  Neumann freezing problem, p=1000 L=3.3e+005 c  m t >  interface  =42000^=2000  problems  34  =0.55 K =2.2 T =5 T =-10 .... p=0.000342 Be  tot  Figure 2.10: A typical Neumann solution to the freezing problem (2.32)(2.34). Here, c is a constant, which is the speed of the interface. So s(t) = ct, and for c > 0, we have an interface moving to the right. We have deliberately not specified boundary or far-field conditions. Now, with parabolic problems either side of an interface on which we give three conditions, we expect to generate a one-parameter family of solutions. The system becomes -cF{  = =  a F{' 0,  £<0,  -cF ' F (0)  = =  a F% 0,  £>0,  PL  =  F (0) }  2  2  ice  water  - %±?A , a  >  (2.37)  J  which has solution Fj(0  =  A (l-e  F (£)  =  A {1 - e ~ ^ ^ )  2  x  - uJ)  £<0,  e  4  2  £ > 0,  (2.38)  Chapter  2. Solutions to prototype interface problems  35  where .  ®ice  I  iCire \  T  i  water A  (2.39)  ®water  We indeed have a one-parameter family of solutions; that is, a family of solutions parameterized by A . Such travelling wave solutions can be constructed on finite domains by specifying the appropriate boundary conditions. O n an infinite or semi-infinite domain, notice that F ( £ ) —> A as £ —* oo, and so A is seen to be the far-field temperature. However, F\{£) is exponentially large as £ --> — o o , and this we note again that this analytical solution serves mainly as a mathematical tool rather than having any physical meaning. Furzeland [31], for example, uses exact solutions such as those presented here for evaluating the performance of numerical methods applied to Stefan problems. A typical profile is shown in Figure 2.11. 2  2  2  2  Travelling wave profile with c=l. A^=3  Figure 2.11: A typical travelling wave solution to the freezing problem (2.37).  Chapter  2.2.3  2. Solutions to prototype interface  problems  36  Formulation of the enthalpy method  In most cases requiring the solution of a two-phase Stefan problem, a numerical solution is sought. Let us consider strategies for the numerical solution of the problem 2.1. One option, as discussed in Chapter 1, is to employ a front tracking method. This is a method which explicitly moves the interface at each time step. We consider a finite difference method which updates the temperature, and moves the interface between two time levels, which we denote n and n + l , and outline this method in Algorithm 2.1.  Algorithm 2.1  (Front tracking) Given an interface location s and temperature at a time level n. n  1. Discretize the interval 0 < z < s , and solve the heat equation in the ice, subject to T = 0 at the interface. n  2. Discretize the interval s < z < D, and solve the heat equation in the water, subject to T = 0 at the interface. n  3. Use the Stefan condition (2.26) to calculate the velocity s. 4- Move the interface according to s  n+1  = s + ks. n  5. Go back to step 1. Here, k is the size of the time step between time levels n and n+l. Notice that with each time step, remeshing takes place in steps 1 and 2. A uniform grid on (0, D) will not suffice for front-tracking, since the interface will, in general, not move by exactly one grid point with each time step. A n alternative is to reselect the size of the time step each time the interface is moved, i n order that it will move by exactly one grid point. Details of such approaches to front-tracking are described by Crank [20]. Clearly, in terms of ease of implementation, a front-capturing method is preferable. Further advantages arise in higher dimensional problems, where tracking methods would encounter additional complication due to the interface geometry. Now we consider how we may reformulate the model problem (2.25) over the fixed domain (0,D). Since the heat equation in each of the disjoint ice and water regions describes energy conservation, we aim to write an equation for energy conservation throughout the entire domain. A measure of energy  Chapter  2. Solutions to prototype  interface  problems  37  over the fixed domain must take into account the specific latent heat of phase change. We define the e n t h a l p y (per unit volume), E, by  E(T)  for T < 0 (ice), e[0,pL] for T = 0 ("mushy"), pCwaterT + pL for T > 0 (water).  =  (2.40)  A t a temperature T = 0, the ice or water may be undergoing phase change, so we refer to it as "mushy". Notice the jump in the enthalpy function at T = 0, as shown in Figure 2.12. To deal with the discontinuous conductivity Enthalpy versus temperature  -100  - 6 0 - 6 0 - 4 0  -20  0 2 temperature T  0  4  0  8  0  8  0  100  Figure 2.12: Enthalpy as a function of temperature, function, we once again make the familiar transformation .  _v = f r Y ( ( ) i  Jo  (2.41)  which, in this context, is known as the Kirchoff transformation; the variable v is often referred to as the Kirchoff temperature. Then the interface problem has now been reduced to the following fixed-domain problem for the enthalpy E: E = v 0 < z < D, t> 0 . E(0,t) = pCiceTbot, E(D,t) = pCwaterTtop + PL, . E(z,0) = E (z). • . t  zz  inil  Chapter  2. Solutions to prototype interface problems  38  The challenge now is to devise an algorithm for solving (2.42) numerically, which carefully implements mappings between E, v and T. After giving an initial temperature, we will step in time to solve (2.42), and then the location of the interface will be recovered from the solution E (or v). Let us consider a finite difference method. Suppose we have a grid Zj = (j - l)h for j = 1 , N + 1, where h = D/N. Let VJ be the approximate value of v((j — l)h,nk), where k is the time step to be used in the finite difference scheme. Then an algorithm for the solution of (2.42) is given below in Algorithm 2.2. Firstly, let us write explicitly the relationship between E and v: 1  v < 0,  Kice  E(v)  Sister  y  +  { e[o, L] P  v>0,  (2.43)  v = 0; E < 0,  PCice  v(E)  pjj  ^i^(E-pL)  E>pL,  (2.44)  0 < E < pL .  0  Algorithm 2.2 (Enthalpy method) Given V : n  1. Compute E  using (2.43).  n  2. Update Ej for j = 2 , N  3. Recover V  n+1 3  by usingfinitedifferences for (2.42). Obtain  using (2.44).  The algorithm shown solves (2.42). In order to solve the original problem for temperature and interface position, we can simply recover T from the solution V , and then interpolate to find the interface location s(t). It remains to decide on the scheme to use for Stage 2 of the algorithm. A n explicit scheme is easy to implement as follows (2.45)  Chapter  2. Solutions to prototype  interface  problems  39  for j = 2 , N , where fi = j^. However, we have the time step restriction k<h  - mml-^-,-^-—). \ "ice water /  2  (2.46)  In situations where this is too restrictive , we may prefer an implicit scheme. However, an implicit scheme will require some clever way of deciding which enthalpy range each grid point should be i n at the next time step. We now describe two ways of achieving this. Nonlinear Gauss-Seidel iteration Consider the P D E E =  v,  t  zz  where  V< 0 v = 0  V/Kice  e(0,pL)  E(v)={  V/K ter  + PL  wa  E <0  K E ice  =  {  (2.47)  V>0  and so v(E)  ,  0  E  K (E  - pL)  water  [0, pL]  e  (2.48)  .  E > pL  Suppose we seek an iterative solve for the implicit scheme Ej  +1  -  ii (v^  -  2v]  + v£})  +1  (2.49)  = Ef:  A Gauss-Seidel solver will solve the j equation for the j unknown using the latest available values of all other unknowns. Let us rewrite ( 2 . 4 9 ) as th  E-  +1  J  + 2»v  th  = E? + p(vpl+v?£)..  +1 j  (2.50)  So if we solve for each j in order, we can iterate to find v™ Seidel as follows:  +1  E"  1  + 2 ^ f  1  = E?  + ii «  J  +  ,  using Gauss-  (2.51)  at the p iteration. Now the right hand side contains only known terms. However, we still have the nonlinear dependence E(v) given by ( 2 . 4 7 ) . Now let th  V? =  ^  + M « i  1  + ^+i)-  (2-52)  Chapter  2. Solutions to prototype interface  So ip? is known at the j  equation of the p  th  problems  iteration for vj  th  +  40  . So (2.51) is  just E  +  +l j  (2.53)  2pv =^. +1  j  Now examine the three possibilities for vy : ice, mushy or water. Ice V  f  <o  1  = vf /k  < o,  l  1  ice  and then (2.53) gives ,P+I  v /K  + 2p^  +1  J  =  +l  ice  ^<0  (2.54)  and hence "  2p +  3  (2.55)  l/k  ice  Mushy  and then (2.53) gives 0 < iff =  1  (2.56)  < pL.  Water ^  p + 1  >0  1  = v] /k  + pL,  +1  water  and then (2.53) gives ^hP = j / water V +1 K  + pL + 2pV  +1 j  > pL,  (2.57)  and hence (2.58) 2p + \j Kwater  Thus, the sign of and size of ip? determines v? . That is +1  o  <e[o,pL]  er  ^  > pL  (2.59)  Chapter  2. Solutions to prototype interface problems  41  This Gauss-Seidel idea is widely used, since it avoids the need for computing derivatives. The convergence of such a scheme may be accelerated using successive over-relaxation (see [3, 20], for example). Smoothing methods arid Newton iteration A n alternative to the trial-and-error methods for temperature recovery is to make use of smoothed functions again. That is, we replace the discontinuous functions E and v from (2.40) and (2.41) with E (T) £  = (1 - H (T))pc T e  ice  and v„(T) = (1 - H (T))K T a  ice  + HeiTXpcvotcr + pL), _ + H (T)K (T), a  water  (2.60)  (2.61)  where H is a smoothed Heaviside, such as those in Examples 2.1 and 2.2, and e, a are smoothing radii in temperature. We then solve (2.42), for E, and the temperature recovery may then be achieved using a Newton iteration. This idea is particularly helpful for implicit time-stepping schemes, where Newton methods are commonly employed. In Figure 2.13, we plot a smoothed enthalpy function, taking e — 5, and using the smoothed Heaviside from Example 2.2. Enthalpy versus temperature  Figure 2.13: Enthalpy and smoothed enthalpy as functions of temperature.  Chapter  2. Solutions to prototype interface  problems  42  We leave the investigation of the effect of such smoothing strategies on solutions to the enthalpy formulation of the Stefan problem as an open problem.  2.2.4  Mathematical justification for the enthalpy method  Here we follow Alexiades [3] and Crank [20] to show that the Stefan condition at the interface is satisfied by the weak solution to the enthalpy problem E = v v(0,t)=v , t  zz  0  E\ =Q = T  z € (0,D), t € ( 0 , r ) , v(D,t) = v u  (2.62)  E(T ). init  Here, if T is temperature, then E(T) is the enthalpy, and v(T) is the Kirchoff variable. Now consider the diagram shown in Figure 2.14. t A  0=  0  0  (f> = <t>(z,0) Figure 2.14: The (z,t)-plane for the enthalpy formulation.  The curve V has z = s(t) being single-valued, and T = 0 with T continuous across the interface. First consider the whole region G = (0, D) x ( 0 , r ) ,  Chapter  43  2. Solutions to prototype interface problems  and test functions </>(z,t) € C°°{G), with 0(0,t) = </>(D,t) = (/>(Z,T) = 0. Then we define a weak solution to (2.62) to be a pair of functions E, v which satisfy J J  (j)E - (pv dz dt = 0, t  zz  Now, integrating by parts, or simply noting that <j>Et = (<j>E) - E<f> and t  <fyv„ = {^>v - vcj) ) + v<j) ,  t  z  z  z  (2.63)  zz  we see that .  J  J  E<p + v(j) dzdt = J J t  {<j)E) - {(j>v - v<j) ) dz dt.  zz  t  z  z  (2.64)  z  A p p l y i n g the boundary and initial conditions, we find j LE4>t + v<b dzdt  -  zz  - f° <j){z,0)E{T ) dz init  (2.65) + J ]v <t> {D,t)-v <t) (^t)dt, and work with this as our weak formulation. It remains to show that a solution of (2.65) satisfies the Stefan condition on the moving interface. To do this, we integrate in a similar fashion over the two domains G\, G , shown i n Figure 2.14. The curve F divides G into Gi and G , which we consider to be ice and water respectively. c  1  z  0  z  2  2  In G i , we have f J  f  f f  E<}> + v<t) dzdt= t  zz  JGi  {<f>E) -(4>v -v<l> ) dzdt. ' t  z  z  z  (2.66)  JGi  J  Now, the integration this time requires use of Green's Theorem in the Plane (or, for higher dimensions, the divergence theorem). Recall  L  IL ^- i -  pdx+Qdv=  9  9 dxdy  where dTZ is positively oriented in the ( x , y ) plane. So the right hand side of (2.66) is given by RHS  = =  - §  Ci  <t>E dz + {(j>v - v<f) ) dt, z  z  -tf <t>(z,0)E(T )dz-j:v <t> (0,t)dt o)  inU  Q  z  - f .(<f>E)dz + {<{>v ) dt, T  z  Chapter  2. Solutions to prototype interface  44  problems  since v = 0 on T. Thus we have, for G\, / /  Gi  E<Pt + v<j) dzdt  .=  zz  - /  <f>(z, 0)E(T ) dz  S ( 0 ) 0  init  -fivoM0,t)  dt -f -((j>E) r  dz + (4>v ) dt. (2.67) z  Now, for the water region G . Arguing i n the same way, we get 2  / J  G2  E<j> + v<f>„ dzdt  =  t  - f°  <j)(z, 0)E(T )  Q)  dz  init  + Jlv <t> (D,t) dt + j {(j>E) dz + {(j)v ) dt. (2.68) x  g  r+  z  Now, adding (2.67) and (2.68), and letting r~, T+ - » F, we see that JJ E<lH GlUGa  + v<l> dzdf =  -J  zz  <f>(z,0)E(T t) dz  D 0  ini  + J v f) (D,t)-vo<f> (Q,t)dt  (2.69)  T  1(  2  z  + f [<l>E]+ dz + [<j)v )t dt. T  z  Now, comparing (2.65) and (2.69), we have [((>E]t dz + [</nj ]Z dt = 0, z  and hence, since <j> and its derivatives are continuous across the interface,  L  dz + [v }+ dt)-=0.  (2.70)  z  Finally, since <j> is an arbitrary test function, and since z = s(t) on T, we have that fc-J*£'  dt~  which is precisely the Stefan condition.  [£]+'  •  (271) {  l  )  Chapter  2.2.5  2. Solutions to prototype interface  45  problems  Numerical results and convergence study  Now, we illustrate some interesting and important features of the enthalpy solution to the Stefan problem, which will help us to explain the behaviour of capturing methods for more complex problems. Firstly, we consider the problem (2.28) with initially all liquid at T = 3, and subject to temperatures T = — 2 at z = 0 and T = = 3 a t z = D = l . Thus, a freezing front propagates through the domain. That is, our boundary and initial conditions are Tbot — ~2;  Ttop — 3,  Ti it n  - Ttop,  s  ( 0 ) = 0.  The physical parameters we use are shown i n Table 2.1. Table 2.1: Constants for the Stefan problem Symbol  Interpretation  T y p i c a l value  Units  P  density specific heat capacity of water specific heat capacity of ice thermal conductivity of water thermal conductivity of ice specific latent heat of freezing  1000 4200 2000 0.55 2.2 3.3 x 10  kg/m? JK^kgJK^kg-  Cwater  Kwater  Rice  L  (SI) 1  1  Wm^KJkg-  1  5  1  For the results shown here, we have used Forward Euler time-stepping, with N = 20 grid points. In Figure 2.15, we plot a succession of temperature profiles i n increasing time, together with the steady-state solution to the problem, which has the interface at s = 8/11 « 0.73. Evolution towards the steady-state (the dots) is apparent. In Figure 2.16, we plot the interface location as a function of time. The stepwise behaviour of the location is due to the definition that we choose. For simplicity, we take the interface location to be the grid point to the left of the first positive value of temperature. So each jump is by precisely one grid point. We have also plotted the interface location for the associated Neumann freezing problem. While our problem is on a finite domain, we expect the influence of the boundary condition at z = D to be small for small times, and thus the interface location to behave like that of.the Neumann solution for small times. Indeed, Figure 2.16 shows good agreement between our computed s(t) and the Neumann solution s(t) = /3^/i, for small t.  Chapter  2. Solutions to prototype  interface  problems  46  Temperature profiles with increasing time  "0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1  Figure 2.15: Evolution of numerical results using the enthalpy method. In Figure 2.17, we show an important feature of enthalpy method solutions to Stefan problems, namely the stepwise temperature history at a point. This an non-physical effect introduced by the rapid adjustments and subsequent relaxations of the temperature each time the interface advances by one grid point. This is explained in more detail in [3] (Chapter 4). Finally, in Figure 2.18, we plot the interface location versus time for two travelling wave solutions to the Stefan problem. For ,each of the two wave speeds c = 1,3, we plot solutions obtained using N = 80,160 grid points, with num — 400,1600 time-steps, respectively. Convergence to the correct travelling wave solution is suggested by the diagram, but a numerical convergence study is necessary to demonstrate this convincingly. In Tables 2.2-2.3, we plot the computed errors between exact travelling wave solutions with speeds c = 1,3, and the enthalpy method solution, and denote the errors in temperature and interface location by H-ETHI and \\E \\i respectively. In each case, we use explicit time-stepping, with p = k/h fixed, as we use a second order scheme. The number of grid points is N, and num is the number of time steps. We quantify the errors by calculating numerically the norms of the errors made in the temperature T , and the interface location s. Specifically, we compute the time-averaged quantities s  2  Chapter  2. Solutions to prototype  0.9i  problems  47  Interlace position calculated using enthalpy method 1 1 r n  1  0.5  .0  interface  1.5  1  2.5  2  timet  3.5  3  ,,„>  Figure 2.16: Evolution of the interface using the enthalpy method, together with associated Neumann result. Table 2.2: Errors for the enthalpy method, Forward Euler time-stepping, with c = 1. N  num  20 40 80 160  25 100 400 1600  factor  H3rlli 0.8604 0.4259 0.2057 0.1011  2.02 2.07 2.04  factor  \\E \\l 6.3693E-8 2.7862E-8 1.3431E-8 6.3753E-9 s  2.29 2.07 2.11  defined by 1  ]  l  E  T  h  =  ^  1  N  num N  £  E l (***») j=l  ^(zj,t )\,  T  n=l  ^  \\E \\i = num  n  num  V  s  T  \s(t ) n  ct \. n  f n=l  We notice that, despite a second order spatial discretization, the we do not have second order accuracy. The errors decrease by a factor of about 2 each time the grid spacing is halved, rather than a factor of 4. This is to be  Chapter 2. Solutions to prototype interface problems  48  Temperature history at a point, by the enthalpy method 3  1  1  i  » OJS  0  - 1  0  0.5  1  .1.5  2  2.5  3  3.5  Figure 2.17: Temperature history at the point z = 0.5. Table 2.3: Errors for the enthalpy method, Forward Euler time-stepping, with c = 3. N  num  20 40 80 160  25 100 400 1600  II-fir Hi 1.4917 0.7783 0.4283 0.2264  factor  ll^lli  factor  1.92 1.82 1.89  4.2443E-8 2.0526E-8 1.0002E-8 4.8610E-9  2.06 2.05 2.06  expected; due to the stepwise temperature histories, and the fact that we have not dealt with the interface explicitly, lower order errors are introduced near the interface. Further details are available in [3]. Methods are available for improving the accuracy of the computed interface location [65], based upon the known phase-change temperature, and extrapolation methods. Second order accuracy in the interface location may be achieved, but not in the temperature. In our final convergence study in Table 2.4, we show the errors calculated using the enthalpy method for the Neumann freezing problem, where we compute on z 6 (0,1), giving the exact Neumann solution as the boundary condition at z = 1. The factor by which the temperature error decreases is now around 2.6.  Chapter  2. Solutions to prototype iQ-  7  x  6,  0  ,  interface  problems  49  Interface location versus time for exact and numerical travelling waves , 1 1 1 , 1  0.2  0.4  0.6  0.8  timet  1  1.2  1.4  1.6 x 1 0  1.8 -'  Figure 2.18: Computed interface location for travelling wave conditions, with c = 1 and c = 3. Table 2.4: Errors for the enthalpy method, Forward Euler time-stepping, for a Neumann problem.  2.3  N  num  II ET  20 40 80 160  200 800 3200 12800  0.0595 0.0229 0.0087 0.0034  ||I  factor  ||-E's||i  factor  2.60 2.63 2.56  0.0145 0.0068 0.0033 0.0016  2.13 2.06 2.06  The Porous Medium Equation  In this section, we discuss some results and features of the Porous Medium Equation. This is a nonlinear diffusion equation which arises in fluid flow and other applications, which may result in the appearance of a moving interface.  2.3.1  Examples and applications  Consider the flow of an isothermal, ideal gas through a porous medium. The flow obeys Darcy's Law (see, for example [8]), which says that the volumetric flow rate of the gas through the medium is proportional to the pressure  Chapter 2. Solutions to prototype interface problems  50  gradient in the gas. Specifically, the volumetric flow rate, known as the D a r c y v e l o c i t y , is given by U  = _-Vp, A*  (2.72)  where p is the gas pressure, K is the permeability of the porous medium, and fi is the viscosity of the gas. For illustration in this section, we shall only be concerned with one-dimensional problems, so that u = --p ,  (2.73)  z  where z is the space variable. Now, conservation of mass requires that {<t>p) + (pu) t  z  = 0,  (2.74)  where <fi is the porosity of the medium, which we will assume to be constant, and p is the density of the gas. Now, the ideal gas law relates pressure and density such that P = ^ °  '  T  (2-75)  where R and M are the universal gas constant and the molar mass of the gas respectively, and T is the temperature. Thus, for isothermal gas flow through a porous medium, we find that Pt = c(pp ) , z  z  (2.76)  for a constant c. This is a nonlinear diffusion equation with a variable diffusion coefficient p. Notice that it will be of parabolic type, provided that p > 0. Next, we consider problems where the diffusion coefficient may vanish. Consider the motion of a thin liquid droplet on a solid substrate, shown in Figure 2.19. In the absence of any air flow driving the motion, and assuming that surface tension effects are negligible, then the droplet will spread under the effect of gravity alone. The free surface of the droplet is denoted y = h(x,t), and this surface meets the substrate at x = L(t) on the right hand side. The thin film approximation of the Navier-Stokes equations gives an equation for the height h of the droplet (see [53]):  ht - {h h ) = 0. 3  x  x  (2.77)  Chapter  2. Solutions to prototype interface  problems  51  y y = h(x,t)  x - L(t)  Figure 2.19: A thin droplet spreading under gravity over a solid substrate. Here, finding the position of the wetting front L(t) will be part of the problem. We will see similar equations for conservation of liquid mass flowing through porous media i n later chapters. Notice here that the equation ceases to be of parabolic type at the point x = L(t), where h = 0. This type of degenerate diffusion problem can lead to solutions with singular gradients. B o t h (2.76) and (2.77) are examples of the P o r o u s M e d i u m E q u a t i o n , the general form for which is Ut =  (u u ) Xl n  x  (2.78)  where n > 0. In this thesis, we are particularly interested in problems with n = 3. The spreading droplet example above suggests that equations of this type admit solutions with compact support. Next, we will describe some well known analytical solutions.  2.3.2  Analytical solutions  Here, we will mostly concentrate on the Porous Medium Equation (2.78) with n = 3. First of all, let us consider the steady-state problem. That is, 0. Writing this as 0, we readily obtain a general solution u(x) = A(L  (2.79)  Chapter  2. Solutions to prototype interface  52  problems  for constants A and L. Further, the solution A(L-x  u(x) = |  f'  A  0  x<L, x > L,  (2.80)  satisfies (2.79) in a classical sense everywhere except the free interface x = L. A sketch is shown in Figure 2.20. The derivative u is undefined at x = L, and we regard this solution as a weak, or generalized solution. Weak solutions satisfy an integral form of the equation, in the same way as weak solutions to the Stefan problem. We will define a weak solution of the timedependent problem shortly. For now, suppose we want to satisfy a boundary condition u(0) = UQ, and that we must have a total "mass" W in the system. That is, x  I L  r  u(x) dx = W.  Then the constants A and L Jo are easily found: L =  1/4  5W 4  M  0  Figure 2.20: A solution of the steady-state Porous Medium Equation. Now let us consider the time-dependent problem, u = {u u ) , 3  t  x  x  (2.81)  and the sketch of a solution profile in Figure 2.21. Again, we seek solutions with compact support. A n y such solutions with discontinuous derivatives at the moving interface must be regarded as weak solutions. Such weak  Chapter  2. Solutions to prototype  interface  53  problems  u(x,t)  x = L(t) Figure 2.21: A solution of the time-dependent Porous Medium Equation. solutions are discussed in detail by Elliott and Ockendon [23], and we will not reproduce their work here. We simply state, following the summary given in [43], that a weak solution of the problem (2.81) is a continuous, non-negative function u(x, t) for all £ > 0 and for all x, such that r  T  /  r°° u / -  f°° dx dt + / ip(x, 0)u(x, 0) dx = 0,  z  + wp  t  (2.82)  for all test functions ip which vanish at infinity and £ = T , and which have continuous first derivatives. We shall not discuss the theory of weak solutions further. Rather, we now seek to construct solutions. Assuming that we can find a solution with compact support, such as that shown i n Figure 2.21, what structure will it have near the moving interface x = L(£)? One way to answer this is to suppose that the solution to equation (2.81) takes the form u~f{t)(x-L) ,  r>0,  r  as i - > L ( t ) " .  (2.83)  Substituting, we have f'.{x - L) + f.r{x - LY' r  1  L(t)  ~ / . r ( 4 r - l)(x - L ) ~ . 4  4 r  2  Then, for a nonzero, finite speed L(t), we require the balance r — 1 = 4r — 2,  r =]-.  Alternatively, a more formal approach considers conservation of "mass" at the interface (see [54], for example), which gives L(t) = -  l i m (u%.  x—>L(t)~  (2.84)  Chapter 2. Solutions to prototype interface problems 54 Substituting the form (2.83) again gives r = | for a nonzero, finite speed. So the solution needs to have this certain singularity structure i n order for the interface to move.  Barenblatt-Pattle "Spreading blob" solutions The Barenblatt-Pattle solution of the Porous Medium Equation is a similarity solution of the form  u(x,t)=t F( ), a  V  V  = ^.  '  A . number of authors have investigated the behaviour of such solutions; here, we follow the introduction given by Lacey et al [43]. The values of a and fi are found by the conditions that the solution be self-similar, and that the total mass contained within a "blob" is conserved. For the general equation (2.78), we find . ' ' ^  For the case n = 3, the solution is „(x, )=( t  ^)" '3  , / 5  (« -*) 2  , / S  -  I 0,  W<«  ( 1 / 5  .  (2.85)  \x\ > at ? , 1  5  where a is a constant depending on the total mass i n the system. In Figure 2.22, we show some solutions with a = 1, for various times, and the spreading blob shape is clear.  Travelling wave solutions As for the Stefan problem, we can also construct a travelling wave solution. For a travelling wave speed c, we seek solutions to (2.81) of the form u(x,i) = / ( 0 ,  i =  x-ct.  The problem reduces to the O D E -of  =  (fry,  which has a solution / ( 0 = ("3c£)  1 / 3  .  Chapter  2. Solutions to prototype interface  problems  55  Barenttatt-Paltle "spreading blob" solutions oJ the PME. with a=l  Figure 2.22: Barenblatt-Pattle solutions of the Porous Medium Equation. Thus, our travelling wave solution is {-Zc(x - ct)f' , x < ct, z  «.<*•*)-{o7  X >  ct.  (2.86)  Now, as before, we continue by describing numerical methods for the solution of Porous Medium Equation interface problems, and will return to our analytical solutions to evaluate the performance of these methods.  2.3.3  Numerical methods  The moving interface problems arising from the Porous Medium Equation require special care when developing numerical schemes for their solution. The singular gradients at the interface may introduce problems, while conservative schemes are key to finding a solution with the correct interface velocity. Consider the equation (2.81). The gradient u becomes singular at the interface where u = 0. To resolve the sharp interface, adaptive gridding may be employed to add resolution near the interface. However, we aim to develop a capturing method which, like the enthalpy method for the Stefan problem, does not require any explicit handling of the interface. One approach, which has appeared for similar problems in the fuel cell literature (see [50, 51], for example), is a computational regularisation. This amounts to replacing equation (2.81) with u = {{v? + rf)u ) , (2.87) x  t  x  x  Chapter  2. Solutions to prototype interface  56  problems  where rj > 0 is a computational parameter. The effect of this is to remove the degeneracy and singularity from the problem at u = 0, as the diffusion coefficient remains positive. Thus, the problem remains strictly parabolic, and the sharp interface is smeared out. A s an alternative, we consider the method described by Evje and Karlsen [25]. Consider the problem u '=(u u ) , 0<x<l,t>0, u(x,0) = u (x), ~u u \ =Q , 3  t  x  x  .  0  ,  2  g  g  ,  s  x {0tt)  L  - u u | ( M ) = QR, 3  x  where Q L , R are fluxes at the left and right hand end points. Now consider a finite difference discretization of the problem, with time-step A;, h = 1/N and Uj is the approximation for u((j — l)h, (n — l)k), for j = 1, ..,N + 1. In [25], it is stated that a naive discretization of the kind  (2.89)  j = 2,..,N, will not necessarily result in a conservative scheme. presented in [25] is to rewrite the problem as u = \{u% , 0 < x < 1, u(x,0) = u (x), t  T h e alternative idea  t>0,  x  ,  0  2  g  o  .  —^C" )*l(o,*) = Q L , 4  ~i( )x\{i,t) u4  = QR,  and then discretize, using standard centered differencing for {u ) . We now show how a conservative scheme results, using a ghost point method for the boundary fluxes. For illustration, we use Forward Euler time-stepping, but the same argument will follow for implicit schemes. Step i n time using the following scheme: 4  xx  t^+i = V» + t {{U^Y - 2(C/;) + (U^f) 4  where fi = k/h . point method: 2  ,  j = 2,.., N,  (2.91)  Now derive discrete boundary conditions, using a ghost  57  Chapter 2. Solutions to prototype interface problems  U?  =  +1  U» + % ((t/ ") - (£/») ) + 2 | g 4  4  2  L  (2.92)  ffti = ^ i +f ( ( ^ ) - ( ^ i ) ) - 2 ! Q « . 4  4  +  +  The scheme must conserve "mass" at each time step. That is, we must satisfy / u(x,nk) dx Jo  =  / u(x, (n + l)fc) dx, Jo  in a discrete sense. Suppose we require that our scheme conserves mass using trapezoidal approximation for the integral. Then we have mass  n+1  =  IK+2^  +  1  m) -(u?)*)} 4  {v? + i {(u?-i)  A  + 7 ^ {(C/?) - ( l ? ) + 4  That is,  + r a  mass  n+1  4  - ( ^ + l )  4  \v? + (j2  ij  h  U  +  ~  2  ( ") + (^JVi) )} t7  4  4  E f = ((^"-i) 2  4  4  - 2(£/;) + ( t ^ ) ) 4  4  } -  +  - Q*j+S,  (2-93)  where 5  =  hi { (C/ ") - ([/?),4 4  2  + E f = ((f^i) 4  2  m f + ( ^ i )  4  ) + ( ^ )  4  - (^+i) } • 4  (2.94)  Chapter  2. Solutions to prototype interface  problems  58.  We recognize the sum in (2.94) as a telescoping sum, and find that 5 = 0. Hence, equation (2.93) gives mass  =. mass  n+1  + k(Q  n  L  - QR),  (2.95)  which says that mass is conserved. Profiles obtained using conservative and nonconservative schemes for u,=(uu ) 3  0  :—1  1  1  0 .1  0.2  0.3  1  1  0.2  0.3 l  0.4  0.3  0.4  T—:  1  0.4  '  1  1  1  I  0.5  0.6  0.7 i  0.8  0.9  1  i  0.5  0.6  0.7  0.8  0.9  1  1  l  <  I  0.7  0.8  0.9  l  l  1  0  0.1  0  0.1  0.2  l  '*>ooot 0.5 0.6 . 1 i  -  = 0.5 0  0  1  0.1  0.2  0 .1  0.2  0.4  0.3  0.3  0.5  0.6  1  1  i  0.7 l  0.4  0.5  0.6  0.7  0.8  0.9  1  > conservative x nonconservative — exact 0.8  0.9  1  X  Figure 2.23: Evolution of profiles given by conservative and non-conservative schemes applied to u = (u u ) . 3  t  x  x  In Figure 2.23, we show a succession of profiles obtained numerically using a conservative scheme, and a non-conservative scheme, applied to the same problem. While both give solutions with compact support which "look reasonable" , the solution obtained using the non-conservative scheme is clearly wrong. For problems with degenerate diffusion-type terms in later chapters, we shall use Eyje's spatial discretiaztion to help ensure conservative schemes.  2.3.4  Numerical results and convergence study  In Tables 2.5-2.6, we show numerical convergence studies for solutions to (2.81) on —0.5 < x < 0.5, with travelling wave solutions from (2.86) as initial and  Chapter  2. Solutions to prototype interface problems  59  boundary conditions. We use Evje's discretization (2.91). Specifically, we compute the time-averaged errors in computed values of u, and the interface location L, defined by ^ ll-^ulll  ^ num N jTT /  =  num N  ^\ ( j,t ) ' *r-i u  . /  z  n  — U t(Zj,t )\ exac  n  ,  n=l ]=1 1  num  \\E \\I  =—y;iL(i )-ci |.  L  n  num  n  f n=l  As with the solutions of the Stefan problem, we have used a second order capturing scheme to approximate solutions of a moving interface problem, in which we do not treat the interface explicitly at all. Again, we do not achieve second order accuracy, but a factor of around 2.4-2.6 decrease for each halving of the grid space is clear. Table 2.5: Errors for Evje's method, Forward Euler time-stepping, for the Porous Medium Equation on —0.5 < x < 0.5, with p = 0.1 and c = 3 N num factor factor II^LIII IlKlli 2.122E-2 20 50 1.348E-2 1.165E-2 1.82 5.217E-3 2.58 40 200 5.689E-3 2.05 2.034E-3 2.56 80 800 1.84 2.45 3.087E-3 160 3200 8.297E-4  Table 2.6: Errors for Evje's method, Forward Euler time-stepping, for the Porous Medium Equation on —0.5 < x < 0.5, with p = 0.5 and c = 0.5 N  num  20 40 80 160  100 400 1600 6400  ll^lli 7.058E-3 2.855E-3 1.140E-3 4.652E-4  factor 2.47 2.50 2.45  ll^illi 2.071E-2 1.151E-2 5.723-3 3.040E-3  factor 1.80 2.01 1.88  60  Chapter 3 Nonisothermal, steady-state phase change in porous media In this chapter, we describe the steady-state model problem of phase change in a sand pack, as presented by Udell [64]. His experimental results are obtained by partially saturating a sand-filled tube with water, and then heating from the top, while cooling from the bottom, and measuring temperatures along the height of the pack. At steady-state, the system supports fluid in  Z=D  4-  two-phase  liquid  Figure 3.1: Udell's experiment. the pore space in one of three configurations: 1. Only vapour throughout the entire pack. 2. A vapour-only zone above a two-phase zone in which both liquid and vapour exist. 3. A vapour-only zone above a two-phase zone, with a liquid-only zone underneath.  Chapter  3. Nonisothermal,  steady-state phase change in porous media 61  A three-zone system is shown in Figure 3.1. Udell [64] presents experimental results for the three-zone system, where the two-phase zone is identified as the almost isothermal region between two single-phase regions with linear temperature profiles. He then presents an analysis of the two-phase zone only, under the assumption that it is isothermal. We can use Udell's steady-state two-phase zone model as a starting point, which can then be extended to the free interface problem for a two or three zone system. A disjoint-domain computational method for locating the interfaces in a onedimensional, steady-state, three-zone (vapour/two-phase/liquid) system is described by Bridge et al [15]. It is worth noting that their model relaxes the popular assumptions of isothermal two-phase zones and phase change only at the boundaries, and we wish to keep temperature and condensation rate effects in this work. Motivated by the fuel cell setting, where regions of liquid water only do not occur [12, 21, 56, 58, 70], we will concentrate here on a two-zone system which consists of a two-phase zone and a vapour-only zone. We present a computational method, based on that in [15], to compute the interface location in such a steady two-zone system. The numerical results obtained from this method will be used as benchmarks with which we will compare the results from the unsteady computations to be described in the next chapter. Also in this chapter, we consider the method of Residual Velocities, proposed by Donaldson [22] for the steady-state Stefan problem and described briefly here in Chapter 2, and show how this method may be used to compute steady-state solutions to our model phase-change problem.  3.1  Mathematical formulation of the model problem  In Figure 3.2, we show a schematic of the basic setup, as used by Udell [64]. The porous layer is initially saturated with a certain amount of liquid water, then heated from above and cooled from below. A steady-state is realised, with two distinct zones appearing. In the two-phase zone, z G (0,1/), liquid and vapour coexist in the pore space. In this region, the liquid is driven upwards by capillary pressure, while the vapour is driven downwards by a vapour pressure gradient. In the vapour-only region, z G (L,D), the water vapour is stationary. A primary goal of the works by Udell [64] and Bridge  Chapter 3. Nonisothermal, steady-state phase change in porous media 62  z=D vapour only  3 two-phase  liquid  t  cool  Figure 3.2: A two-zone system for Udell's experiment. et al [15] is to locate the interface z = L. We note that the model presented in [15] assumes a constant vapour density throughout the two-phase and vapour-only regions. Here, we allow for compressible vapour, in accordance with the Ideal Gas Law. The method we use is then adapted from that used in [15]. A list of physical constants and parameters used in our model appears in Table 3.1. In order to locate the interface, we consider the saturation s through the porous layer. The saturation is defined by s  volume of pore space occupied by liquid water  (3.1)  volume of pore space  Therefore, a liquid-only region (which we do not consider here) has s = 1, a vapour-only region has s = 0, and a two-phase region of vapour and liquid has 0 < s < 1. In this model, we assume that s is continuous throughout the porous layer, and that at the interface z = L, we have s —• 0 as z —* L~. Now, we formulate a model for the saturation s (the liquid volume fraction), and the temperature T , and their variations in z, the height up the porous layer. We need to consider energy conservation and mass conservation in each of the two zones. First, we consider the two-phase region 0 < z <,L. +  Chapter 3. Nonisothermal, steady-state phase change in porous media 63  Table 3.1: Constants for Udell problem  Symbol  Interpretation  T y p i c a l value  U n i t s (SI)  0  porosity permeability liquid density specific heat of vapor specific heat of liquid water viscosity of water vapour viscosity of liquid water  0.38 6.4 x 1 0 10 10 4.2 x 10 2.2 x 10~ 2.5 x 10~  m kg/m J/kgK J/kgK kg/ms kg/ms  thermal conductivity of vapor saturated medium  1.0  W/mK  mass averaged density heat capacity product  10  J/Km  latent heat (water liquid-vapor) capillary pressure scaling universal gas constant molar mass of water characteristic vapour pressure characteristic temperature heat flux  2.5 x 10 1.7642 x 10 8.31 18 x 1 0 0.19743 0.03525 ~10  K  Pl Cy  Pv  tie k  v  pc  5 R M a b Q  - . 2  12  3  3  3  3  5  4  3  5  6  3  3  4  J/kg Pa J/mole K kg/mole Pa K' W/m 1  2  Conservation of liquid mass gives (AW) = r. z  (3.2)  Here, p is density, u is the superficial Darcy velocity, and the subscript I denotes liquid. The source term T is the rate of production of liquid water, so we shall refer to this as the condensation rate. In the same way, conservation of vapour mass i n the two-phase zone reads (p u ) v  v  z  = -r,  (3-3)  Chapter  3. Npnisothermal,  steady-state phase change in porous media 64  where subscript v now denotes vapour. The term on the right hand side is the rate at which vapour mass is produced. That is, — T is the evaporation rate. Now, conservation of energy, neglecting convective effects as in [15], is given by 0=  (kT^  + KapY,  (3.4)  where K is the effective thermal conductivity of the liquid-vapour saturated porous medium, and h is the specific heat of vaporization of water. A s in [15], we will neglect saturation effects on the thermal conductivity, and assume a value for this mass-averaged quantity. The Darcy velocities ui, u are the average volumetric flow rates of liquid and vapour, respectively, through the pore space. Darcy's law gives relations between the velocities and the pressure gradients in the two fluid phases as follows: vap  v  Ul  P>i  u  v  {(Pi)z + Pi9),  (3.5)  (p + Pv9).  (3.6)  z  AV  Here, K is the permeability of the porous medium, pi is the liquid pressure, p is the vapour pressure, and g is the acceleration due to gravity. The quantities K i and K are the relative permeabilities of liquid and vapour respectively. These relative permeabilities account for the decrease in mobility of one phase due to the presence of another, and hence depend on the saturation s. In particular, we require that K i is unity when only liquid is present, and zero in pore space occupied by vapour only. That is, K I should be an increasing function of s. B y similar reasoning, we require that n is a decreasing function of s. Here, we will use the cubic forms suggested by Udell [64], following the empirical results of Fatt and Klickoff [28]: r  rv  r  T  rv  Kri Kl r  =  ^  (3.7)  ,  = (1 ~  S f  (3-8)  We now seek relationships between between the pressures and the primary variables s and T . A major assumption of the model presented i n [15] is that the vapour i n the two-phase region is fully saturated. In keeping with this model, we will assume here that the vapour is fully saturated, and that the  Chapter  3. Nonisothermal,  steady-state phase change in porous media 65  temperature and saturation pressure, p t approximately obey the exponential relation p (T) = ae , (3.9) sa  bT  3at  where the constants a and b are fitted to data from saturated steam tables in [1], for example. Then, in the two-phase zone, we take p = p (T), as explained by Baggio et al [7]. Now we present a constitutive Jaw for the liquid pressure. A t the pore scale, the interfacial tension between liquid and vapour phases gives rise to capillary effects. The capillary pressure, p , is defined as the difference between the vapour and liquid pressures, sat  c  Pc = P~Pi-  (3-10)  As shown by Leverett [45], 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. Udell [64] correlates this Leverett function to write the capillary pressure as c  p (s) = <5 J(s),  c  (3.11)  c  where J(s) is given by J(s) = 1.417(1 - s) - 2.120(1 - s) + 1.263(1 - s) , 2  3  (3.12)  and  '"iffwhere o is the vapour-liquid interfacial tension, and (fi is the porosity of the porous medium. Also, we should note that the function J(s) given in (3.11) is an empirical relationship which will depend on the particular porous medium being used. The correlation given in [64] 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, in the two-phase zone (and indeed throughout the entire system), we assume that the water vapour obeys the ideal gas law, namely P=?fPvT,  (3.13)  where R and M are the universal gas constant and the molar mass of water respectively. The work presented by Bridge et al [15] assumes a constant  Chapter  3. Nonisothermal,  steady-state phase change in porous media 66  vapour density, but we wish to include compressibility due to thermal effects here. In summary, the two-phase zone conservation equations (3.2)-(3.4), together with all the constitutive relations and empirical laws, form a second order system i n three unknowns s, T and T. Furthermore, we can eliminate the condensation rate T by summing (3.2) and (3.3), arriving at the system (f/{j- {Psa {T)-5J(s)) z  + g)  t  Pl  (3.14) [ K -  +^  —  -  ^  -  (  _)  l  ( (T))  8  Psat  +  r ^ — j r - ) ) = * •  (3.15) Equations (3.14) and (3.15), form a coupled, second order, nonlinear system for the two-phase zone variables s and T, i n the region 0 < z < L. T h e system clearly has a singularity and degeneracy at s = 0. In the case of constant vapour density, Bridge [14] shows that s = 0(L — z ) / as z L~, and as such, the degenerate diffusion type term requires careful treatment when constructing a numerical solution. The vapour-only zone L < z < D is more simple. T h e temperature is harmonic, and the saturation is zero everywhere i n this zone. Conservation of mass reads 1  {pvUy) = z  4  0,  (3.16)  while conservation of energy reads 0= T,  (3.17)  zz  since there is no phase change i n the vapour-only region. Again, we assume that the vapour is ideal, and the system (3.16)-(3.17) can be cast as a system for the temperature T and vapour density p i n the region L < z < D. For this one-dimensional problem, exact solutions are available in terms of the interface location L, the fixed boundary temperature T(D) = - T i , and the temperature and pressure at z = L , which we will denote T and p respectively. T h e temperature i n the vapour region is given by v  +  T(z)=T  +  +  Ti -  +  T  +  +  -^- (z-L),  1  T  or  T(z)=T  +  +  4-(z-L),  (3.18)  Chapter  3. Nonisothermal,  steady-state phase change in porous media 67  where K is the thermal conductivity of the vapour-saturated porous medium, and q is the constant heat flux through the porous layer. Notice that, if the vapour in this zone is stationary, then the vapour pressure in the vapour-only region is the solution of a separable first order equation, and is given by v  p(z)=P [^-j  •  +  (3.19)  The vapour density in this region may be written _ .  .  x  Pv(z) = P [^J +  ft j  Makv\  •  -(3.20)  These exact solutions will be useful when implementing both the disjoint domain method described in [15], and the Residual Velocity Method proposed by Donaldson [22], where we can iteratively solve problems parameterized by heat flux. Now we examine the boundary and interface conditions required to close the model problem. We have elliptic problems in two variables in the two regions, so we shall specify two conditions on the two physical boundaries z = 0 and z = D, and four conditions on the interface z = L. Since the interface location L is also an unknown in the problem, we require a fifth condition at the interface. A t the boundary z = D, we impose a temperature T i , say. Also, since we consider a closed porous pack, there can be no mass flux across z = D. Hence, at z = D, we have T = Ti,  (3.21)  u = 0.  (3.22)  v  Now, at the free interface, we require five conditions. First of all, the saturation s is zero. Temperature and vapour pressure should be continuous. Finally, mass and energy should both be conserved across the interface, so that the heat flux is continuous, and there is no mass flux across the interface.  Chapter  steady-state phase change in porous media 68  3. Nonisothermal,  In summary, at z = L, we have the following five conditions: s = 0,  (3.23)  [T]+ = 0 ,  (3.24) (3.25)  ]p]t = 0,  (piui + p u )  = (p u )  r  v  (kT -h z  v  v  v  = (KT y  uY v  vapPv  (3.26)  ,  +  x  .  (3.27)  We note here that condition (3.27) corresponds to a singular evaporation rate at the interface. A t the boundary z = 0, we have an imposed temperature and zero mass flux, which give (3.28)  T — Jo,  Piui + p u v  v  (3.29)  = 0.  We have specified five conditions at the interface, but we see that a unique solution will not be available. Imposing the boundary condition (3.29) leaves the interface condition (3.26) redundant, and a further condition is required. In order to determine the vapour pressure uniquely throughout the porous layer, we now impose a global integral constraint on the system. Experimentally, the total water mass i n the system can be controlled, and will be known. Suppose that the fixed water mass per cross sectional area of the saturated porous pack is W. Then we have the integral constraint /  Jo  (s-pi + (1 - s)p ) v  dz  +  /  p dz v  =.W,\  (3.30)  J L  which closes the system. It is clear that the three control parameters for the experiment are the temperatures To, T i and the water mass W. In the next section, we will describe a numerical method which, given these three parameters, will calculate the location of the free interface z = L.  3.2  An iterative disjoint-domain method  The system we are trying to solve is summarized in Figure 3.3. In [15], an algorithm is described for solution of the steady-state problem with a liquid zone, a two-phase zone, and a vapour zone, where there are two interfaces  Chapter  3. Nonisothermal,  steady-state phase change in porous media 69  to find. Here, there is just one interface. There is an additional nonlinearity in the model here, as we have allowed for compressible vapour. Now we will describe an adaptation of the existing method to solve our problem for the interface location L. T = T i , piui + p u v  = D  = 0  v  S  mass equation (3.16) vapour, (p, T)  energy equation (3.17)  z = L  interface conditions 3.23), (3.24), (3.25), 3.27) mass equation (3.14) energy equation (3.15)  two-phase, (s, T)  z = 0  J_ T = T , piui + p u 0  v  v  = 0  W i t h integral constraint (3.30). Figure 3.3: System for steady-state solution of the free interface problem. Our solution method is parameterized by a heat flux q and the saturation at z = 0, which we denote s . Now, if the heat flux across the boundary 2 = 0 is q, then we have . 0  KT  Z  - h  u  vapPv  v  = q,  at z = 0.  Notice that this then allows us to integrate the energy equation (3.15) once, leaving KT — h p u = q, for 0 < z < L. ' (3.31) Z  vap  v  v  Also, taking the zero mass flux condition at z — 0 allows us to integrate the  Chapter  3. Nonisothermal,  steady-state phase change in porous media 70  mass equation (3.14) once, leaving Pm + p u v  = 0,  v  for 0 < z < L.  (3.32)  Thus, given q, we reduce the problem in the two-phase zone to a coupled system of ordinary differential equations, f/  (j-(p (T)-5J{s))  +  sat  g)  Pl  (3.33) ir M.h dT  M  P**22  \*(  d  h  0  (n  ^  ,  Mg (T)\ Psat  (3.34) which will be solved as an initial value problem for 5 and T , with initial values So, T . This initial value problem in z is solved numerically until s = 0. This requires some care, since the coefficient s in the first term of (3.33) makes the problem singular at s = 0. One approach is to reformulate and solve for z(s) rather than s(z), which is the method we have used here. Another approach is to make the change of variables w = s , which leaves a regular problem for w. This idea has already been discussed with respect to numerical methods for degenerate diffusion problems, and will feature again in the next chapter. Once we have solved until s = 0, then we stop, and we have found the interface location L. Then it remains to solve the problem for pressure and temperature in the vapour region, subject to the remaining boundary, interface and global conditions. Given continuity of temperature (3.24) and pressure (3.25), we find the values of temperature and pressure just above the interface, namely T and p . Also applying the heat balance (3.27) at the interface, and no mass flux at z = D, we find that the problem (3.16)-(3.17) has the analytical solution (3.18)-(3.19). So, given q and So, it is relatively straightforward to solve the system of equations and find the interface L. However, the condition T = T i at z = 0, and the global integral constraint on the mass W (3.30) have yet to be satisfied. The idea is to iterate on q and So until these two conditions are satisfied. Suppose that, for given q and SQ, we have computed a solution using the algorithm described here, which has T = T at z = D, and which has a total water mass (per unit area) of W. Then, defining 0  3  4  +  +  Chapter  3. Nonisothermal,  steady-state phase change in porous media 71  we seek zeros of the function Q. Following [15], we compute the zeros using a quasi-Newton method which uses centered differences for the numerical derivatives.  3.2.1  Numerical results and discussion  Our Newton iteration is seen to be quite sensitive to initial guesses. That is, we require a good initial guess in order for the iterations to converge. Good initial guesses in many of our computations have been generated by a method of continuation, similar to that described by Zwillinger in [75] (Chapter 168). Suppose we have base values of the control parameters T , T i and W, for which we have a numerical solution with a two-zone system, and "target" values of these parameters, for which we wish to solve a new problem. We take the values of q and SQ corresponding to the base values as initial guesses for a problem with slightly adjusted values of the control parameters. The numerical solution of the resultant problem generates new values of q and SQ, which are then used as initial guesses for the next slightly adjusted problem. This process is repeated as we gradually adjust the control parameters to their target values. In practice, we have adjusted the control parameters one at a time. That is, we have continued in the parameters TQ,T\ and W independently. In Figure 3.4, we plot the numerical solutions obtained using the iterative method, with T = 320, T i = 450 and W = 30. Our computations have used the parameter values given in Table 3.1, together with the height of the sandpack used by Udell [64], D = 0.254m. Also, we have neglected gravity effects, so that g = 0, and we will continue to do so throughout the remainder of this work. Notice that, by (3.19), this gives a constant pressure throughout the vapour-only region. The interface L « 0.12 is clear, as is the singularity in the saturation gradient as we approach the interface z — > L ~ . We make note of the vivid variation in the vapour density p , which is a feature not included in the previous work [15]. Also, the temperature in the two-phase region is not constant. Two-phase regions in porous media are often modelled as being isothermal, but we are particularly interested i n capturing methods for nonisothermal phase change problems. To underline the effect of not making the isothermal assumption, in Figure 3.5, we plot a close-up of the temperature profile, showing a variation of about 7 K over the two-phase region. The structure o f the temperature profile in this region is the same as that described in [15]. 0  0  v  Chapter  3. Nonisothermal,  steady-state phase change in porous media 72  Steady-state, compressible Udell problem  0  0.05  0.1  0.15 z  0.2  0.25  T1=450, T0=320, q=928.4, s0=0.78, L=0.12  0  0.05  0.1  0.15  0.2  0.25  z  Figure 3.4: Solution profiles for free boundary problem using the iterative disjoint-domain method. In computing the solutions and interface location, we have eliminated the condensation rate T. The structure of this function is of interest. Using (3^4), we calculate the condensation rate for To = 320, T\ = 450 and W = 30, and plot T in the two-phase zone in Figure 3.6. The condensation rate is zero in the vapour-only region. In the two-phase region, far from the interface, T is positive, and increasing towards the cold lower boundary, as we would expect. A s we approach the interface L, we know that the condensation rate becomes singular and negative, signifying an evaporation front. The liquid water that is driven upwards towards the interface evaporates very quickly, and the resulting vapour is then driven downwards. A s in [15], the phase change is concentrated in a layer near the cold boundary and at a front at the interface between the two-phase and vapour zones, but is nonzero throughout the two-phase region.  Chapter 3. Nonisothermal, steady-state phase change in porous media 73 Structure of temperature profile  Figure 3.5: Close-up of temperature profile for free boundary problem using the iterative disjoint-domain method. In Table 3.2, we show further numerical results for various values of the control parameters T , T i , and W. In the first three rows, we have T = 340 and W = 20, while we vary the upper boundary temperature Ti. Decreasing Tj corresponds to an decrease in the magnitude of heat flux through the porous layer. This results in an increase in the length of the two-phase region, as described in [14]. 0  0  Table 3.2: The effect of varying boundary temperature. To  Ti  W  340 340 340 375 375 375  600 550 400 670 550 500  20 20 20 36 36 36  Q 1595 1356 542 2520 1858 1541  so  L  0.51 0.41 0.21 0.61 0.42 0.35  0.092 0.100 0.145 0.140 0.160 0.174  In Table 3.3, we show results of another numerical experiment. This time,  Chapter  3. Nonisothermal,  steady-state phase change in porous media 74 Condensation rate In the two-pnase zone  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09  0.1  0.11  Figure 3.6: Condensation rate in the two-phase zone. we fix the boundary temperatures, and vary the total water mass W. A s W increases, the boundary saturation s increases, and so does the length L of the two-phase zone, as expected. 0  Table 3.3: The effect of increasing mass. To  T\  W  320 320 320 320  450 450 450 450  15 20 25 30  Q 790 848 893 928  so 0.29 0.42 0.61 0.78  L 0.093 0.105 0.113 0.122  Finally, we note that the method may fail, if at the stage of computing solutions to the two-phase zone initial value problem (3.33)-(3.34), the computed interface location L is greater than the height of the layer, D. This signifies that a two-phase zone is supported by the heat flux, without an accompanying vapour-only region. Indeed, given a heat flux and initial saturation, (3.33)-(3.34) may be used to easily solve for the saturation and temperature in a single-zone, two-phase system. A n y numerical capturing method we develop should general enough so that it can deal with both single and two-zone systems.  Chapter  3.3  3. Nonisothermal,  steady-state phase change in porous media 75  The method of residual velocities  In this section, we extend the method of residual velocities described in the previous chapter to our steady-state phase change problem in porous media. Consider again the problem shown in Figure 3.3. We follow the same idea as for the Stefan problem, by giving an initial interface position L, and solving the problem subject to all but one of the interface conditions. The computed residual in this interface condition will be used as the interface "velocity", and we step in "time" to evolve the interface to the correct steady-state position. There are four interface conditions to satisfy. For illustration here, we describe the method, using the residual in interface condition (3.24) for the interface velocity. For a given L, we solve the boundary value problems on either side of the interface by using an iterative method which will ensure that the global constraint (3.30) is satisfied. The method we describe below is again based on one-dimensional integration which allows us to solve a series of initial value problems. First, make a guess for the heat flux q and the temperature just below the interface T~. Then s and T in the two-phase region come from the solution of the initial value problem given by the ordinary differential equations (3.33)(3.34), together with the initial conditions T ( L ) = T~, s(L) = 0, and solving on 0 < z < L. Using the continuity of pressure (3.25), we have p = p (T~), and then the solution of the vapour-only problem is given by +  SBtt  T(«)=T  1  + -|;(z-D),  p,  p  +  0 < z < D  (3.36)  so that D). Bearing in mind that we are using (3.24) for the interface velocity, the two conditions left to satisfy are T(0) = T , and the mass constraint (3.30). So we iterate on q and T~ until these are satisfied, using a quasi-Newton method again. Once all the other conditions are satisfied, we define the interface "velocbot  ity" to be L(t) = T~-T , +  and step in time using a numerical integrator.  (3.37)  Chapter  3.3.1  3. Nonisothermal,  steady-state phase change in porous media 76  Numerical results and discussion Evolution to correct steady-state using method of residual velocities  Figure 3.7: Convergence to correct interface location, using method of residual velocities, for T = 320, T i =450, W = 30. 0  In Figure 3.7, we plot the interface location as a function of "time", as given by the numerical solution of equation (3.37), for the problem with T = 320, T i = 450 and W = 30, for three different initial values of the interface position. The interface location is seen to evolve and converge to the correct value L = 0.12, for all three of the initial values. In Figure 3.8, we plot the computed L ( i ) for the problem with T = 375, T i = 500 and W = 36, for which we have previously computed an interface at L = 0.174. 0  0  3.3.2  A model problem in higher dimensions  The extension of this method to higher dimensions will require solutions of boundary value problems on irregular domains, and a consideration of the interface geometry. A s discussed in earlier chapters, front-tracking methods, even such as this one which solves elliptic problems in the disjoint domains, often require significant computational effort in remeshing and updating the interface position. We will not attempt to apply this method to our phase change problem in higher dimensions. Instead, we propose a related model problem which falls into a class of free interface problems which have been studied by Chen and Wetton [19]. They discuss a class of problems which  Chapter  3. Nonisothermal,  0.21.  0 1 2  l 0  1  , 0.5  steady-state phase change in porous media 77  Evolution to correct steady-state using method of residual velocities 1 1 1 i 1 1  ,  1  , 15  , 2 Mme-t  . 2.5  _. 3  1  3.5 x  1 4 , -a 0  Figure 3.8: Convergence to correct interface location, using method of residual velocities, for T = 375, 7\ = 500, W = 36. 0  have a free interface which lies between two regions, in each of which a vector Laplacian type problem must be satisfied. To start with, consider a generalization of our current model problem to two dimensions, as shown in Figure 3.9. The five conditions at the free interface are now 8 = 0,  (3.38)  [T\± = 0,  (3.39)  [p]t = 0, (piui + p u )~ v  v  (3.40)  n = (p u )  +  v  v  .n,  (KWT - hvoppviij) . n = ( x V T ) . n . +  (3.41) (3.42)  where n is the unit normal to the interface. Now we consider a reduction' of this model problem, as we seek a problem of vector Laplacian type in each subdomain i l i , and interface conditions which are linear in the dependent variables. This class of problems is amenable to the linear analysis presented by Donaldson [22] for scalar problems, and extended to vector problems in [19]. In the following reduction of the physical problem given in Figure 3.9 and (3.38)-(3.42), we replace most of the physical constants with unity, and remove the nonlinearity in the governing equations by assumptions and simplifications which leave a related, ) 2  Chapter  3. Nonisothermal,  fl  +  steady-state phase change in porous media 78  (vapour)  V.(p u ) v  =0  v  AT = 0  V.(piui + p u ) v  V.(KVT  -h p u ) vap  v  t  = 0  v  v  = 0  interface conditions (3.38)-(3.42)  f T (two — phase)  Figure 3.9: Steady-state system in higher dimensions. but not physical, model. Also, the reduced model lacks the singularity and degeneracy associated with the physical problem. However, it is with a view to future work in extending the residual velocity method to our nonlinear boundary value problems that we present this loosely related model. In the vapour region f2 we already have Laplace's equation for the temperature. Furthermore, if we assume that the vapour density is almost constant, then we approximate the mass equation by 1(  A p = 0. In the two-phase region f i , we seek a saturation equation. Taking a constant relative permeability, assuming a linear function for the Leverett function J(s), and just keeping the highest order term in s, we replace mass conservation with . A s = 0. 2  Finally, in f^, we replace the energy equation with A T = 0,  Chapter  3. Nonisothermal,  steady-state phase change in porous media 79  which is equivalent to neglecting phase change effects in this two-phase zone. We make similar reductions to the interface conditions (3.38)-(3.42), such that they still loosely represent continuity of saturation, temperature, pressure, mass flux and heat flux. The conditions we propose are: = 0,  (3.43)  [T]t = 0,  (3.44)  p =T~,  (3.45)  8  +  (VT + Vs)  -  .n = ( V p )  +  O J V T - V s ) " .n = ( V T )  .n,  (3.46)  .n.  (3.47)  +  A similar model problem has been studied, in [19], and we leave the extension of this model to include the nonlinear, singular features of the physical problem as future work.  3.4  Benchmark solutions  The steady-state solutions to the one-dimensional Udell problem that we have described in this chapter will be useful in validating computed solutions to the time-dependent problem that we will study i n the next chapter. We shall henceforth refer to solutions generated using the disjoint-domain method here as benchmark solutions. We will check that solutions of timedependent problems with the same values of To, T i and W, with no mass flux across the physical boundaries, evolve to the correct steady-state. The implementation of the interface capturing method that we develop in the next chapter benefits from a relatively straightforward extension from one dimension to two dimensions, and our benchmark solutions will also be used to validate solutions to the two-dimensional problem, in cases where we give two-dimensional initial data, but one-dimensional boundary data.  80  Chapter 4 The M 2 mixture method for the transient Udell problem In this chapter, we describe the time-dependent extension of the steady, onedimensional phase-change problem presented in Chapter 3. We carefully derive the conditions on the moving interface, and show that these will be difficult to implement in disjoint-domain numerical solution methods which involve front tracking. Clearly, front capturing methods appeal for this type of problem. We shall describe a reformulation of the problem over a fixed domain, in which the interface conditions are not explicitly imposed. This formulation is based in part on the mixture formulation presented by Wang and Beckermann [68], who implement a numerical solution method for the case of an isothermal two-phase region. The convergence of such numerical schemes for nonisothermal phase change in porous media has not been well demonstrated in the literature, owing to the lack of exact solutions of such problems. We describe the finite volume solution of the nonisothermal mixture problem, and demonstrate the validity of this method by establishing exact similarity solutions for reduced model problems, and presenting a numerical convergence study. We shall name our new reformulation and solution procedure a the " M 2 " mixture method.  4.1  Mathematical formulation of the model problem  Here, we present a model based on the time-dependent extensions to the steady-state conservation equations (3.2)-(3.4), and (3.16)-(3.17). Firstly, we consider the two-phase zone 0 < z < L(t), where the interface location is now a function of time. In this region, we have conservation of liquid mass. T h e mass (per unit volume) of liquid i n a control volume is 4>pis, where <j> is the porosity of the medium, pi is the liquid density, and 5 is the saturation,  Chapter  4.  The M2 mixture method for the transient Udell problem  81  as before. Then, conservation of liquid mass is given by {4> s\ + Pl  { m) = r, Pl  (4.i)  t  where, as before, u represents the Darcy velocity, and T is the condensation rate. Similarly, conservation of mass for the vapour phase is given by  (<j>p (l-s)) v  + (p u )  t  v  v  (4.2)  = -T.  z  The energy equation is now (pc) T = KT t  (4.3).  + h F,  ZZ  vap  with a mass averaged product of density and heat capacity appearing. We assume here that the dominant density-heat capacity product is that of the porous medium, such that we can neglect variations in this quantity with saturation. Certainly, this will be true for small values of saturation near the interface. A s in the steady-state problem, we can eliminate the condensation rate between the three conservation equations to give conservation of mass as  (j)(pis + p (l-s)) v  + (piui + p u )  t  v  v  z  (4.4)  = 0,  and (pc) T = KT t  ZZ  - h  vap  ( (</>p (l - s)) + (p u ) v  t  v  v  z  ).  (4.5)  Again, we have a coupled system for the two unknowns s and T in the twophase zone 0 < z < L(t). In the vapour-only region L(t) < z < D, there is no condensation, and conservation of mass gives = 0,  ((f>p ) + (p u ) v  t  v  v  z  (4.6)  and conservation of energy is given by the heat equation with no source term: (pc)T = KT . t  zz  (4.7)  We now have parabolic systems in two unknowns on either side of the moving interface z = L(t), and, as such, we require five conditions to be specified on the interface. The first three conditions are that saturation is zero, the  Chapter  4.  The M2 mixture method for the transient Udell problem  82  temperature is continuous, and the vapour pressure is continuous, giving the conditions = 0,  (4.8)  [T]t = 0,  (4.9)  [p]i = 0,  (4.10)  5  which are exactly the same as conditions (3.23)-(3.25). The two remaining conditions again come from conservation of mass and energy across the interface. These conditions require careful consideration of the effect of the nonzero velocity of the interface, and we derive these conditions in the following subsection.  4.1.1  Modified Stefan conditions at the interface z = L(t)  The well known Stefan condition describes conservation of energy across a freezing/melting interface moving through a body of water, say, separating regions of liquid and solid. Following Crank [20], we have presented a derivation of this condition in Chapter 2. The major assumptions which are made are that the water density is the same in either phase, and that the water in each phase is stationary. Clearly, the Stefan condition must be modified for problems in which there is an interface between liquid and vapour phases, where the fluid on either side of the interface may be moving, and in cases where additional heat sources or sinks exist. W i t h this in mind, we now proceed to find conditions for mass and energy conservation across the interface, in terms of the interface velocity L(t). Firstly, we consider conservation of mass across a general moving interface which separates two fluids, possibly different phases, which may have different densities and velocities. The argument is presented for a one-dimensional problem, which is shown in Figure 4.1. The fluid to the left of the interface has density pi, and is moving with velocity Ui, while the fluid to the right of the interface has density pn, and is moving with velocity un. Suppose that the interface moves a distance 5z in short time 5t, and that it moves with velocity L(t). Mass conservation requires that the difference between the mass in the control volume [L(t),L(t) + 5z] over the time interval [t,t + St] is due to the net mass flux into the control volume during that interval. For  Chapter 4. The M2 mixture method for the transient Udell problem 83 Sz, St <—-*>  •  fluid density p\ ' velocity u\  L(t)  fluid density pn velocity uu z = L(t)  Figure 4.1: Mass conservation at the moving interface. the problem shown i n Figure 4.1, this gives (Pi - Pii)Sz = (Piui ~ P\iuu)5t. Taking the limit as 6t —* 0, we arrive at a condition on the velocity of the interface, that is (pi ~ P\\)L(t) = pim - p u . (4.11) u  u  Further, i n the case of flow through a porous medium which has porosity cf), the interface velocity will be given by (4.12)  0(pi - Pn)L(t) = piui - puuu,  where u\ and u\\ are now understood to be the Darcy velocities. Now, for the Udell phase change problem, suppose we take region I to be the two-phase region, and region II to be the vapour-only region. T h e condition (4.12) holds, but we must consider mixture quantities i n region I. That is, we consider quantities associated with the mixture of liquid and vapour i n the two-phase region. Let us define pi — pis + p (l — s), the mixture density, = piui + p u , the mixture mass flux.  (4.13)  v  piU\  v  v  Then equation (4.12) gives 4> ( (pis + p (l v  - s))i - (p )u v  ) L(t) = (p ui + p u )i L  v  v  -  (p u ) . v  v  n  Chapter  4. The M2 mixture method for the transient Udell problem  84  Now, in view of the continuity of temperature (4.9) and of vapour pressure (4.10), we have that (p„)i — (Pv)n = 0, and hence <t> ((pi ~ Pv)s\ L(t) = (pm + p u )i - (p u ) . v  v  v  v  (4.14)  n  To explicitly find the interface velocity, we must be careful. Specifically, the saturation condition at the interface (4.8) is s = 0 at z = L~. The velocity L(t) must be considered as a limit of an indeterminate form. In particular, the velocity is explicitly given by L  (  t  )  =  ] i m  (w  +w)i-(ft«.)n  (j)(pi -  z-»L(t)-  (  4  1  5  )  p )iS v  We note that such an indeterminate form for the velocity of a moving interface, resulting from a mass conservation argument, also arises when considering the porous medium equation, as we have seen in Chapter 2. The common feature here is the degeneracy as s —* 0. Now, for the Udell phase-change problem, we also require an energy balance across the moving interface, which we will also see to be a limit of an indeterminate form. The argument we present here is an extension of Crank's [20] derivation of the Stefan condition in one dimension. Our problem includes an extra term due to motion and phase change in the two-phase region, and we refer to the energy balance as a Modified Stefan Condition. Consider the problem shown i n Figure 4.2. The interface between the twophase region and the vapourronly region in the Udell problem moves a small distance Sz in the small time interval St. The heat which flows into the control volume during the time interval is  while the heat which flows out of the control volume during the time interval is . K^=-\ dz)  St. x  Now, the heat required to evaporate the liquid which moves up to the interface is (h piUi)St. vap  These three quantities are exactly the same as for the steady-state problem. We now consider the extra heat released or required when the interface moves.  Chapter  4. The M2 mixture method for the transient Udell problem 85 5z, 5t  region I (two-phase)  L(t)  z = L(t)  Figure 4.2: Temperature profile and the moving interface. First consider the case 5z,L > 0. The additional mass which appears in the two-phase zone after the interface has moved is (pspiSz. This must be exactly the amount of mass from region I I which has condensed, and the heat released upon condensation is then h 4>spi5z. vap  Given these four terms, we see that the energy balance across the interface is given by K  dz J i  5t +'h (j)spi8z vap  =  h piui5t. vap  Now, we take the limit as 5z, 8t —» 0 to give h (f)spiL(t) vap  = h piui  (4.16)  —  vap  To find L(t), we note that s —> 0 as z —> L(t) , and evaluate the limit +  .nil  Lit)  =  lim <  h PlUl  dz J i  vap  h (j)spi vap  (4.17)  Chapter  4. The M2 mixture method for the transient Udell problem  86  Now, if we repeat this argument for the case 8z,L < 0, we get exactly the same condition. Notice that, if we take (4.17) as the indeterminate modified Stefan condition which determines the interface velocity, then the mass balance (4.15) becomes - K% ^ |  (pi - Pv)i\h (piui)i vap  (piui+p u )i-(p u )n v  v  v  =  v  :  —j—  ——,  :  as z -> L .  tlvapPl  (4.18) We note here that the indeterminate modified Stefan conditions (4.15) and (4.17) can be thought of as limiting cases of the Rankine-Hugoniot conditions for our system. T h e Rankine-Hugoniot condition states that, for a conservation law of the form A + B = 0, t  (4.19)  z  where B is the flux of the quantity A, then for any smooth space-time curve z = L(t), the jumps i n A and B are related to the velocity L(t) by [A]_L(t) =  for z = L(t),  (4.20)  which expresses the conservation of A across the curve (see, for example [3]). For our problem, this may be easier to see after we have reformulated the system as a mixture problem, so the careful derivations of (4.15) and (4.17) are valuable here.  4.1.2  Model summary  In summary, we have a two-phase region 0 < z < L(i), i n which (j)(pis  + p (l-s)) v  + (piui + p u )  t  v  v  z  = 0,  (4.21)  and (pc) T = KT t  ZZ  - h  vap  ( ( # „ ( 1 - s)) + (p u ) ). t  v  v  z  (4.22)'  The vapour-only region L(t) < z < D has ((j>p ) + ( u ) v  t  Pv  v  = 0,  z  (4.23)  and (pc)T = kT . t  zz  (4.24)  Chapter  4.  The M2 mixture method for the transient Udell problem  87  A t the moving interface z = L(t), the five conditions are s  0  (4.25)  [T} -  0,  (4.26)  \P) -  0,  (4.27)  +  +  <t>((pi-  Pv)s) L(t)  (Piui +  <t>spiL{t)  KapplUl —  1  PvUv)l - (PvV>v)lI>  K  (4.28) (4.29)  A t the upper boundary z = D, we impose a temperature T = T i , and have no mass flux, so that u = 0. A t the lower boundary z = 0, we impose a temperature T = T (< T i ) , and have no mass flux, so that piui + p u = 0. Finally, we give initial profiles of saturation, temperature and pressure throughout the entire system at time t = 0. v  0  4.2  v  v  Fixed domain, mixture formulation  A n y numerical method which is based on solutions in the disjoint domains 0 < z < L(t) and L(t) < z < D will require an implementation of the interface conditions. A natural approach to take is that of front tracking, which requires that the interface velocity be imposed explicitly. Suppose we were to take the interface condition (4.28) as the condition which defines the velocity L(t). A n y explicit implementation of this will require evaluation of the indeterminate form given in (4.15). Thus, any front tracking scheme requires not only the explicit computation of interface location at each time step, but an accurate numerical evaluation of this limit. Furthermore, extension of a front tracking scheme to higher dimensions would also require consideration of the interface geometry, and solutions to problems on irregular domains. Clearly, front tracking is not feasible for this model problem. A n alternative is to reformulate the problem over the fixed domain 0 < z < D, thereby avoiding the need for explicit consideration of the complex interface physics. The interface location can be recovered, a posteriori from the solution of the transient, fixed domain problem. That is, we aim to develop an interface capturing method. Here, we present a reformulation based in part on the mixture model described by Wang and Beckermann [68]. The reformulation is in terms of a density over the entire pack, rather than saturation in just the two-phase zone. The main point is that if we  Chapter  4.  88  The M2 mixture method for the transient Udell problem  consider the water anywhere in the porous pack to be a mixture of vapour and liquid, then the density of this mixture must be continuous, even as we cross the interface. Suppose we. define the m i x t u r e d e n s i t y p by p = p s + p {l-s). l  (4.30)  v  Then equations (4.21)-(4.22) reduce to the system 4>Pt (pc)T  t  = =  ~(piui + p u )z v  1  v  KT -h (((j)p (l-s)) zz  vav  Now, the variables u ,ui,p v  v  Ui =  v  + {p u ) )  t  v  v  z  > J  (4.31)  are all functions of s and T, and in particular  —K , ( d —s (J-z(pv-SJ(s)) 6  +  (  py  4  3  2  )  (4.33)  l9  and • J(s) = 1.417(1 -s)  - 2 . 1 2 0 ( 1 -sf  + 1.263(1 - s ) , 3  (4.34)  where p is the vapour pressure. So if p ,p ,s can be found as functions of p and T , then (4.31) is a system in the two dependent variables p and T. Furthermore, given that p is the density of the liquid-vapour mixture, the system (4.31) is valid over the entire porous pack, rather than just the two-phase zone. Thus, we seek solutions to (4.31), from which we can recover the position of the interface between the two-phase zone and the vapour-only zone. Now that we are considering a system valid over the entire domain, we can not take the vapour pressure to be equal to the saturation pressure. Rather, the vapour pressure at a point will be equal to either the saturation pressure (in which case, the point is in the two-phase zone) or the pressure given by the ideal gas law for the vapour-only zone. Given values of p and T , a comparison of these two pressures, namely v  v  p (T) sat  = ae , bT  v  p*( ,T) p  =  ^  ,  determines whether or not the vapour is fully saturated, and hence the values of p ,p and 5. If p* < Psat, then the vapour is undersaturated, so must v  v  Chapter 4. The M2 mixture method for the transient Udell problem 89 be i n the vapour-only region. In this case, the vapour pressure p = p*. If Psat < P*i then the vapour is fully saturated, and hence p = p t- We see that equations (4.31), together with the algebraic constraint sa  p„ =  mm^ (T),^T)  (4.35)  o t  form a differential-algebraic system for the two variables p, T, which is valid over the entire domain. In the next section, we describe the numerical solution of this model problem, with p, T as the primary variables.  4.3  Computational method  Equations (4.30)-(4.35) give a coupled parabolic system on the domain 0 < z < D. Now we consider the discretization of the problem, written as  *-.iM*-s)-r  <4 36)  Tar + w(p,T) = ^- Q(p,T,^-\ t  T  (4.37)  z  where the fluxes q, Q are given by  * (?' ' ti) i (Me- ") T  =  1  ^SA(P.T)  + f (p,T) B  where f (p,T)=*s* A  Pl  + ^(l-s)\ Pv  9A(P,T)=P,  f (p,T) B  =  ^S,  Pi g {p,T)  = i>{s),  fc(p,T)  =  B  9c(p,T)=p.  p (l-s) , 3  v  y SB{p,T)\ Z  ,  (4.38)  Chapter  4.  The M2 mixture method for the transient Udell problem  90  and the function w is given by (4.40)  w(p,T) = h <f>p (l - s). vap  v  The time-dependent terms in the energy equation have been grouped together, leaving the time derivative of an enthalpy-type quantity. Also, in keeping with formulations of degenerate diffusion problems for numerical computation [25], we have rewritten the degenerate liquid velocity as (4.41)  U j = - — (s p + 5(i>(s)) ) , 3  z  z  where the function ip is given by iP(s)  =  -  'S  (ffffidt  =  0.2415s - 0.6676s + 0.6315s . 4  5  6  (4.42)  This change-of-variables idea is used in preference to regularisations of the type K i = s + n, such as those seen in the fuel cell literature [50, 51]. The numerical convenience offered by such regularisations is due to the fact that they smear out the sharp interface. Also, we have the boundary conditions 3  r  q\z=o = Qbot(t),  q\ =D = Qtapit), z  T|  z = = 0  = T (t), bot  -T\  Z=D  = T (t). top  (4.43)  For the closed sand pack, the boundary fluxes are zero, and we will take the boundary temperatures to be constant in time. Now we wish to implement a numerical scheme for the solution of the system (4.36)-(4.43). Two implementation options are available for explicit time-stepping. One option is to use a quasi-enthalpy method. The density p may be explicitly updated from the mass equation (4.36). Then the e n t h a l p y E, defined by E(p,T)  = pcT +  w(p,T),  (4.44)  may be explicitly updated from the energy equation (4.37). The temperature T must then be recovered from the enthalpy by way of a nonlinear solver. A second option is to form a mass matrix, and write the system (4.31) as (4.45)  Chapter  4. The M2 mixture method for the transient Udell problem  ghost cell  ghost cell  I "II  Pr 1, T  X  91  P' 2 . I *—,I — I * — 2  R T,  B ,T  Q , . ^  Q '^2  V  2  • Vi I W l | II II i•'•t — I 1*—I  Pj-1'Vl  , *—I  T  I  I  R ,T  2  M  Q  2  I *—i  P  I I I  R., Tj  H  i-r H  Q  q  N r N 1, N 2' N 2 T  +  P  +  R N +  r i q  T  +  +  *— x  q  r  T N +  1  N+V«N I +  Figure 4.3: G r i d and staggered grid. where O H = (p,  a  12  = 0,  together with a-21 =u h J.(n (p I (1 - \sP»j — d  vap  - p T^\ I , d  s  V  and a explicit = pc +methods h 4> ^(1 - s ) the ^ -computation Pv^j Either of these two requires of derivatives of all quantities with respect to the primary variables p, T, as would an implicit method. Given the inherent stiffness i n the problem, here we shall use a fully implicit scheme. Now we discretize the parabolic problem (4.36)(4.43) by finite volumes, in order to conserve the total mass, given by vap  2 2  f  p(z,t) dz,  10  Jo in a discrete sense. Given that mass flux is given at the boundaries (4.43), we develop a scheme which computes a discrete approximation to the solution of (4.36)(4.37) with mass fluxes on a grid which coincides with the boundary points, and density p on the corresponding staggered grid. T h a t is, we use a finite volume scheme for updating pj, which is the cell average density for cell j. Also, we use cell averaged temperatures. Consider Figure 4.3. Let Tj and Pj be the average values of temperature, T , and density, p respectively, over grid cell j. We introduce vectors T and R to represent the interpolated values of temperature T a n d density p falling on the grid. So we have Tj = T  j+1/2  =  2±2±i  " .for j = l , ..,#.+1, (4.46)  Rj=p  j+1/2  =  forj = l,..,AT + l .  Chapter 4. The M2 mixture method for the transient Udell problem 92 Then a fully implicit (Backward Euler), conservative, finite volume scheme for the mass equation (4.36) is given by the discretization o  — ri? 1  n+1  F, =  P j  > -  (  P k  Qj  -  q j  - )  n + 1  x  = 0,  j = 2,.., N+l,  (4.47)  where the discrete fluxes are given by +f (R ,T ) B  j  ^ B i P i + u r ^ - g B i p , , ^  j  ^ ^  (4.48) Next, we discretize the energy equation (4.37) in a similar manner. We have W  -  P°  k~  +  ,k  - ^ ( Q J - Q J ^  1  = 0,  j  =  2,..,N + l,  (4 49)  where  Q j  . _  hygpK  j r > T j + i - T j  h  f  / D  rp \ I  gc(Pj+1 iT>+l)~9C(pji j)\ T  T  ^ ^ J ' ^ V j = l,..,iV+l.  /»  j '  (4 50)  Here, fc, are the chosen time step and grid spacing, and superscripts denote the time level. To close the system for the 2N + 4 unknowns, Tj,pj, j = 1 : N + 2, we need four more equations. The boundary conditions give =  qi — qbot =  Gi = T i  — Tbot  o,  = o, (4.51)  FN+2 GN+2  = QN+I — qtop = = TN+I  Ttop  o,  = o,  Note that we cannot find the ghost values pi,pw+ explicitly in terms of interior cell values, so we just include them in the system to be solved. That is, we have 2(AT + 2) nonlinear equations 2  Fj = 0,  Gj = 0,  j = l,..,N  + 2,  Chapter  4.  The M2 mixture method for the transient Udell problem  93  for the unknowns Pj,Tj, j = 1,..N + 2. We use a Newton method with analytical Jacobian. In order to avoid "if" statements in the code, we replace the nonsmooth map (4.35) with the smoothed minimum p =  min (p ,p*) £  sat  =  H (p* - p )p e  sat  sat  + (1 - H (p* - Psat))p*, e  (4.52)  where the smoothed Heaviside function is given by  H (X) £  = ^ l  +  tanh(^pjY  (4.53)  for a smoothing radius e, and p* = -j^pT. In practice, the smoothing radius e can be taken arbitrarily small, and thus, the singularity and degeneracy in the saturation equation are retained.  4.4  Numerical results and discussion  In Figure 4.4, we show a typical result, starting from an initial uniform density, with no mass flux across the boundaries. The initial temperature is uniform, and the system is subject to sudden heating up to a fixed temperature at the upper boundary. A number of density, temperature and saturation profiles are shown, for increasing time, and we note that the profiles shown approach the correct steady-state profiles. The transition to an interface s = 0 is captured, and the interface between the two-phase and vapour regions is clearly moving to the left. In Figure 4.5, the interface position, L(t) is plotted. The stepwise behaviour is due to the discretization, and the fact that we take L(t) as the first grid point to the right of s = 0. Also, we plot L(t) for three grid sizes: N = 20,40,80. Convergence to a base numerical solution is clear as the grid is refined. The question remains, however, of whether the numerical method indeed computes the interface velocity accurately. This will be addressed in the next section. Figure 4.6 shows the temperature history at a point. A n initial increase in temperature is numerically the result of the sudden-heating boundary condition. After this, the temperature history has a wavy, stepwise nature, in common with the enthalpy-method solution of the Stefan problem [3]. For the results shown in Figures 4.4-4.6, small time steps are initially required,'in order for the Newton iteration to converge to within the specified tolerance. In fact, initially, we use Forward Euler time-stepping to get past  Chapter  4. The M2 mixture method for the transient Udell problem  94  Evolution to correct steady-state profiles, starting with uniform density and temperature  0  0.05  0.1  0.15  0.2  0.25  Figure 4.4: Evolution to correct steady-state. the initial transients, and then employ the fully implicit scheme. However, the time step required for convergence to within a fixed tolerance decreases as the saturation profile steepens before the interface moves by a grid point. In Figure 4.7, we demonstrate the stiffness in the problem. A n experiment was performed using implicit time stepping with adaptive time-step. Whenever the nonlinear solver is unable to converge to within the required tolerance, the time-step is halved. The figure clearly shows the relationship between the time-step and the interface advance.  Chapter  4. The M2 mixture method for the transient Udell problem  95  Moving interface location, with T0=360 T1=367.9727 W=13.8153 L=0.175  Figure 4.5: Interface location L(i), with grid refinement.  4.5  Similarity solutions and numerical convergence study  Numerical experiments using our mixture-based capturing method for the time-dependent problem show that initial distributions evolve to the correct steady-state solutions, as given by the disjoint-domain method of Chapter 3. In this section, we carefully examine the dynamics of the system to ensure that the method indeed captures the moving interface accurately in time. Typically, the convergence of numerical methods for moving boundary problems is demonstrated using initial and boundary data consistent with known similarity solutions. Similarity solutions are available for the two-phase Stefan problem and the porous medium equation. In the literature, analytical convergence studies of numerical schemes have appeared for reduced, scalar models of two-phase flow in porous media [2, 49]. Also, analytical solutions have been found for reduced and scalar models of phase change in porous media [9, 13, 34, 41, 42, 59, 60, 74]. Wang and Cheng [69] summarize analytical solutions to model problems presented by a number of authors, for both steady-state and transient problems. These have largely been constructed by neglecting temperature and capillary effects. Wang and Cheng [69] propose  Chapter  4.  The M2 mixture method for the transient Udell problem  364.51  360  1  0  1  1  Temperature history atz=0.218 1 1 1—  1  ' 1.5'  1  0.5  96  1  1  1  2 time t (seconds)  L  2.5  :  3  1  1  3.5 x  1Q  «  4  Figure 4.6: Temperature history at a point. a similarity solution for a full model including temperature and capillary effects, and underline the need for further analytical solutions to such problems. Here, we reduce our model problem only by simplifying the coefficients, but we retain the full, nonlinear, vector problem, and construct a travelling wave solution.  4.5.1  A reduced model problem  Let us consider a model problem with reduced Darc'y velocities, and the majority of coefficients equal to one. Specifically, we consider liquid flux = —pis (p + s ), 3  z  z  and  vapour flux = —p p , v  z  where we have simplified the Leverett function and the relative permeability, and have taken absolute permeability and viscosity to be unity. We note that removing the relative permeability of vapour, namely (1 — s ) , will not affect the singularity structure at s = 0, but it will allow for solutions with saturation s > 1. In this sense, the solutions cease to be physically meaningful. 3  Chapter  4. The M2 mixture method for the transient Udell problem M2 interface evolution and timestep  10'  T~I  1  r  1  1  r  97  r  i  10  Vi 10  E  10  10  8  10  9 x10  Figure 4.7: Timestep and interface location. Conservation of liquid mass reads (pis)t-  (pis (s +p ))  = T.  3  z  z  z  (4.54)  Conservation of vapour mass reads (pv)t - {pvPz)  z  = -r,  (4.55)  where the weighting (1 — s) has been taken as 1. Again, this will not affect the structure near s = 0. Conservation of energy is T =T t  zz  + T.  (4.56)  In the vapour region, we have T = 0,  p = p T (ideal gas), v  p <p t(T) sa  (undersaturated),  s = 0,  Chapter  4.  The M2 mixture method for the transient Udell problem  98  leaving a problem for p , T (or p, T). In the two-phase region, we have v  p = p T (ideal gas),  p = p t(T)  v  (fully saturated),  sa  leaving a problem for s,T and T. A t the interface z = L(t), we have five conditions, corresponding to (4.25)-(4.29):  (Mass) (Energy)  s = 0,  (4.57)  [T}_ = 0,  (4.58)  \p£ = 0,  (4.59)  pia-L{t)  = - ( s {s +p )  + p )~  3  Pl  z  z  p,s"L(t) = (-pis (s  +p ))~  3  z  + (p(pT) )  +  vPz  z  - [T }_ .  z  z  Note that, in view of the ideal gas law p = p T, further imply continuity of vapour density v  ,  (4.60) (4.61)  conditions (4.58)-(4.59)  [Pv]t = 0.  (4.62)  Also, for convenience later, we note that (4.60)-(4.6l) give 0 = -(p )vPz  + (p(pT) )  + [T }_.  +  z  z  (4.63)  Then, the five conditions (4.57),(4.58),(4.62),(4.61), (4.63) can be used to solve the problem. For the fixed domain problem on 0 < z < D, first we define the mixture density p = ps + p , l  (4.64)  v  and then use the "M2-map" p =  mm(pT,p (T)), 3at  n  =  £ •  r>v  —  j< i  (4.65)  Chapter  4. The M2 mixture method for the transient Udell problem 99  leaving the coupled system  (4.66)  Now, we consider forms of the saturation pressure function, p t(T) yield semi-analytical travelling wave solutions. sa  4.5.2  Travelling wave solution for case p t(T) sa  =  which  otT  Suppose the saturation pressure is a linear function of temperature, such that PsatiT) = aT. Then p = a in the two-phase region. Also suppose that we have an infinite porous medium —oo < z < oo, with far-field temperatures not yet specified. At time t, let there be an interface at z = at separating a two-phase region z < ct, and a vapour-only region z > ct. So the interface initially is at z = 0, that is, L(0) = 0, and we seek travelling wave solutions with speed c. Then, in the two-phase region, mass and energy conservation give v  (s (s 3  z  + aT )) z  (l + a )T 2  2 2  +  z  -Tzz 2  .  (4.67)  (4.68)  In the vapour region, mass and energy conservation give p =  (p( T) )  t  P  T =' Tzzt  z  z  z  (4.69) (4.70)  Finally, at the interface z = L(t), the 5 conditions (4.57), (4.58), (4.62),  Chapter 4. The M2 mixture method for the transient Udell problem 100 (4.61), (4.63) give s= 0  (4.71)  [T]! = 0  '  . • M t =0,  (so pt = <*),  3  Pl  (4-73)  + aTz)Y - [T ) _ ,  s-c = (- s (s  Pl  (4.74)  +  z  z  0 = - (a T )~ + ( ( T) ) 2  + [T } _ .  +  z  P  P  (4.72)  +  Z  z  (4.75)  We seek travelling wave solutions to the system (4.67)-(4.70), which satisfy the interface conditions (4.71)-(4.75). In the two-phase zone, seek solutions s(z,t) = G(Q,  T(z,t) = F (£),  (4.76)  T(z,t) = F (0,  (4.77)  1  and, in the vapour zone, seek solutions '  p(z,t) = R(t),  2  where i = z-ct.  (4.78)  W i t h c > 0 we have wavefronts moving to the right with speed c. Solutions of the two linear heat equations are F (0  = A + B e-^S,  1  1  (4.79)  1  and F (0 2  = A + B e-'*. 2  (4.80)  2  The constants A ,B ,A , B must be consistent with the interface conditions. Continuity of temperature (4.72) gives 1  1  2  2  A +B = A +B. 1  1  2  (4.81)  2  Now, the saturation equation (4.67) gives -cG'=(G (G' 3  + aFl))' + -F[', Pl  Chapter  4. The M2 mixture method for the transient Udell problem  and hence  a + —F[ + A , Pi  101  2  -cG = G {G'-r-aF[) 3  or G G' = -(cG  + a (G + j j  3  (4.82)  3  Fi + A ^j .  3  (4.83)  3  The constant A must be consistent with the interface conditions. Condition (4.74) gives 3  ss 3  z  G G" = - (cG + < * G ( — \ 1+ a Comparing with (4.83), we find 3  3  2  A  - ( T - ^ ) ^ - -cB ) . Pi 1 + a Pi / 2  2  = -(B -B ).  3  1  (4.84)  2  Then the travelling wave saturation profile G(£) may be found by solving the O D E (4.83), subject to G(0) = 0, for f < 0. The constant A is given by (4.84), and (4.79) gives 3  _c _ t  Notice that we can find an analytical solution to (4.83) if B\ = 0, which amounts to an isothermal two-phase region, with no variation in vapour pressure, and hence no vapour flow. In this special case, the exact solution G(£) is given implicitly by G + B  3  -G - -B G +B G-B \og  l  3 l  2  2  3  3  B  +c£ = 0,  3  where B = -—, Pi 3  (4.85)  and G(£) may be accurately computed using a root finder. In general, though, an analytical solution to (4.83) will not be available, and we compute &(£) using an O D E solver. However, equation (4.83) is singular at G = 0. Letting w = G , we have the regular problem 4  ^w' = - (cw  1/A  + a (w * 3/  + j^j  F[ + A ^j, 3  w(0) = 0.  (4.86)  Chapter  4. The M2 mixture method for the transient Udell problem  102  It remains to find the density p in the vapour-only region. Equation (4.69) becomes -cR' = (R(RF )')', 2  and hence -cR = R{RF )'  + A.i.  2  (4.87)  Now, to find the constant 4 , we use condition (4.75) to give that, at £ = 0, 4  R(RF )'  =  2  {l + a )F[-F^  =  2  -cB +cB , x  2  and condition (4.73) to give —cR = —ca. Substituting in (4.87), we find A  = c(B  4  From (4.87), we can write . . R' = -  -B -  1  a).  2  (4.88)  .  ^ { { c + RF! )R + AA,  1  )  (4.89)  1  and then solve for i?(£), in general, using an O D E solver. Also, we note that, for the special case Bi = B + a, then A\ = 0, and the exact solution R(£) is given explicitly by 2  fl(fl=  ~ ^ * ^  5  ,  where A = a(A 6  1  + B ).  (4.90)  1  Solution - Summary We have constructed the following travelling wave solution. F (£)=A 1  F (0 2  + B e-&*,  1  (4.91)  1  = A + B e"*: 2  . ' (4.92)  2  and G ( £ ) = i v / ^ ) , where w solves 1  4  - w' = - [cw '  l A  A  + a (w '  3 4  + j^j F[ + A^j,  and finally, i?(£) solves R' = - ^ 7-A(c + RF>)R + A ), 1  7  i  £<0,  w(0) = 0, (4.93)  £>0,  R(0) = a.  (4.94)  Chapter  4. The M2 mixture method for the transient Udell problem  103  The constants A . A? A4 are given by 2  A  2  n  = AX + BX-B2,  A = - (B, - B ), 3  A = c(B -B -a).  2  4  1  (4.95)  2  Pi  We have a three-parameter family of solutions to the model Udell problem. That is, a family of solutions parameterized by Ai,Bi,B . These can be carefully chosen in order to have temperature increasing in both regions, and the density decreasing to the right of the interface. A typical result is shown in Figure 4.8. Now we are i n a position to compare numerical results with the travelling wave solution. 2  Travelling wave profiles for o=0.7 p =5 A,=1 B^-0.5 B =-0.83557 2  -0.5  -0.4  -0.3  -0.2  -0.1  0  0.1  0.2  0.3  0.4  0.5  rj=z-ct  Figure 4.8: Travelling wave profiles for mixture density, temperature and saturation. N o t e : W i t h B\ < 0 to ensure that F\ is increasing, we have F\ —> —oo as £ — » —co. Clearly, the temperature i n the two-phase region will become negative for some value of t, and remain negative as t increases. Our code takes the M 2 comparison as p = min(p, a). This gets around the fact that the vapour pressure becomes negative as temperature does. So, while comparison with this travelling wave solution serves the purpose of evaluating the performance of the numerical method, we note that the v  Chapter  4.  The M2 mixture method for the transient Udell problem  104  t r a v e l l i n g wave solutions t o this reduced p r o b l e m o n l y represent a n y t h i n g close t o the p h y s i c a l p r o b l e m u p to the t i m e at w h i c h T a n d p become negative. T h e s o l u t i o n s t r u c t u r e of the p h y s i c a l p r o b l e m at the interface is retained, though.  4.5.3  Numerical results and convergence study  In order t o evaluate the performance of o u r c a p t u r i n g m e t h o d a p p l i e d t o the reduced U d e l l p r o b l e m , we s i m p l y give our code i n i t i a l c o n d i t i o n s a n d b o u n d a r y c o n d i t i o n s consistent w i t h the exact t r a v e l l i n g wave s o l u t i o n given b y (4.91)-(4.95), for the chosen values of the parameters Ai, Bi, B . I n F i g ure 4.9, we show c o m p u t e d interface l o c a t i o n , n o t i n g close agreement w i t h the exact t r a v e l l i n g wave solutions, a n d apparent convergence w i t h g r i d refinement. N o w we w i l l d e m o n s t r a t e the convergence of o u r schemes n u m e r i c a l l y . 2  Interface position for reduced Udell problem with travelling wave solutions 0.12  1  — —  0.1  :  1  :  1  :  :  1  :  r  N=160 N=320 exact TW solution  c=4  -1 0.08  I  0.04  -  0.02  -  0  0.005  0.01  time, t  0.015  0.02  0.025  F i g u r e 4.9: C o m p u t e d a n d exact interface l o c a t i o n for the reduced t r a v e l l i n g wave p r o b l e m .  Chapter  4. The M2 mixture method for the transient Udell problem  105  Although we have a second order scheme and an exact solution, we don't expect to be able to show second order accuracy. Firstly, we need an accurate measure of interface location. A s mentioned in Chapter 2, methods are available for the enthalpy solution of the Stefan problem [65], but are based on a known phase-change temperature, and will not be applicable here. In the results shown, the interface location is simply the next grid point to the right of s = 0. Also, the stepwise behaviour of not only the interface location, but of the temperature history at a point, introduces errors, which result in less than second order accuracy (see [3], for example). Even i f novel methods are used to capture the interface position accurately, the stepwise temperature histories associated with the enthalpy method can give rise to large pointwise errors in temperature. Indeed, the norm of the temperature error is a poor measure of accuracy, and though the L\ norm is a better measure, convergence in this norm is at first order. Bearing in mind these results, we expect similar, or worse, behaviour for our capturing method applied to the reduced Udell problem, and no better than first order accuracy in the solutions. In the Tables .4.1-4.3, we show the results of a numerical convergence study for the finite volume, capturing method, applied to the reduced Udell problem. Table 4.1: Errors for reduced Udell problem, implicit ( B E ) time-stepping with p, = 0.2, c = 4.  N  num  20 40 80 160 320  50 200 800 3200 12800  ll-Erlli 5.2534 1.8412 6.7618 2.6595 1.1249  E-3 E-3 E-4 E-4 E-4  factor 2.85 2.72 2.54 2.36  \\E h P  9.6053 E-2 4.1025 E-2 1.7308 E-2 7.2590 E-3 3.0324 E-3  factor 2.34 2.37 2.38 2.39  \\EL\U 1.4120 E-2 7.3350 E-3 3.8450 E-3 2.0283 E-3 1.0625 E-3  factor 1.93 1.91 1.90 1.90  We have exact solutions given by the travelling waves, with speeds c = 1, 2 and 4. In each case, we use implicit time-stepping, with p, = k/h? fixed. The number of grid points is AT, and num is the number of time steps. We quantify the errors by calculating numerically the norms of the errors made in the temperature T , the mixture density p, and the interface location L.  Chapter  4. The M2 mixture method for the transient Udell problem  106  Table 4.2: Errors for reduced Udell problem, implicit ( B E ) time-stepping with /j 0.2, c = 2. N  num  20 40 80 160 320  50 200 800 3200 12800  || F r l l i 1.9986 E-3 7.0790 E-4 2.5842  E-4  9.7712  E-5  3.8206 E-5  factor 2.82 2.74 2.64 2.56  \\E h P  8.1047 3.3881 1.4097 5.8601 2.4351  E-2 E-2 E-2 E-3 E-3  factor  II EL ||I  factor  2.39 2.40 2.41 2.41  1.3620 E-2 7.3350 E-3 3.8450 E-3 2.0283 E-3 1.0625 E-3  1.91 1.90 1.89 1.91  Table 4.3: Errors for reduced Udell problem, implicit ( B E ) time-stepping with fj, — 0.2, c = 1. N 20 40 80 160 320  II F r | | i  num 50  1.1646 E - 3  200  3.6319 E - 4  800 3200 12800  1.2909 E-4 4.5586 E-5 1.6463 E-5  factor  \\E \\i  factor  3.21 2.81 2.83 2.77  5.9305 E-2 2.6721 E-2 1.1080 E-2 4.6023 E-3 1.9151 E-3  2.22 2.41 2.41 2.40  P  \\EL\\I  1.5220 E-2 7.0375 E-3 3.7072 E-3 1.9445 E-3 1.0119 E-3  factor 2.16 1.90 1.91 1.92  Specifically, we compute the time-averaged quantities defined by ..  - num N n=l j=l  ^ l l ^ "  1  ^ num N S Yl \P^ J^) 71=1 ] = 1 Z  =  ibLLlIb i v  -  Pexact(Zj,  t )\, n  num |£L||I =  — mum num  JZlHQ-ct, n=l  In all three cases, we see convergence at first order. We note here that, when computing solutions to the two-phase Stefan problem with travelling wave  Chapter  4. The M2 mixture method for the transient Udell problem  107  solutions, using the enthalpy method described in Chapter 2, we found the temperature error to decrease by a factor of around 1.9-2.5 as the grid spacing is halved. Also, applying the discretization used by Evje and Karlsen [25], applied to the one-dimensional porous medium equation u = (u u ) , we found the error in u to decrease by a factor of 2.4-2.5 as the grid spacing is halved. Our capturing method for this coupled, vector, nonlinear problem clearly gives convergence rates comparable with those for these scalar prototype problems. z  t  4.5.4  Other choices for  x  x  p t(T) sa  We briefly mention two more forms for the saturation pressure which yield semi-analytical travelling wave solutions which may be of interest. C a s e (T) = aT + (3 The temperature in the two-phase region, jF\(f) may be found analytically, but it is given implicitly. C a s e p (T) = aT + (3T In this case, we have p {T) = a + (3T, which is increasing if a, (3 > 0. The travelling wave solutions in the vapour region remain the same. In the twophase region, an analytical solution is again available for the temperature profile. For T(z, t) = F \ ( z — ct), we find that Fi(£) satisfies Psat  2  sat  v  -c{l  + (3)F = (l + (a + (3F )(a + 2(3F ))F{ + A , 1  1  1  1  (4.96)  with Ai yet to be specified. Thus Fi(£) is given implicitly by ^{Fr + df + (a - 2d)(F + d) + (d - ad.+ b) log |*\ + d| + 2  x  £=  B  u  (4.97) for some Bi, where the constants a, b, d are given by _ 3a ~2jV  a  4.6  1 + o? ~2W'  ,  A  1  c{l + (3)\  (4.98)  Computations in higher dimensions  One of the distinct advantages of a fixed-domain, front-capturing method over front-tracking is that extension to higher dimensions is much more  Chapter  4. The M2 mixture method for the transient Udell problem 108  straightforward. In this section, we return to the full mathematical model for the physical problem, and present some computational results i n two dimensions.  4.6.1  Mathematical model  If we generalize the governing equations (4.1)-(4.3) to higher dimensions, we have ( c M s ) + V.(flu,)-=r, (4.99) t  where, as before, u represents the Darcy velocity, and T is the condensation rate. Similarly, conservation of mass for the vapour phase is given by (# (l-s)) l )  t  + V.(p„u,) = - r .  (4.100)  The energy equation is now (p-c)T = V.(kVT) t  + h r,  (4.101)  vap  where the Darcy velocities are functions of the saturation and the vapour pressure, and are given by u„ = - ^ V p ,  (4.102)  fly  and u, = - ^ ( V p - W ' ( s ) V s ) ,  (4.103)  Pi  where p is viscosity, K is the permeability of the medium, K ^ are the relative permeabilities, and J(s) is the Leverett function. For our sample 2D computations, on Q = (0, D) x (0, D), we take periodic conditions in the x direction. A t the bottom boundary y = 0, take a given temperature, and no mass flux, so that R  T = T (x, t), 0  (pmi + p u ).n v  RV  = 0.  v  (4.104)  A t the top boundary y = 1, we take a higher given temperature, and no mass flux, so that T = T {x,t)(>T (x,t)), x  0  (piui + p u ).n v  v  = 0.  (4.105)  Now, given that there is no mass flux across the boundary, our initial condition fixes the total water mass W. Now, i n order to compute solutions to this problem, we simply extend the finite volume scheme described in Sections 4.2-4.3 to the two dimensional problem.  Chapter  4.6.2  4. The M2 mixture method for the transient Udell problem  109  Numerical results and discussion T1=513.3, T0=360, q=1500.s0=0.5, L=0.15  0  0.05  0.1  0.15  0.2  0.25  Figure 4.10: Steady-state, one-dimensional profiles. The first results we show use the steady-state results for the one-dimensional problem shown i n Figure 4.10. The total water mass (per unit area) in the system is W = 35.7, and the interface is L = 0.15. Now, we consider the two-dimensional, transient problem with uniform boundary temperatures at y = 0, D corresponding to those at z = 0, D for the one-dimensional problem. Now we give an initial condition with a fixed water mass (per unit length) of WD. Imposing no mass flux across the lower and upper boundaries, and giving periodic conditions in x, we expect the solutions to evolve towards the one-dimensional, steady-state results shown. Giving a one-dimensional initial condition gives the required result, showing agreement with the one-dimensional transient solutions at all times. In Figures 4.11 and 4.12, we show saturation, temperature and mixture density plots at increasing times, for a two-dimensional initial condition which has liquid water at a uniform saturation concentrated near the lower boundary in a block as shown, and vapour in the block above it. The temperature is taken initially uniform, so that T(x, y, 0) = Tb t = 360, and then the system is subject to sudden heating at y = D. A fully implicit scheme is used, with a 20 x 20 grid. 0  Chapter  4. The M2 mixture method for the transient Udell problem  110  We note the evolution of the system towards the correct steady-state. In Figure 4.12, we can see that after a time of about 20 minutes, the interface is approximately at y = 0.15. The influence of the initial condition is still apparent, but the saturation and temperature plots have roughly the familiar shape of the steady-state one-dimensional profiles. In Figure 4.13, we plot the liquid and vapour fluxes at early and long times, on a saturation contour plot. At t = 16, the liquid fluxes are small everywhere except near the interface. At the later time, there is a larger liquid flux at the lower boundary, where the vapour condenses, and then is driven upwards by the capillary pressure gradient. The vapour flux in the vapour region is initially large, but decreases in time as a steady-state is approached. The one-dimensional steady-state solution has stationary vapour in the vapour region. As a steady-state is approached, the vapour flux becomes essentially one-dimensional. In Figure 4.14, we plot contours of saturation at increasing time, for an initial condition with a "blob" of two-phase fluid at uniform saturation in the centre of the domain. Again, we have sudden heating at the top, closed upper and lower boundaries, and periodic conditions in x. We see  Chapter  4.  The M2 mixture method for the transient Udell problem  1=0  111  t=1115  Figure 4.12: Saturation and temperature profiles - initial condition and long time. the two-phase region spreading, and migrating towards the lower boundary. There is initially a large downward vapour flux in the vapour region (outside the blob), and condensation occurs at the cold boundary y = 0. A second two-phase region thus appears near the lower boundary, creating a second interface which is clearly visible at t = 189. Liquid accumulates in the lower two-phase region, and the blob migrates downwards through the surrounding vapour region, until the two separate two-phase regions coalesce, leaving just one interface once more. Our capturing method encounters no problems upon these topological changes. By t = 6709, it appears that the solution is almost one-dimensional. The total mass (per unit length) in the system in Figure 4.14 is 3.7425, and the uniform boundary temperatures are T = 360 and Ti = 512.7. Thus, the steady-state to which the system evolves should correspond to the onedimensional steady-state with mass (per unit area) W = 3.7425/D = 14.734. Using our disjoint-domain method from Chapter 3, this gives s — 0.2, q = 1000, and the interface position L = 0.104. In Figure 4.15, we plot saturation and temperature profiles for a fixed x in the two-dimensional problem, at a large time. Clearly, the system has almost reached the correct one0  0  Chapter  4. The M2 mixture method for the transient Udell problem 112  Saturation contours and liquid water flux at time t-16 0.25 f ' '. . .  Saturation contours and vapoui water flux at time t=16 0.25  0.2 0.15 0.1  TTTTTTTTt \ % i t i i i i i t i i .  0.05  0.1  0.2  Saturation contours and liquid water flux at time t=1115 Saturation contours and vapour water flux at time t=1115 0.251 ' .. .1 0.25 r  0.15 0.1 0.05  0.1  0.1  0.2  0.2  Figure 4.13: Saturation contours with liquid and vapour fluxes. dimensional steady-state. In our final numerical experiment, we again give periodic conditions i n x, no mass flux across the lower and upper boundaries, an initially uniform temperature distribution, at To = 360, and sudden heating at the upper boundary with a uniform temperature T i = 627.8. This time, we have an initial condition with total water mass of 7.5466, which corresponds to a onedimensional problem w i t h W = 7.5466/Z) « 29.7. The corresponding onedimensional steady-state solution has s = 0.58, 'q = 2000 and L = 0.12. The initial condition this time has two two-phase regions at uniform saturation, one i n the shape of a rectangle, and one i n a circle. 0  In Figure 4.16, we plot contours of saturation at increasing time. Our method again captures the topological changes with no difficulty. The two two-phase regions spread, while a third two-phase region develops through condensation at the cold lower boundary. The two initial two-phase regions  Chapter  4. The M2 mixture method for the transient Udell problem 113 t=o  t=64  t=189  t=498  t=1435  t=5612 0.2  t=6709 0.2  0.1  -_:.-_.ilI.I.^-.::. .'::^-J :  0 0  0.1 0.2 x  wm—mm  0 0  0.1 0.2 x  F i g u r e 4.14: S a t u r a t i o n contours for a n i n i t i a l " b l o b " of two-phase fluid i n the centre o f t h e d o m a i n . coalesce, a n d also t h e i n i t i a l l y c i r c u l a r region coalesces w i t h t h e two-phase region near t h e lower b o u n d a r y .  E v e n t u a l l y , t h e r e m a i n i n g v a p o u r region  near t h e lower r i g h t h a n d corner of the d o m a i n disappears, l e a v i n g just one two-phase region, a n d a o n e - d i m e n s i o n a l steady-state is approached. I n F i g ure 4.17, w e p l o t t h e s a t u r a t i o n a n d t e m p e r a t u r e profiles for a fixed x i n this t w o - d i m e n s i o n a l p r o b l e m , at a large t i m e . C l e a r l y , this system has also almost reached i t s correct one-dimensional steady-state.  Chapter  4. The M2 mixture method for the transient Udell problem  0.25  r  Comparison ol blob solution after long time, and 1 -D steady state 1 ; n 1 1  114  r  y  Figure 4.15: Long time saturation and temperature profiles at fixed x for the spreading, migrating blob, together with one-dimensional, steady-state solution.  Figure 4.16: Saturation contours for an initial condition with multiple twophase regions in the centre of the domain.  Chapter  4. The M2 mixture method for the transient Udell problem  Comparison of circle and rectangle solution after long time, and 1-0 steady state 1 1 1 -i  0  0.05  0.1  0.15  0.2  116  r  0.25  y  Figure 4.17: Long time saturation and temperature profiles at fixed x for initial multiple two-phase zone problem, together with one-dimensional, steadystate solution.  117  Chapter 5 Conclusions and future work 5.1  Summary of results  In this thesis, we have described a number of free and moving interface problems, related to phase change processes in porous media. Throughout the thesis, we have made note of mathematical points of interest, with a view to motivating further work. The first result of this work is an asymptotic analysis of smoothing methods applied to the one-dimensional, steady-state, two-phase Stefan problem. This appears to be new work. Smoothing techniques are often implemented in capturing methods for moving interface problems, with little mention of the effect they have on computational results. In fact, Crank [20] mentioned that very little analysis of smoothing for the enthalpy function used in the enthalpy method had appeared in the literature, and this still appears to be the case. Also in Chapter 2, we have presented numerical convergence studies for the well established enthalpy method for the Stefan problem, as well as a more recently described method for the Porous Medium Equation [25]. This has been done with a view to drawing a direct comparison between the convergence rates for the capturing methods for these scalar problems, and the capturing method that we develop for the vector problem of phase change in porous media in Chapter 4. In Chapter 3, we have modified a recently developed method for solving the one-dimensional, steady-state free boundary problem of a three zone system to be applied to a two-zone system of particular interest. The model problem has been extended to allow for compressible vapour. Numerical results have been generated using this iterative disjoint-domain method. Also in Chapter 3, the method of residual velocities [22] for solving steady free interface problems with linear, scalar elliptic problems on either side of a free interface, has been extended to solve the nonlinear, vector problem of phase-change in porous media. Numerical results show agreement with the  Chapter  5. Conclusions and future work  118  results from the iterative disjoint domain method. In Chapter 4, we have developed a numerical interface capturing method for a model problem of time-dependent two-phase flow with phase change in porous media. We have allowed for a nonisothermal two-phase region, rather than making the popular assumption of an isothermal region. This will allow for computations in particular settings where thermal effects are important, such as fuel cell electrodes. Numerical solutions found using this capturing method are seen to converge in time to the correct steady-state solutions generated using the method of Chapter 3. These benchmark steady-state solutions have been particularly useful in v a l i d a t i n g the results of the capturing method. We have found similarity solutions to a model vector problem for phase change in porous media. The singularity, degeneracy and coupling in the mathematical model have been retained, and we have demonstrated convergence of our numerical results to the exact travelling wave solutions, thus showing that the capturing method recovers an interface moving at the correct velocity. To our knowledge, this is the first such convergence study for a full, coupled model. Convergence rates are comparable with those presented in Chapter 2 for the simpler, scalar problems. The implementation of the capturing method has been extended to the two-dimensional problem, and we have demonstrated that our method can easily deal w i t h m u l t i p l e interfaces, and topological changes such as disappearing interfaces. This is a distinct advantage over front-tracking methods. Also, no grid refinement near the interface is required by our method, which is an advantage over a recently described numerical capturing method for a related problem in the fuel cell literature [10]. Convergence to the correct steady-state solutions is once again demonstrated using our benchmark solutions.  5.2  Future work  We suggest using the ideas presented in Chapter 2 as part of an asmymptotic analysis of smoothing techniques applied to the time-dependent Stefan problem. This will involve the analysis of smoothing strategeies applied to the discontinuous enthalpy functions, and possibly an accompanying analysis of the relationship between the smoothing parameters and the computational grid spacing.  Chapter  119  5. Conclusions and future work  • In Chapter 3, the method of residual velocities has been used for the onedimensional porous medium problem. Future work which we have suggested includes a n i m p l e m e n t a t i o n of the method applied to the two-dimensional. The analysis of the residual velocity method has so far been applied to problems of Laplacian type [19, 22], and we suggest using the model problem described in Subsection 3.3.2 as a starting framework for the analysis of the full Udell problem. The capturing method that we have developed in Chapter 4 appears to be robust, and no grid refinement has been necessary in our computations. We have, however, used an adaptive time-stepping method. Future work will involve looking at strategies for this adaptive time stepping. Our computational method has been implemented for the so-called Udell problem of phase change in porous media. The model problem which we have considered is quite general, and future work will include applying this capturing method to specific physical and industrial problems of interest. In particular, we will look at model problems of phase change in fuel cell electrodes and oil reservoirs. Also, we will consider modified models of phase change which relax the saturation assumption which we have used throughout this thesis. We briefly describe such a model problem below, and some ideas concerning its solution.  5.2.1  The " B i g - H " regularisation  We take this idea from the fuel cell modelling literature (see, for example, [10, 21, 50]), a n d refer the reader to Wang [67], who discusses the computational convenience offered by this regularisation. The vapour is not assumed to be at saturation in order for phase change to occur. Instead, any oversaturated vapour is supposed to condense at a very large rate, proportional to the degree of oversaturation. Undersaturated vapour may not condense. In the two-phase region, if vapour is undersaturated, then evaporation of liquid water occurs at a large rate, proportional to the degree of undersaturation. O n the other hand, i f vapour is undersaturated and no liquid is present, then no phase change occurs. We take the condensation rate to have the following form: / H (p - p (T)) p > p (T) , \ H-s(p -p (T)) p < sat(T) • -> +  v  sat  v  sat  (  {  v  sat  v  P  Here, H and H~ are large constants, which account for the large condensation rates. The factor of s precludes the possibility of evaporation if there +  Chapter  5. Conclusions and future work  120  is no liquid present. Now the governing equations for the model problem are the system (4.31), together with the condensation rate (5.1). These constitute three equations for three unknowns s, T and p on the fixed domain 0 < z < D, rather than the two equations in two unknowns when we use the saturation assumption. This idea has been used in the fuel cell literature, but computations have only been performed under simplifying assumptions, such as isothermal conditions, or with computational regularisations or local grid refinement [10]. The method we have described in this thesis does not require any of these simplifications. The relation between our model and this " B i g - H " model in the limit of large H , H~ is of interest to us, and in particular, how the size of these parameters affects the computations. We propose to implement this method for the reduced model problem described in Subsection 4.5.1. A n asymptotic analysis of this method will appear in future work, where we take H — H~ = l/e. Preliminary analysis looking for travelling wave solutions suggests a corner layer will be introduced which will smear out the sharp interface, but that the sharp interface solution will be approached as e —> 0. This is continuing work. +  +  121  Bibliography [1] M . M . Abbott and H . C . Van Ness. Theory and Problems of Thermodynamics (2nd edition), (1989), McGraw-Hill. [2] M . Afif and B . Amaziane, O n convergence of finite volume schemes for one-dimensional two-phase flow in porous media. J.. Comp Appl Math, 145, (2002), pp. 31-48. [3] V . Alexiades and A . D . Solomon, "Mathematical Modeling of Melting and Freezing Processes". Hemisphere Publishing, (1993). [4] M . B . Allen III, Numerical modelling of multiphase flow in porous media. Adv. Water Resources, 8, (1985), pp..162-187. [5] M . B . Allen III, G . A . Behie and J . A . Trangenstein, "Multiphase Flow in Porous Media". Springer-Verlag. [6] K.^Aziz and A . Settari, "Petroleum Reservoir Simulation". Applied Science Publishers, London, (1979): [7] P. Baggio, C . Bonacina and B . A . Schrefier. Some considerations on modeling heat and mass transfer in porous media. Transport in Porous Media, 28, (1997), pp. 233-251. [8] J . Bear, "Dynamics of Fluids in Porous Media". American Elsevier, (1993). [9] J . Benard, R^ Eymard, X . Nicolas and C . Chavant, Boiling in porous media: model and simulations. Transport in Porous Media, 60, (2005), pp. 1-31. [10] E . Birgersson, M . Noponen and M . Vynnycky, Analysis of a two-phase non-isothermal model for a P E F C . J. Electrochem Soc, 152 (5), (2005), pp. A1021-A1034.  122 [11] R. Bradean and D . B . Ingham, Stefan Problem in Blow Moulding Operations. Math. Heat Transfer, (1998), pp. 89-96. [12] 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). [13] L . Brevdo, R. Helmig, M . Haragus-Courcelle and K . Kirchgassner, Permanent fronts i n two-phase flows in a porous medium. Transport in Porous Media, 44, (2001), pp. 507-537. [14] L . Bridge, Condensation in a Porous Medium. M S c thesis, Dept. of Mathematics and Institute of Applied Mathematics, University of British Columbia, Vancouver, Canada (2002) 51pp. [15] L . Bridge, R. Bradean, M . J . Ward and B . R. Wetton, The analysis of a two-phase zone with condensation in a porous medium. J. Eng Math, 45, (2003), pp. 247-268. [16] J . Braining, D . Marchesin and C . J . van Duijn, Steam injection into water-saturated porous rock. Comp and Appl Math, 22 (3), (2003), pp. 359-395. [17] G . Caginalp, Stefan and Hele-Shaw type models as asymptotic limits of the phase-field equations. Physical Review A, 39 (11), (1989), pp. 58875896. [18] S. Chen, B . Merriman, S. Osher and P. Smereka, A simple Level Set method for solving Stefan problems. J. Comp Phys, 135, (1997), pp. 829. [19] W . Chen and B . Wetton, Residual velocities in steady free boundary problems of vector Laplacian type. P R E P R I N T (2006). [20] J . Crank, "Free and Moving Boundary Problems". Oxford University Press, (1984). [21] J . Divisek, J . Fuhrmann, K . Gartner and R. Jung, Performance modeling of a direct methanol fuel cell. J. Electrochemical Soc, 150 (6), (2003), pp. A811-A825.  123 [22] R. D . Donaldson, Generalised Stefan Problems - Linear analysis and computation. M S c thesis, University of British Columbia, (2003). [23] C . M . Elliott and J . R. Ockendon, "Weak and variational methods for moving boundary problems", Pitman, (1982). [24] A . Esen and S. Kutluay, A numerical solution of the Stefan problem with a Neumann-type boundary condition by enthalpy method. Appl. Math. Comput, 148, (2004), pp. 321-329. [25] S. Evje and K . H . Karlsen, Monotone difference approximations of B V solutions to degenerate convection-diffusion equations. SIAM J. Num. Anal, 37 (6), (2000), pp. 1838-1860. [26] S. Evje, K . H . Karlsen, K . - A . Lie and N . H . Risebro, Front tracking and operator splitting for nonlinear degenerate convection-diffusion equations. Parallel Solution of Partial Differential Equations, IMA Volumes in Mathematics and its Applications 120, (2000), pp. 209-228 [27] M . Fabbri and V . R. Voller, The phase-field method in the sharpinterface limit: A comparison between model potentials. J. Comp Phys, 130, (1997), pp. 256-265. [28] I. Fatt and W . A . Klikoff. Efffect of fractional wettability on multiphase flow through porous media. A I M E technical note #2043. AIME Transactions, 216, (1959), pp. 246. [29] C . Figus, Y . Le Bray, S. Bories and M . Prat, Heat and mass transfer with phase change in a porous structure partially heated: continuum model and pore network simulations. Int. J. Heat and Mass Trnsfr, 42, (1999), pp. 2557-2569. [30] A . C . Fowier, "Mathematical Models in the Applied Sciences". C a m bridge University Press, (1998). [31] R. M . Furzeland, A comparative study of numerical methods for moving boundary problems. J. Inst. Maths Applies, 26, (1980), pp.411-429. [32] J . Gratton and C . Vigo, Evolution of self-similarity, and other properties of waiting-time solutions of the porous medium equation: the case of viscous gravity currents. Euro J Appl Math, 9, (1998), pp. 327-350.  124 [33] K . Hanamura and M . Kaviany, Propagation of condensation front in steam injection into dry porous media. Int. J. Heat and Mass Trnsfr, 38, (1995), pp. 1377-1386. [34] K . T . Harris, A . Haji-Sheikh and A . G . A g w u Nnanna, Phase-change phenomena in porous media - a non-local thermal equilibrium model. Int J. Heat and Mass Trnsfr, 44, (2001), pp. 1619-1625. [35] W . He, J . S. Y i and T . V a n Nguyen, Two-phase flow model of the cathode of P E M fuel cells using interdigitated flow fields. AIChE Journal, 46 (10), (2000), pp. 2053-2064. [36] T . Y . H o u , Z . L i , S. Osher and H . Zhao, A hybrid method for moving interface problems with application to the Hele-Shaw flow. J. Comp Phys, 134, (1997), pp. 236-252. [37] R. Juanes, A variational multiscale finite element method for multiphase flow in porous media. Finite Elements in Analysis and Design, 4 1 , (2005), pp. 763-777. [38] K . Hyistendahl Karlsen, K . - A . Lie and N . H . Risebro, A fast marching method for reservoir simulation. Comp Geosci, 4 (2), (2000), pp. 185206. [39] S. K i m , S. Anghaie and G . Chen, A fixed-grid two-phase numerical model for convection-dominated melting and solidification. P R E P R I N T . [40] D . A . K n o l l , D . B . Kothe and B . Lally, A new nonlinear solution method for phase-change problems. Num. Heat Trnsfr, Part B , 35, (1999), pp. 439-459. [41] I. Kocabas, Thermal transients during nonisothermal fluid injection into oil reservoirs. J. Petroleum Sci Eng, 42, (2004), pp. 133-144. [42] D . Kulasiri and I. Woodhead, O n modelling the drying of porous materials: analytical solutions to coupled partial differential equations governing heat and moisture transfer. Mathematical Problems in Eng, 2005 (3), (2005), pp. 275-291. [43] A . A . Lacey, J . R . Ockendon and A . B . Tayler, "Waiting-time" solutions of a nonlinear diffusion equation. SI AM J Appl Math, 42 (6), (1982), pp. 1252-1264.  125 [44] R. J . LeVeque, "Numerical Methods Birkhauser, (1990).  for  Conservation  [45] M . C . Leverett. Capillary behaviour in porous solids. AIME tions, 142, (1941), pp. 152-169.  Laws".  Transac-  [46] Z. L i , Immersed interface methods for moving boundary problems. Numerical Algorithms, 14, (1997), pp. 269-293. [47] W . L i u , S. Shen and S. B . Riffat, Heat transfer and phase change of liquid i n an inclined enclosure packed with unsaturated porous media. Int. J. Heat and Mass Trnsfr, 45 (2002), pp. 5209-5219. [48] I. Lunati and P. Jenny, Multiscale finite-volume method for compressible multiphase flow in porous media. J Comp. Phys, 216 (2), (2006), pp. 616-636. [49] R. J . Mackinnon and G . F . Carey, Positivity-preserving, flux-limited finite-difference and finite-element methods for reactive transport. Int J. Num Mthds in Fluids, 4 1 , (2003), pp. 151-183. [50] S. Mazumder and J . V . Cole, Rigorous 3-D mathematical modeling of P E M fuel cells. J. Electrochem Soc, 150 (11) (2003), pp. A1510-A1517. [51] D . Natarajan and T . V . Nguyen, A two-dimensional, two-phase, multicomponent, transient model for the cathode of a proton exchange membrane fuel cell using conventional gas distributors. J. Electrochem Soc, 148 (12), (2001), pp. A1324-A1335. [52] H . N i , A . K . D a t t a and K . E . Torrance, Moisture transport in intensive microwave heating of biomaterials. Int. J. Heat and Mass Trnsfr, 42, (1999), pp. 1501-1512. [53] H . Ockendon and J.R. Ockendon, "Viscous flow". Cambridge University Press (1995). [54] J . Ockendon, S. Howison, A . Lacey and A . Movchan, "Applied Partial Differential Equations". Oxford University Press (1999). [55] S. Osher and J . A . Sethian, Fronts propagating with curvaturedependent speed: algorithms based on Hamilton-Jacobi formulations. J. Comp Phys, 79 (1), (1988), pp. 12-49.  126 [56] U . Pasaogullari and C . - Y . Wang, Two-phase transport and the role of micro-porous layer in polymer electrolyte fuel cells. Electrochimica Acta, 49, (2004), pp. 4359-4369. [57] D . W . Peaceman, "Fundamentals of Numerical Reservoir Simulation". Elsevier (1977). [58] K . Promislow, J . M . Stockie and B . R. Wetton, Multiphase transport and condensation in a hydrophobic P E M fuel cell cathode. P R E P R I N T (June 2003). [59] K . Pruess, C. Calore, R. Celati and Y . S. W u , A n analytical solution for heat transfer at a boiling front moving through a porous medium. Int. J. Heat and Mass Trnsfr, 30 (12), (1987), pp. 2595-2602. [60] E . A . Santillan Marcus and D . A . Tarzia, Exact solutions for drying with coupled phase-change in a porous medium with a heat flux condition on the surface. Comp and Appl Math, 22 (3), (2003), pp. 293-311. [61] J . A . Sethian, "Level Set Methods", Cambridge University Press (1996). [62] A . B . Tayler, "Mathematical Models in Applied Mechanics". Clarendon Press (1986). [63] K . E . Torrance. Phase-change heat transfer in porous media. Heat Transfer 1986, 1, (1986), pp. 181-188. [64] K . S. Udell, Heat transfer in porous media heated from above with evaporation, condensation, and capillary effects. J. Heat Trnsfr, 105, (1983), pp. 485-492. [65] V . Voller and M . Cross, A n explicit numerical method to track a moving phase change front. Int. J> Heat and Mass Trnsfr, 26 (1), (1983), pp. 147-150. [66] C . - Y . Wang, A fixed-grid numerical algorithm for two-phase flow and heat transfer in porous media. Num Heat Trnsfr, Part B , 3 1 , (1997), pp. 85-105. [67] C . - Y . Wang, Fundamental models for fuel cell engineering. Reviews, 104, (2004), pp. 4727-4766.  Chemical  127 [68] C . - Y . Wang and C . Beckermann, A two-phase mixture model of liquidgas flow and -heat transfer in capillary porous media I. Int. J. Heat and Mass Trnsfr, 36 (11), (1993), pp. 2747-2758. [69] C . - Y . Wang and P. Cheng, Multiphase flow and heat transfer in porous media. Advances in Heat Transfer, 30, (1997), pp. 93-196. [70] Z. H . Wang, C . Y . Wang and K . S. Chen, Two-phase flow and transport in the air cathode of proton exchange membrane fuel cells. J. Power Sources, 94, (2001), pp. 40-50. [71] Y . Wang and C . Y . Wang, A nonisothermal, two-phase model for polymer electrolyte fuel cells. J . Electrochem. S o c , 153 (6), (2006), pp. A1193-A1200. [72] S. Whitaker, Simultaneous heat, mass and momentum transfer in porous media: a theory of drying. Advances in Heat Transfer, 13, (1977), pp. 119-203. [73] A . W . Woods, Liquid and vapour flow in superheated rock. Ann. Fluid Mech, 31, (1999), pp.171-199.  Rev.  [74] A . W . Woods and S. D . Fitzgerald, The vaporization of a liquid front moving through a hot porous rock. J. Fluid Mech, 251 , (1993), pp. 563579. [75] D . Zwillinger, "Handbook of Differential Equations". Academic Press, (1997).  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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-0080046/manifest

Comment

Related Items