Citation
Source-and-sink method of solution of two and three dimensional Stefan problems

Material Information

Title:
Source-and-sink method of solution of two and three dimensional Stefan problems
Creator:
Li, Hongjun, 1961-
Publication Date:
Language:
English
Physical Description:
xv, 180 leaves : ill. ; 29 cm.

Subjects

Subjects / Keywords:
Boundary conditions ( jstor )
Freezing ( jstor )
Greens function ( jstor )
Heat transfer ( jstor )
Melting ( jstor )
Solidification ( jstor )
Subroutines ( jstor )
Wall temperature ( jstor )
Water temperature ( jstor )
Waxes ( jstor )
Dissertations, Academic -- Mechanical Engineering -- UF
Mechanical Engineering thesis Ph. D
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1993.
Bibliography:
Includes bibliographical references (leaves 174-179).
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Hongjun Li.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. §107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
Resource Identifier:
030431838 ( ALEPH )
31209724 ( OCLC )

Downloads

This item has the following downloads:


Full Text












SOURCE-AND-SINK METHOD OF SOLUTION OF
TWO AND THREE DIMENSIONAL STEFAN PROBLEMS










By

HONGJUN LI


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


1993



















TO MY WIFE YU LUO,
MY DAUGHTER ANGELA,
AND MY PARENTS.















ACKNOWLEDGEMENTS


I would like to express my sincere gratitude and appreciation to Dr. C. K. Hsieh, chairman of the supervisory committee, for his guidance and encouragement throughout my entire graduate study, especially for his patience in revising this manuscript. His enthusiasm and great scholarly attainments have contributed substantially to the quality of this work. I also want to extend my appreciation to Dr. Y. D. Goswami, cochairman of the supervisory committee, for his continuous support, encouragement and guidance during the research. I would like to thank Dr. W. Shyy, Dr. J. F. Klausner, and Dr. R. Mei for serving on my supervisory committee and for their fruitful comments and suggestions.

I acknowledge my colleagues and fellow graduate students for their friendship and valuable help.

Finally, my deepest gratitude goes to my wife, Yu Luo, for her unselfish love, patience, and understanding, and to my parents for their love, encouragement and support.















TABLE OF CONTENTS
Page

ACKNOWLEDGEMENTS .................................................... iii

LIST OF TABLES ...................................................... vii

LIST OF FIGURES .................................................... viii

NOMENCLATURE ......................................................... xi

ABSTRACT .......................................................... xiv

CHAPTERS

1 INTRODUCTION ................................................. 1

2 REVIEW OF LITERATURE ........................................ 6

2.1 Analytical Solution of the Stefan Problems ........... 8
2.1.1 Exact solutions .................................. 8
2.1.2 Approximate analytical solutions .............. 9
2.1.3 Other approximate solutions .................... 12
2.1.4 Analytical solutions of two dimensional
Stefan problems .. ............................. 14
2.2 Numerical Solutions of the Stefan Problems ........... 15
2.2.1 Classical formulation ........................... 15
2.2.2 Weak formulation ................................ 17
2.2.3 The boundary element method .................... 18

3 NUMERICAL ANALYSIS ......................................... 20

3.1 Problem Formulation .................................... 20
3.1.1 Governing equations ............................. 20
3.1.2 Brief review of the SSM ........................ 22
3.1.3 Time marching method ........................... 25
3.2 SSM With Superposition Techniques ..................... 26
3.2.1 The validity of the superposition method ...... 26 3.2.2 Solution method ................................. 30
3.3 Solution of One dimensional Stefan Problems .......... 35
3.3.1 Solutions in one dimensional Cartesian system . 35 3.3.2 Solutions in cylindrical coordinate system .... 39 3.3.3 Solutions in spherical coordinate system ...... 40 3.3.4 Multiple interfaces ............................. 42
3.3.5 Stability analysis .. ............................ 43










3.4 Extension to Multi-Dimensional Stefan Problems ........ 49
3.4.1 Solutions in two dimensional Cartesian system . 50 3.4.2 Solutions in cylindrical coordinate system .... 57
3.4.3 Solutions in three dimensional
Cartesian system ............................... 62
3.4.4 Stability analysis .............................. 65
3.5 Application to Latent Heat Thermal Energy
Storage Systems ...................................... 66
3.5.1 Problem formulation ............................. 67
3.5.2 Solution method ................................. 71
3.5.3 Energy conservation test ....................... 73
3.6 Application to Encapsulated PCM TES Systems .......... 75
3.6.1 Problem formulation .............................. 75
3.6.2 Solution method ................................. 79
3.6.3 Energy conservation test ....................... 81

4 TWO DIMENSIONAL PHASE-CHANGE EXPERIMENT ................... 84

4.1 Objective and Basic Considerations .................... 84
4.2 Experimental Setup and Measurement ..................... 86
4.2.1 The setup ....................................... 92
4.2.2 Test procedure .................................. 92
4.3 Sample Results and Experimental Uncertainty .......... 95

5 RESULTS AND DISCUSSION ..................................... 95

5.1 Numerical Results of One Dimensional Problems ........ 95
5.1.1 Test of accuracy and stability ................ 95
5.1.2 Subcooling effect on the interface position ... 105 5.1.3 Solution with oscillating interfaces .......... 107
5.1.4 Solution with multiple interfaces ............. 109
5.1.5 Melting of a sphere ............................ 114
5.2 Numerical Results of Two Dimensional Problems ........ 116
5.2.1 Test of accuracy and stability ................ 118
5.2.2 Long-time solutions ............................ 121
5.2.3 Comparison with experimental results .......... 124
5.2.4 Test for coalescing fronts .................... 126
5.3 Application to Two Dimensional Latent Heat
Thermal Energy Storage Systems ..................... 128
5.4 Application to Encapsulated PCM TES Systems .......... 134
5.5 Visualization of Numerical Results ................... 141

6 CONCLUSIONS AND RECOMMENDATIONS ........................... 144


APPENDICES

A DERIVATION OF EQUATIONS (3.52) AND (3.53) ................. 148










B DERIVATION OF EQUATION (3.68) ............................. 152

C FORTRAN PROGRAM FOR TWO DIMENSIONAL STEFAN PROBLEMS ...... 153 REFERENCES .......................................................... 174

BIOGRAPHICAL SKETCH ................................................. 180













LIST OF TABLES


Table Page

3.1 Expression to account for the effects of boundary conditions ..................................... 23

4.1 Thermophysical properties of paraffin wax ................ 85

5.1 Cases tested for one dimensional Stefan problems ......... 96 5.2 Properties of PCMs analyzed ................................ 97

5.3 Cases tested for two dimensional Stefan problems ......... 117 5.4 Input data for TES analysis ............................... 129

5.5 Specifications of encapsulated TES unit .................. 134















LIST OF FIGURES


Figure Page

2.1 Solution methods of Stefan problems ........................ 7

3.1 H versus x plot ............................................ 48

3.2 System analyzed in a two dimensional Cartesian coordinate system ....................................... 51

3.3 System analyzed in cylindrical coordinates: Interface moves in z-direction .......................... 58

3.4 System analyzed in cylindrical coordinates: Interface moves in r-direction .......................... 58

3.5 A schematic drawing for the thermal energy storage unit .. 68 3.6 Simplified system for analysis ............................. 68

3.7 A schematic drawing for the encapsulated latent heat TES unit .................................... 76

3.8 Control volume for fluid analysis ......................... 77

3.9 Matching conditions for PCM/HTF analysis ................. 77

4.1 Test setup .................................................. 87

4.2 Assembly drawing of the test cell .......................... 88

4.3 Exploded view of the test cell ............................. 89

4.4 Exploded view of the heater ................................ 91

4.5 Sample pictures for interface position measurements ...... 93 5.1 Temperature profile of the superposition method .......... 98

5.2 Test of accuracy and stability using different At and Ax . 99 5.3 Errors in the interface positions for different At ....... 101

5.4 Interface positions for time-variant temperature condition .................................. 102
viii










5.5 Interface positions for time-variant flux condition ...... 102 5.6 Interface positions in finite domain ..................... 104

5.7 Temperature history at the insulated wall (x=H) .......... 104

5.8 Interface positions for different subcooled conditions .... 106 5.9 Test of oscillating interfaces ............................ 108

5.10 Sketch of double interfaces generated by conditions
imposed on two boundaries .............................. 110

5.11 Temperature profiles of a double interface problem ....... 111 5.12 Interface positions of a double interface problem ........ 112

5.13 Sketch of double interfaces generated by conditions
impose on one boundary ................................. 113

5.14 Interface positions of a melting-freezing problem ........ 113 5.15 Interface positions of sphere melting .................... 115

5.16 System used in two dimensional stability test ............ 119

5.17 Interface positions at different times at x=O and x=H/2 .. 120 5.18 Interface position curves at t=3000s for different At .... 120 5.19 System used for long time solution ....................... 122

5.20 Comparison of numerical and analytical solutions for
a two dimensional problem at long time ................ 123

5.21 Comparison between numerical results and experimental
data for two dimensional phase change ................. 125
5.22 Interface position curves in a two dimensional
coalescing-front problem ............................... 127

5.23 Test of energy conservation ............................... 131

5.24 Water exit temperature at different flow rates ........... 132

5.25 Isotherms in a half PCM plate ............................. 133

5.26 Test of energy conservation ............................... 136

5.27 Water exit temperature at different flow rates ........... 137









5.28 Water exit temperature for different ball diameters ...... 138 5.29 Water exit temperature for different shell thickness ..... 139 5.30 Temperature profiles in the PCM and the shell ............ 140

5.31 Post-processor for a two dimensional TES unit ............ 142

5.32 Post-processor for an encapsulated TES unit .............. 143














NOMENCLATURE


A Amplification factor, area b Thickness of TES containing wall B,C Constants D Constant, dimension of PCM in y-direction Dh Hydraulic diameter c Specific heat d Weighting function, half width of channel E Energy f Friction factor G Green's function g Function defined in Eq.(3.110) etc. H Dimension of PCM in x-direction h Convection heat transfer coefficient k Thermal conductivity L Latent heat Nu Nusselt number n Normal direction J Bessel function Pe Peclet number Pr Prandtl number QVolume flow rate











q

R

1" R"t

Re

S Ste

T

t

U

V

W Xlylz


Greek


Thermal diffusivity Eigenvalues in the Green's function Dirac delta function Small constant Parameter in finite difference method Void fraction Function defined in Eg.(3.86) Density

Dummy variable of time Solution domain


Heat flux Radius

Coordinate in cylindrical and spherical systems Thermal resistance Reynolds number Interface position Stefan number Temperature Time

Decomposed temperature, mean velocity Decomposed temperature, volume PCM width Coordinates in Cartesian system











Subscripts

b cV f i

j 1



n P 8


xiii


and Superscripts

Boundary value Control volume

Fluid

Boundary index, spatial index

Spatial index

Liquid phase

Melting point

Time level index

Phase regions, PCM

Solid phase, interface position, sphere

Wall, shell













Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy



SOURCE-AND-SINK METHOD OF SOLUTION OF
TWO AND THREE DIMENSIONAL STEFAN PROBLEMS By

Hongjun Li

December 1993



Chairman: Dr. Chung K. Hsieh Cochairman: Dr. Yogi D. Goswami Major Department: Mechanical Engineering

A source-and-sink method has been developed for the solution of conduction heat transfer with phase change in multiple dimensions. In this method, the heat transfer in one direction is decoupled from the other directions by changing the second partial differential in these directions where the phase change is less dominant to a finite difference form and solving the problem in the other direction with an analytical solution that accounts for the motion of the interface in a phase change problem. The solution developed in this work is thus independent of the equations used to represent the interface as well as the conditions imposed on the boundaries.

The method has been applied to the tracking of single, double, and coalescing interface fronts formed by different phases assuming


xiv









equal properties. The method has been demonstrated to be accurate, convergent, and stable by numerical computation as well as experimental measurements. The method has also been applied to the solution of conjugate heat transfer including convection in a fluid medium and conduction in the phase change material, a problem applicable to the analysis of latent heat thermal energy storage systems.
















CHAPTER 1
INTRODUCTION



Phase-change problems involving melting and freezing have found numerous applications in science and engineering. Phase-change processes are encountered in manufacturing and material processing, and in annealing and continuous casting. In the important area of food processing, phase change has been an integral part of a freezing and thawing process in which one substance can be separated from another by means of a phase change. The latent heat absorbed and released by a substance during a phase change also works favorably in energy storage. Since the latent heat is usually large for some substances, a small volume-to-energy ratio can be achieved. This makes the process particularly attractive for solar energy applications. Recently, with the increased use of lasers in industrial processes such as cutting, drilling, surface hardening and densification, the phase change has added another dimension to contemporary research.

Phase-change problems have been commonly referred to as a class of moving boundary problems. In contrast to many conventional physical problems in which the solutions are sought in a fixed domain, the phase-change problems are solved in a variable domain.











There are two phases in the medium which are separated by a phasechange interface. Under the continuous heat transfer with the surroundings through the boundary, this interface position can change with time. The interface position is unknown a priori and must be determined by solving nonlinear equations. The phase-change problems thus pose a real challenge for analysis.

Phase-change problems that have been solved in the literature have often been treated with various degrees of simplification. Earliest phase-change analyses have been confined strictly to a heat diffusion study in a domain imposed with boundary conditions that are invariant with position and time. Later models account for the property variations in different phase regions and for a medium imposed with time-variant conditions. While these are more realistic physically, they also add complexity to the analysis to the extent that these problems can no longer be solved exactly. Particularly noteworthy is the recent development of methods for the solution of phase change coupled with convection and radiation. These problems can only be solved numerically. These numerical methods have recently been further extended for the solution of phase change due to a moving heat front, a problem of great importance in welding and annealing. A detailed literature review is provided in Chapter 2.

Phase-change heat transfer problems can be classified according to the configuration and the formation of the interface in the medium. Thus, in some analyses, the phase change is assumed to occur at a distinct temperature which does not change with time. In such cases,











the solid and liquid regions are separated by a sharp front whose position can be determined rigorously. This permits the interface position to be verified experimentally. In some materials, however, the phase change occurs over a narrow temperature range. Then, the phase change takes place in a mushy zone where the phase change is smoothly dispersed. Probably most intriguing of all is the phase change in alloys. In the alloys, the phase-change temperature is not distinct. There is a crystalline-like structure at the interface which is not necessarily smooth or even continuous. In practice, materials with crystalline structure and mushy zone can be analyzed more conveniently with computational approximations. This is not so, however, for a material with a distinct phase-change temperature. The interface position in such a material must be determined more rigorously. This interface position can also serve as a basis for evaluation of the accuracy of the analysis.

Phase-change problems that have been solved in the literature are mostly for domains in one dimension. Solutions of two dimensional problems have rarely been attempted. Although an enthalpy method can be used to solve such problems conveniently, that method only works well if the phase change occurs in a mushy zone.


It is thus the purpose of this work to develop a hybrid scheme that employs a finite-difference method in one direction and a sourceand-sink method in another direction. The former will be used in the direction where the phase change is less dominant, while the latter will be used in the other direction where the phase-change interface









4

moves rapidly. The source-and-sink method is selected primarily for its convenience and accuracy as demonstrated in eight previous studies (see Hsieh, 1993; Hsieh and Choi, 1992a, 1992b; Choi and Hsieh, 1992; Choi, 1991). As will be shown, the source-and-sink method only requires one set of equations to solve for temperatures in 'all' phase regions. By contrast, it is necessary to use a set of equations for Each' phase region in the traditional method of solution by the Neumann-Stefan formulation. The work involved in the solution can thus be reduced considerably for the source-and-sink method. In addition, in the traditional method of solution, the interface Stefan condition must be given as a separate condition, whereas in the source-and-sink method, this condition has been incorporated into the governing heat diffusion equation. The source-and-sink method is thus highly efficient in the solution of the phase-change problems.

In the present work, the source-and-sink method is subjected to a rigorous test for its feasibility for superposition. The method is then used to develop strategies for the solution of the problems in Cartesian, cylindrical, and spherical coordinates. For problems in Cartesian coordinates, the method has been tested analytically for its accuracy for long-time solutions. An apparatus has also been built for experimentally testing its accuracy for short-time solutions. In both efforts, the interface positions are used for comparison. Finally, the method has been applied for the analysis of two latent heat, thermal energy storage systems.








5

This dissertation is divided into six chapters. A brief literature review is provided in Chapter 2. Analytical methods that are related to the solution of phase change in multiple dimensions are compiled in Chapter 3. An experimental setup is described in Chapter 4. The analysis developed in this dissertation is used in a wide range of case studies and their results are presented in Chapter 5. Finally, conclusions and recommendations are given in Chapter 6 to close this dissertation. To facilitate use of the analysis developed in this work, computer programs written in FORTRAN language are compiled in the appendices. It is hoped that the insight gained from this study will contribute to the solution of phase change in multiple dimensions in the future.















CHAPTER 2
REVIEW OF LITERATURE



Study of phase-change problems dates back to 1831 when Lame and Clapeyron (see Rubinstein, 1971) first studied the thickness of the solid crust formed by cooling a liquid under a constant temperature condition. Ever since, the problems have drawn attention from various quarters, and today the phase-change principles have been greatly extended and applied to study various scientific phenomena. The phasechange problems have been in studied for over one hundred and sixty years, and the phase-change literature is spread over several disciplines. It is thus a daunting task even for a brief review in this chapter. Since the present work strictly deals with the development of a two dimensional, diffusion-driven, phase-change solution scheme that combines an analytical solution in one direction and a numerical solution in another direction, only mathematical methods for diffusion-driven phase change are relevant here. It is thus the purpose of this chapter to provide a brief review of the mathematical techniques popular for the solution of phase-change problems by mechanical engineers; see a summary of methods diagramed in Figure 2.1. Interested readers are referred to books and monographs (e.g., Crank, 1984; Yao and Prusa, 1989) for a detailed treatise on the subject.





















Method of solution of Phase-Change Problems


Analytical Solutions








Numerical Solutions


Exact Solutions (Neumann Solution)



Approximate
-Analytical Solutions


- Power Series Solution Methods
- Solution of Integrodifferential Equations
(Source-and-Sink Methods,
Embedding Methods, etc.)
- Solution of Integral Equations
- Others (Perturbation, Transforms, etc.)


- Classical Formulation (Solution in Temperature Domain) - Finite Difference Methods

- Finite Element Methods Weak Formulation T Boundary Element Methods
- (Solution in Enthalpy Domain)


Figure 2.1 Solution methods of Stefan problems









8

2.1 Analytical Solutions of the Stefan Problems



2.1.1 Exact solutions

In the heat transfer literature, phase-change problems have often been referred to as Stefan problems in recognition of Stefan who first reported the solution of such problems in the open literature in 1889 (Carslaw and Jaeger, 1959). However, it has now been fully settled that the first exact solution of the problems was attributed to Neumann some thirty years earlier. Neumann's work was concealed in a series of lecture notes but was complete in the sense that it provided derivation not only of the square root time law that governs the motion of the phase-change interface as predicted earlier by Lame and Clapeyron but also of the value of the proportionality constant which was missing in Lame and Clapeyron's work. Thus, for the first time in history, Neumann solved the Stefan problems completely.

Neumann's work is important for the fact that, despite innumerable efforts made by others in search of additional exact solutions after Neumann, Stefan problems that can be solved exactly today remain those that were originally solved by Neumann some one hundred and thirty years ago. Indeed, it has now been firmly established that, for a Stefan problem that can be solved exactly, the problem must be amenable to a similarity transformation. This limits the problems to be solved in an infinite domain, composed of a material of constant phase properties, and imposed with a constant temperature boundary condition.








9

2.1.2 Approximate analytical solutions

With the great limitation imposed by the exact solution, it is not surprising to see a wide range of approximate solutions developed and documented in the literature. One of the early approximate methods expands the interface position and the temperature in the medium in terms of a power series or error and complementary error functions. These series are then substituted into the governing equations and boundary conditions to determine the coefficients in the series expansion. Tao used this method in a series of studies. He also extended the method for study of phase change imposed with arbitrary initial and boundary conditions (Tao, 1978, 1980, 1981).

Power series methods find great usefulness in the determination of short-time solutions. At short time, the series converge rapidly; only a few terms are thus sufficient for an accurate solution. However, at long time, more terms are necessary. The methods become tedious and the determination of the coefficients becomes time consuming. They lose their appeal.

Phase-change problems can also be solved by using integrodifferential equations. Lightfoot (1929) developed a moving heat source method for the solution of solidification imposed with a constant temperature boundary condition. In his method, liberation of latent heat due to solidification is taken as a moving heat source. The temperature in the medium is expressed in terms of a Green's function. Since Lightfoot used this method for solving a problem that is possible for an exact solution, the integrodifferential equations








10

obtained in his work can be solved by a similarity transformation. This leads to results expressed in error and complementary error functions (Crank, 1984). Thus, Lightfoot's integrodifferential solution effectively retrieves Neumann's exact solution.

Using a moving heat source front for solidification naturally leads to the use of a moving heat sink for melting. Lightfoot's method can thus be taken as a source-and-sink method. Chuang and Szekely (1971) have used this method in the development of a numerical scheme for the solution of Stefan problems. Recently, the source-and-sink method has been greatly extended for the solution of Stefan problems imposed with temperature and flux conditions that are either constant or variable with time (Hsieh and Choi, 1992a,1992b; Choi and Hsieh, 1992). In these efforts, the method has been shown to be highly accurate particularly for solution at long time. It has thus been used to assess the errors in some other solution methods frequently used in the literature. The method has also been applied to the solution of phase change subjected to cyclic temperature and flux conditions (Hsieh and Choi, 1992a). In these problems, both melting and freezing fronts appear simultaneously in the medium. Recently, the source-andsink method has been shown to be equivalent to the boundary element method (Hsieh, 1993). This discovery carries promise of broadening the utility of both methods for the solution of heat transfer problems.

An important feature of the source-and-sink method is that the method uses only one set of equations to solve temperatures in 'all' phase regions. By contrast, in the conventional methods, one set of








11

equations is required to solve temperatures in 'each' phase region. Furthermore, in the source-and-sink method, the interface flux condition has been incorporated into the governing equation. The number of equations solved in this method is thus far fewer than with the conventional methods. However, the method works best for a medium of constant and equal properties in different phase regions. For a medium with different property values, double sources and sinks or a combination of the source-and-sink and the boundary element can be used as suggested by Kolodner (1956) and Hsieh (1993). However, the equations are lengthened considerably. The method loses its attractiveness.

Integrodifferential equations have also been used by Boley (1974) in the solution of Stefan problems. Following a different approach, Boley embedded the time-variant phase regions into a large domain enclosed by a fixed boundary. The problem can then be solved in a fixed domain. However, since the conditions imposed on the fixed boundary are intended to reproduce the conditions existing on the real physical boundaries including the interface, the conditions on the fixed boundary become unknown. In one case study, Boley (1974) introduced a fictitious flux condition at the fixed boundary. This results in the solution of a system of integrodifferential equations for the unknown flux and the interface position.


Stefan problems can also be solved by a heat balance integral method (Goodman, 1958). Essentially a branch of a large class that is commonly known as the method of weighted residuals, the heat integral








12

method is convenient to use particularly in the solution of phase change in a medium that is initially at the phase-change temperature (i.e., one-phase problem) and imposed with a time-invariant boundary condition. In this method, the governing equation is changed to an integral form and a trial function (usually polynomials) is used to represent the temperature of the medium. This temperature function is then substituted into the integral equation to solve for the relation of the interface position with time. Since the temperature function must be assumed a priori, the accuracy of the solution is related to the appropriateness of this function in the solution. Goodman and Shea (1960) and Yuen (1981) have also applied this method to solve phasechange probems with two phase regions.

It has been widely accepted that, in the heat balance integral method, using a higher-order polynomial for temperature in the trial function provides no assurance for improvement of accuracy in the solution (Goodman, 1964; Langford, 1973). However, a spatial subdivision has been found to be useful as reported by Noble (1975) and Bell (1978). Bell (1979) has further applied this method to the analysis of solidification around a cylindrical pipe by a Runge-Kutta numerical integration algorithm.



2.1.3 Other approximate solutions

Stefan problems can also be solved by using quasi-steady and quasi-stationary approximations. In the quasi-steady approximation, the unsteady terms in the governing equations are dropped. This








13

results in a temperature that is good for long-time solution. As for the quasi-stationary approximation, the interface velocity term in the governing equations and the interface condition are all discarded. As such, the interface is rendered stationary. Since the initial condition must still be satisfied, the temperature is good for short-time solutions.

In the asymptotic method of solution, these quasi-steady and quasi-stationary limits are used for inner and outer expansions in the construction of temperature and interface position. Duda and Vrentas (1969) have demonstrated that the quasi-stationary approximation is actually the zeroth order solution of a regular perturbation, which is, indeed, good for a short-time solution. On the other hand, Weinbaum and Jiji (1977) have shown that the lowest order term of a long-time perturbation actually leads to a quasi-steady limit and is thus good for a long-time solution. These perturbation methods work well for small Stefan numbers. Since the mathematics involved in the asymptotic solution is prohibitively tedious, the method is practically useful only for first-order approximations. Solomon (1978) developed an approximate solution for a phase change imposed with a constant heat flux condition. His method follows that of the integral method; however, the energy balance is not preserved in the solution.

Phase-change problems have also been solved by a Fourier transform (Ockendon, 1974). In this method, the Stefan problems are reduced to an initial-value problem. In a similar fashion, a Laplace transform can be used to change the problems into boundary-value








14

problems. While the transforms are relatively easy to apply, the inversion poses a problem. The inversion process often leads to the solution of integrodifferential equations.



2.1.4 Analytical solutions of two dimensional Stefan problems

Stefan problems that have been solved in the literature are mostly one dimensional. There has been a lack of solutions for two dimensional problems. Poots (1962) extended the heat balance integral method for the solution of phase change in a square prism imposed with a uniform temperature condition. He was able to show that his analytical results were in reasonable agreement with the numerical results obtained earlier by Allen and Severn (1962). Two dimensional Stefan problems can also be solved by using an embedding method as shown by Sikarski and Boley (1965). Following a totally different approach, Rathjen and Jiji (1971) were able to combine Lightfoot's moving heat source method and Patel's treatment of the interface condition (Patel, 1968). They developed an analytical solution for solidification in a two dimensional corner. Budhia and Kreith (1973) followed essentially the same approach in the analysis of phase change in a planar wedge. In these efforts, the interface position must be assumed a function a priori, and the constants in this function are determined by substituting the function into integrodifferential equations. The success of the method thus hinges on the use of the proper function in the early stage of solution.













2.2 Numerical Solutions of the Stefan Problems

It is expected that, in the numerical solution, the phasechange problems can be solved by using a finite difference method (FDM), a finite element method (FEM), and a boundary element method (BEM). It is noted that, in the study of heat diffusion without phase change, the problems are usually analyzed in a temperature domain. This is not so, however, for diffusion with phase change as in a Stefan problem. In the Stefan problem, because the phase-change interface is unknown a priori and must be determined as a part of the solution, there is an alternative to the solution in the temperature domain. In fact, the problem can be solved more conveniently in an enthalpy domain. The former (solution in temperature domain) requires a classical formulation, while the latter (solution in enthalpy domain) requires a weak formulation. They are briefly reviewed in the sections that follow.



2.2.1 Classical formulation

In the classical formulation, temperature is targeted for solution. The governing equations, the initial and boundary conditions are all expressed in terms of temperature and they are used together in the solution of the interface position as well as the temperature in the medium. This method works well if the interface position is the focus of investigation. As such, the method is useful if the material has a distinct phase-change temperature.











In the solution of the phase change by a classical formulation, it is difficult to track the interface position. Since the predicted interface position may not fall on the grid lines initially set up in the solution, iteration is usually necessary which leads to complexity in the solution. A remedy to this is to use adaptive grids or transformation of variables. In these methods, the interface position is always tied to a grid line or is immobilized in the transformed domain. Furzeland (1980) has made an extensive study for comparison of numerical methods. Some popular methods have been provided in works by Wilson et al. (1978) and Crank (1981).

In the numerical solution by the FEM and FDM, the computational grids can be either fixed or time-variant. In the fixed grid method, an interpolation technique is usually necessary at each time step so that the interface positions can be determined accurately (Murray and Landis, 1959; Hastaoglu, 1986; Hastaoglu and Baah, 1991). By contrast, in the variable grid method, the grid must be modified at each time step (Gupta and Kumar, 1980,1981). Bonnerot and Jamet (1981) solved a one dimensional Stefan problem by using a third-order-accurate finite element method with an adaptive grid. In their work, the computational grid is re-constructed at the end of each time step and the results of the previous time step are mapped onto the new grid by interpolation. These methods are often referred to as front-tracking methods and have been used extensively in the numerical solution of phase-change problems.








17

An alternative to the fixed grid method is to introduce an appropriate set of space coordinates so that the governing equation and boundary conditions can be transformed from the physical domain to a computational domain. In the computational domain, the interface position rests on a grid line. This method fits into the category of front-fixing or boundary-fixing. Furzeland (1980) and Lacroix and Voller (1990) have made extensive studies for comparison of the fixedgrid method and the transformed-grid method. As expected, when the grid transformation is extended to the solution of phase change in two and three dimensions, the choice of the transformation coordinates, known as body-fitted coordinates, becomes particularly important. It is thus necessary that the interface front curve be smooth and continuous and the number of grid lines between the interface and the boundaries be fixed. Hsu et al. (1981) combined the transformation method with an energy-conserving finite-difference scheme in the analysis of freezing outside a cooled tube. Shyy et al. (1992) used body-fitted grids together with an adaptive grid technique in the solution of solidification in continuous casting. Additional numerical solutions by FDM and FEM have been documented in Oberkampf (1976); Saitoh (1978); Pham (1985) and Liang (1993), among others



2.2.2 Weak formulation

In the weak formulation, the problems are solved in an enthalpy domain. Analogous to the Rankine-Hugoniot equation in the shock-wave analysis, the weak formulation can be used to solve phase-change problems without searching for the exact position of the interface. It








18

is thus superior if the interface position is not targeted for investigation (Elliot and Ockendon 1982; Boucheron and Smith 1985).

In the enthalpy method, the phase-change interface is determined a posteriori from the numerical solution carried out in a fixed domain. A survey of the recent developments of the enthalpy method has been given in Tamma and Namburu (1990) and a detailed formulation of the method can be found in Alexiades and Solomon (1993).

Although the enthalpy methods have been used in the solution of phase change in many engineering problems, it is found that those methods work well if the phase change occurs over a mushy zone. For problems with a distinct phase-change temperature, the nonlinearity behavior of the enthalpy near the phase-change temperature presents great numerical difficulties which necessitate the application of an arbitrary smoothing function of the enthalpy (Goodrich, 1978) and suppression of the oscillation in the solution (Voller and Cross, 1983; Pham, 1985, 1986). While these difficulties can be resolved by tracking the interface, it also complicates the solution. Most of all, it defeats the purpose of using the enthalpy method which is intended for bypassing the search for the interface in the first place (Yao and Prusa, 1989).



2.2.3 The boundary element method

Phase-change problems have recently been solved by a boundary element method (BEM) and a combination of BEM and FDM. In the BEM, with the use of a time-dependent fundamental solution, only boundary








19

grids are required in the solution. O'Neill (1983) developed a noniterative boundary element analysis, which permitted the temperature gradients normal to the elements at the phase-change interface to be determined directly in the solution. Zabaras and Mukherjee (1987) developed an iterative implicit algorithm using BEM, which was shown to be free from major time-step limitations. Sadegh et al. (1987) applied BEM to the solution of solidification around a buried pipe in a semi-infinite medium with initial superheating. The planar boundary was thus initiated at the surface and propagated inward with time. The BEM has also been applied to the solution of solidification of metal casts in metal or sand molds (Hong et al. 1984). This problem has also been solved by coupling the boundary element method and the finite difference method in a conjugate heat transfer analysis (Hong et al. 1984). More recently, two and even three dimensional phase-change problems have been analyzed for laser heating, drilling, and cutting (see, for example, Basu and Date, 1990; Sarler et al. 1991; Hsieh et al. 1992). Hsieh (1993) showed that the BEM is closely related to the source-and-sink method and their unification is discussed in a great detail.















CHAPTER 3
NUMERICAL ANALYSIS



In this chapter, the solution of Stefan problems by a sourceand-sink method (SSM) is presented. To ensure both accuracy and efficiency, a superposition technique with a new time marching scheme is introduced. The validity of the superposition technique for the solution of nonlinear Stefan problems is verified. Finally, the SSM with superposition technique is extended to the solution of multidimensional Stefan problems.



3.1 Problem Formulation

3.1.1 Governing equations

Consider a phase-change process that is diffusion driven in a phase change material (PCM) having constant and equal properties in different phase regions. For such a process, the governing equation can be written as




1liTP(F, t)
a t -V2Tp(F,t), p=s,l, t>to, FEQ (3.1)



where the subscripts I and s denote the liquid-phase region and solidphase region, respectively.








21

The initial and boundary conditions are given as follows:



T(F,to) = To(F) (3.2) OT(Fi, t)
inbi +hi[T(9i,t)-T,]=q, i=1,2,... (3.3)



where nbi denotes the outward-drawn normal of the ith boundary. For the sake of generality, the initial condition, Eq.(3.2), starts at to. The boundary condition, Eq.(3.3), is also general in that it reduces to the boundary condition of the first kind when ki and qi are taken to be zero. It reduces to the boundary condition of the second kind when hi is taken to be zero.

In addition to the initial and boundary conditions, the interface between the liquid and solid regions is subjected to the temperature and flux conditions. Known as interface conditions, they are given by the following relations:



T, (S,t) = Ts (S, t) = Tn (3.4)



OTs(S,t) OT(S,t)_piLl dS(t) (35) 0n an k dt


where n is in the normal direction of the interface. L is the latent heat of fusion (L>O for freezing, L









3.1.2 Brief review of the SSM

To describe the source-and-sink method, a one dimensional Stefan problem is considered. In this method, the interface condition, Eq.(3.5), is incorporated into the governing equation to form an equivalent problem (Hsieh and Choi, 1992a,1992b), which, in one dimensional Cartesian system, takes the following form



1lOT= 2 T Ld6
Ot o2 + Ld- 6(x-S(t)), t>t, xEQ (3.6)



The corresponding initial , boundary, and interface conditions are given as follows:


T(x,to) = To(x) (3.7) OT ( x,t)
k� hObit +hi[T(xi,t)-T ] = qi (3.8)


T(S,t) = Tm (3.9) The general solution of this problem can be expressed in terms of a Green's function as



T(x, t) f To(x')G(x, t I x', to)dr'


t
+aJ [BCi]dr
rt


+ S m.(T) I S(r), 7)d-r (3.10) rt
0-t








23

In Eq.(3.10), the first term on the right hand side accounts for the contribution of the initial condition, the second term accounts for the boundary conditions, and the third term accounts for the motion of the interface. The expressions for [BC] in the second term vary with the specific condition imposed at the boundary and are listed in Table

3.1.



Table 3.1 Expressions to account for the effects of boundary conditions

Boundary conditions Expressions for [BC] in Eq.(3.10) i OG(x, t Ixb,7)
T(zb,t) = Tb(t) - T"b nb OT(Xb,t) qb(t) q)(_ , t I ,
n b kb
aT(xb,t) hb~r
bt)+ T(xbt) = -bHb(t) H G(x, t I xb, r)
anb b kbkb


Hsieh and Choi (1992a, 1992b) and Choi and Hsieh (1992) used the SSM to solve one dimensional Stefan problems in a semi-infinite domain. In the problems they solved, the medium was initially at a uniform temperature. The boundaries were imposed with constant or time-variant temperature and flux conditions. To account for possible initial subcooling or superheating of the medium, the solutions were divided into two stages: a pre-melt stage and a phase-change stage. The solutions for these stages were then combined into a single equation as











t
T(z,t) -To(z,t) �-H(t-to)L dS(-r) (, t IS(r),r)dr (3.11)



where to is the time when phase change took place, and H(t-to) is the Heaviside function given by



0 for t< to
H(t-to) = {
1for t >_to


Since the initial temperature was uniform, the first term on the right hand side of (3.11) could be expressed in a compact form as t 2
To(X, t) - [ E(r) e4.(tr)dr (3.12) r=tot_ l2



where

xTb(r) for the first-kind boundary condition
E~){2a t-r


;qb(T) for the second-kind boundary condition



In practice, since the interface position is unknown a priori, Eq. (3.11) must be solved numerically. In this effort, the entire time range is divided into small time intervals, At, over which numerical integration is performed as

t N tn f F(r)dr = F(r)d (3.13) 1-=to n= t-










where tN =t and t"=to+(n-1)At. One of the advantages of this numerical method is that its accuracy does not degrade with time. Large errors only occur over the first few time steps and they diminish rapidly with time. Accurate results can thus be obtained with this method. However, since the time integration always traces back to the initial condition at to, the effort of computation may be considerable, particularly at long time.



3.1.3 Time marching method

Another approach to the solution of Eq.(3.11) is the use of a time marching scheme. In this approach, the solution is expressed as



T(,x, tn ) = f T(z', tn)G(z, tn+1 I a', t")dz'


tnT _n_2

t+l_2
+'g E(T) e4.(tn+l -') dr
,=t ( tn+l-,r)l/2

tn+l
+ k S(r) , tn+1 S(7),r)dr (3.14) 7=tn



The time integration is thus evaluated one step at a time, and the temperature evaluated at each time step serves as the initial condition for the succeeding time step and it moves on with time. For this method, a domain integral emerges as the first term on the right.

In this term, the temperature at time tn is T(z,tn) which is nonuniform. A numerical integration must therefore be performed, giving rise to additional computational effort. It also degrades the accuracy, which








26

is particularly pronounced when the domain for integration extends to infinity.



3.2 SSM With Superposition Techniques

To alleviate the problems stated in the preceding section, a superposition method with a new time marching scheme is developed. Before addressing this method, the validity of the superposition method for the solution of nonlinear Stefan problems is critically analyzed.



3.2.1 The validity of the superposition method

The superposition method has been widely used in the solution of linear boundary value problems. Although this method has also been applied to the solution of nonlinear Stefan problems, its validity has never been established. For the ease of derivation in what follows, it is assumed that the domain is one dimensional. The PCM has constant and equal property values in different phase regions. Thus, the governing equations are all linear. The initial and boundary conditions are also linear if the radiation effect is incorporated into the convection effect by means of an overall heat transfer coefficient. The only nonlinear relation comes from the interface flux condition. It is thus necessary to decompose this condition if the superposition method is to be applied to the solution of the Stefan problem.

To start the derivation, a small constant c is introduced such that, for a melting problem,












T(S+c,t)=Ts(S,t)


As c approaches zero, lim T(S�ct) = Tm By definition, the differential of T at x equal to S can be recast as


dTI _OT(S�,t) �
- dT + -_ 0


(3.16)


where, on the right hand side, the partial differentials are to be evaluated in the same phase region. That is to say, a S+c in OT/ax is to be followed by a S+c in aT/at. Then dS/dt can be expressed as


dS _ aT(S � c, t)/at 7T -aT(S4c,t)/cxLO


(3.17)


Substituting (3.17) into (3.5) changes the interface condition to the following


[ Ts aTl, pILl aT(S�c,t)/at
Lax ox]=s- k aT(s�, t)I o


(3.18)


which is clearly nonlinear.

An alternative form of Eq.(3.17) can be derived by writing it in both solid and liquid regions and equating their results as


dS OT(S+c, t)/at OT(S-c, t)/at 1o
idt OT(S+c, t)/ax o= OT(S-, t)/ax (3.19)


T(S-c,t)=Tt(S't),


(3.15)










which can be recast as


T(S-,t) o O t [OT(S+ct)/tx11 - (3.20) [07'(S+c, t)/OtJJ


Since for the Stefan problem


[aT(S+,t) OT(S-,t) 0 (3.21) SLx ax (+o


because of the flux jump at the interface, one can multiply the numerator and denominator of Eq.(3.17) by the partial differential difference appearing on the left hand side of (3.21) as


dS _ [aT(S+,t)/Ot aT(S+,t)/Ox-OT(S- ,t)/Ox 1
dt -[T(S+ct)---i OT(S+,t)/Ox-OT(S-,t)/ox - (3.22) Then, expanding the product in the numerator of the fractions in this equation and introducing (3.20) and canceling out terms enable the following relation to be derived



dS _ aT(S+c,t)/t-aT(S-,t)/Ot (3.23) dt OT(S+,t)/Ox-T(S-c,t)/x (3LO Equation (3.18) can then be expressed as


[ O T11 _ _pLl[OT(S+ct)/t-OT(S-c,t)/Otl
5x Ox]=s - k [aT(S+,t)/Ox-OT(S-ct)LxJo (3.24)


These equations will now be used to validate the use of the superposition principle in the solution. In this effort, the








29

temperature is decomposed into a linear part U and a nonlinear part V such that T(x,t)=U(x,t)+V(x,t). For the diffusion problem solved, U is a smooth function of x and t. It thus has continuous first derivatives even at the interface, given by the relation



9U(S+,,t) U(S-Ct) _U(S+Ct) OU(S-,t)0 as c-40 (3.25)
Ox Ox - Ot Ot


Equations (3.23) and (3.24) can then be reduced to



dS _ OV(S+E,t/Ot-OV(S-f,t)/Ot (3.26) dt- OV(S+,,t)/Ox-OV(S-,,t)/Ox (O




[OI' aV pILI dS(t) (3.27) [O-x Ox __-XX - k dt


Equations (3.26) and (3.27) show clearly that U can be totally removed from the interface flux condition; only the V problem has a nonlinear condition at the interface. The key point in the superposition method is that the U-problem is linear and smooth at the interface, while the V problem is nonlinear which contains all the features characteristic of the Stefan problem.

It is noted that the derivation above can be easily extended to the solution of multi-dimensional Stefan problems. The derivation also illustrates well that the superposition principle is only valid for situations when the PCM has constant and equal properties in different regions.










3.2.2 Solution method

To describe the solution procedure, a general one dimensional Stefan problem is taken for analysis. For this problem, the governing equation is given by the relation


1 1xdaT) at --f _d O-X)


t > to, z E


where xd is the Sturm-Liouville weighting function, given as follows:


0 lab
d- cylinder
- 2 sphere


The initial, boundary, and are, respectively,


interface temperature and flux conditions


T(x, to) = To(z)


(3.29)


kiaT (xi) +hi[T(xi, t)-T.] =qi
I anbi


T(S, t) = Tm


(3.30)


(3.31)


[aT8 aT pill dS(t) ax -a -L k dt


(3.32)


Application of the superposition as described in the preceding section leads to the solution of two problems in which the linear problem U(x,t) satisfies the following equations


(3.28)












au 1


t > to, x E


U(x, to) = To(x) (3.34)


k, + h.(Uj - T.) = qj (3.35)



In this problem, Eq.(3.25) is satisfied implicitly because of the smoothness of U.

Correspondingly, the V(x,t) problem satisfies the following relations


10V 1
O t -td a t, 5-x '


t> to, XE Q


V(X, to) = 0


[a Ov J [ a Ox -5I1=


_pill dS(x)
k dt


(3.40)


The U-problem, being linear, can be solved exactly by established methods, such as separation of variables or integral transformations. However, these analytical solutions often lead to


(3.33)


(3.36)


ki- ' + hiVi = 0


(3.37)


(3.38) (3.39)


V(S, t) = V,(S, t) = Tr,.- U(S, t)








32

infinite series of trigonometric or Bessel functions, which are slow to converge. In addition, the eigenvalues must often be found numerically. For these reasons, the analytical solutions may be impractical. Use will therefore be made of a numerical solution by a finite difference method (FDM) for the solution of the U-problem.

To solve the V-problem, the source-and-sink method is employed. Following Hsieh and Choi (1992a), the interface condition (3.40) is incorporated into the governing equation (3.36) to form an equivalent problem as follows



1OV=-LaxdOV+ LdS (-S(t)), t>to, xEQ (3.41)



This equation can be solved in terms of the Green's function. In this effort, since the initial and boundary conditions are all homogeneous, (see Eqs.(3.37) and (3.38)), the V problem can be solved as

t

V(X,t) =J d (r7)sd(r)G(XtlS(),r)dr (3.42)
0

Application of the interface condition, Eq.(3.39), gives a relation




rt


which is to be solved implicitly for the interface position. In this effort, it is noted that, if S is located at a grid point, U(S,t) can be directly obtained from the finite difference solution. Otherwise, a cubic spline will be used to determine U at point S.












The solution procedure can now be stated as follows. For each time step after phase change, the T-problem is decomposed into a Uproblem and a V-problem. The U-problem satisfies the nonhomogeneous initial and boundary conditions (see Eqs.(3.33) through (3.35)), while the V-problem satisfies homogeneous initial and boundary conditions together with interface conditions (see Eqs.(3.36) through (3.40)). In the solution, it is convenient to use (3.41) as the governing equation for the V-problem, which incorporates (3.36) and (3.40), a reduction of one equation. Then the U-problem is solved first by FDM. The Vproblem is solved next by SSM. These solutions are then superimposed and this is done for each time step. In this method, the temperature evaluated at the new time level will be used as the initial condition for the U-problem at the next time level and it moves on incrementally with time until the desired time is reached. A summary of the procedure can be given as follows:



(1) Start from an initial condition To(x) and set n=l,

(2) Set to=t,t=t+ and solve Eq.(3.33) for U(x,tn+l),

(3) Solve Eq.(3.43) for S(tn+l),

(4) Calculate V(x,tn"l) from Eq.(3.42),

(5) Find T(x,t0+l) = U(x, tn+l)+v(X, tn+l),

(6) Set To(x)=T(x,tn+l), n=n+l, and repeat (2) to (6).



An important feature in this solution method is that the integral in V(z,t) (see Eq.(3.42)) is evaluated one step at a time and








34

the V-problem always starts with homogeneous initial and boundary conditions. No domain integration is thus involved. Since a time marching scheme is employed, the method is equally efficient and accurate for short time as well as for long time.

In the numerical evaluation of the integrals in Eqs.(3.42) and (3.43), a multi-point Gauss-Quadrature formula is used. With sufficiently small time steps, the interface velocity can be approximated as


dS(r) S(tn~l) - S(tn) ( .4
dr - At


where tn <,r < tn+1, and At = tn+- tn. Then, only one unknown S(tn+l) appears in Eq.(3.43). This unknown can be found by using a rootfinding method. To ensure a fast and convergent solution of this method, a conjugate Newton-Raphson and bisection method is employed. It is expected that, for problems imposed with a sudden change of temperature at the boundary, the interface velocity dS tends to be infinity at time equal to to. A large error may thus result if a finite and constant velocity is assumed in the first time step (see also Hsieh and Choi, 1992a). To reduce this error, a slow-start procedure is used in which the time step is treated as a variable. It starts with a very small value and increases gradually with time until it reaches a preset upper bound. In this way, the errors can be reduced considerably.








35

3.3 Solutions of One dimensional Stefan Problems

The analysis given in the preceding section will now be used to solve one dimensional Stefan-problem examples in Cartesian, cylindrical, and spherical systems. For the examples in Cartesian system, both semi-infinite and finite domain (plan slab) will be solved, whereas in the cylindrical and spherical systems, only finite domain of the cylinders and spheres will be considered. For these examples, the governing equation, initial, boundary, and interface conditions have been given as Eqs.(3.28) through (3.32). The U and V problems are formulated as Eqs.(3.33) through (3.40). The V problem has been solved using the equivalent problem (3.41) with the result given as (3.42). The interface position is to be solved using (3.43). Solutions for these examples are now compiled in the subsections that follow.



3.3.1 Solutions in one dimensional Cartesian system

U-problem:

The governing equation for the U-problem in a Cartesian system can be obtained from Eq.(3.33) by setting d equal to zero. To ensure numerical stability for all time steps, an implicit finite difference method is used with the results given as follows:



-- 1 U+l l+20x)u' +l' n~ (3.45a) or

. 1 + Ol _re (3.45b)
h xj-1 + ( -+=x)t/Ax2. - He2rXej+, th _seXU-i n( refers t X me+a


where 0,= aAt/Ax. Here, the superscript n refers to the time level and








36

the subscript j refers to the nodal point location. Thus,

Uj! = U(X, tn), ujl r+l1 =VU(X� A, 0+l)



Notice that Eq. (3.45a) is second-order accurate in space and firstorder accurate in time and is often referred to as the simple implicit method. Equation (3.45b), however, is second-order accurate in both space and time and is commonly known as the Crank-Nicolson method. Both methods have been well established and are unconditionally stable for all time steps. These equations can be solved by using a tridiagonal system of linear algebraic equations with a simple forward- and backward-substitution, an efficient numerical procedure.

V-problem:

The solution of the V-problem can be directly obtained from Eq.(3.42) by setting d to zero. For simplicity sake, let t0 be to and tn+1 be t1. Then, at the end of the first time step, t1





For the PCM in a semi-infinite domain, the Green's function is given by


( X_ )2 _ (X+X )2
G(x, t I ',r) -e 4(tl ) e 4e (t1-)} (3.47)



where the plus (+) and minus (-) signs inside the brackets are to be used for flux (second-kind) and temperature (first-kind) conditions at the boundary, respectively.








37

A close examination of the Green's function in (3.47) reveals that a singularity is located at r=t1. This singularity may degrade the accuracy of the numerical integration and can be resolved by introducing a transformation


?I= 4c-(tj - r) (3.48)


so that Eq.(3.46) is changed to




'4oaAt -(x.-S)2 -(x+S)2 tt J 2 1 e 72 e 2 ] dq (3.49) 17=0





For the PCM in a finite domain, such as a plan slab (O
00

G(x,t x',nr) = Hjfr ea~n(t)sin(fnx)sin(fnx') (3.50) However,


n-, for a temperature condition at x =H,

fn = (2n-1)
2H-i-' for a flux condition at x=H.



On the other hand, for a flux condition imposed at x=O, the Green's function is given by















Z't) =





(2n-1) r for a temperature condition at z=H, nff for a flux condition at x=H. H',


(3.51)


and


0, for a temperature condition at x=H,


Co = {


1
H7',


for a flux condition at x=H.


In practice, when the Green's functions are substituted into Eq.(3.46), a tedious and time consuming algorithm will be necessary to evaluate the integral because of the repetitive computation of the infinite series. However, by assuming a constant interface velocity from to to tj, Eq.(3.46) can be simplified to the following forms (see Appendix A for derivation):


00
V ( X , t ) = 2 cnL = B(0 t ) s i n(nx ) for a temperature condition imposed at x =0, and Y ) = C L dSAt + 2L 0C


(3.52)


(3.53)


where


On








39

for a flux condition imposed at z=O. In these equations,



B.(t1) = 2 2v (a.)sin(13s)- (f.v)cos(f.s,)



-e-k-n' [(ak.)sin(I.So)-( .v)CoS(1.3So)} (3.54)



and

Cn(t1) = (22 V ) (aI )COS(#.S1) + (I3nv)sin(/3,S1)


- en [/ cn)COS ) + (nnV)sin(,3nSo)] } (3.55)



Here, in both equations, So=S(to), Sl=S(tj), and v=dS/dt.



3.3.2 Solutions in cylindrical coordinate system

U-problem:

The governing equation for the U-problem in a cylindrical system can be obtained from Eq.(3.33) by setting d equal to unity. An implicit finite difference solution for this problem can be derived as


+-)Uj_1+ (1+29_)U '+1= U! if x 0 (3.56)
( 0 ,,,--O' + O1) U j + 1 U 3 , and

(1 + 40x)U,+ _4OxUn+l=U, if x-O (3.57)



Here, O'x=aAt/(2xAx) and Ox has been defined previously. These

equations can then be solved by using the same tridiagonal solver mentioned previously.










V-problem:

The solution of the V-problem can be obtained by setting d equal to unity in Eq.(3.42). Then, t1
v(,t C Srs(r)Xt Irr (3.58)
d-
0


The Green's function for a solid cylinder (O


G(z,tj I z',r) = - E-n - )3.x) (3.59) RC n=1 1),


where J. and J1 are Bessel functions of the first kind of order zero and one, respectively. The eignvalues /3,, are roots of



J0(IO.R) = O, n = 1,2,3... (3.60)



3.3.3 Solutions in spherical coordinate system

U-problem:

The governing equation for the U-problem in a spherical system can be obtained from Eq.(3.33) by setting d equal to two. An implicit finite difference solution for such a problem can be derived as



(-20)U_ + (U' if x0 (3.61) and

(1+60 )U''-6&0U .,., = if z=O (3.62)


J J J


X /








41

where O and 0, have been defined previously. These equations can also be solved by using the tridiagonal solver.

V-problem:

For a problem in spherical system, the solution of V can be obtained directly from Eq.(3.42) by setting d equal to two. This gives t1


V(x,t1) =L dSjT)S2(r)G(XtlIS(r),r)dT (3.63)
r=t0


The Green's function for a solid sphere (O
o

G(z, tI x', r) = x(Ee-ak3n(tl-T)sin(,3nx')sin( nx) (3.64) n=1


For a temperature condition imposed on the spherical surface, 1_2 #nwn N-Rs' =R (3.65) RS


On the other hand, for a mixed condition imposed on the surface,

1 2(_n+K_2)
- =sn +(3.66)



where 0,ns are roots of jncot(nRs) + K = 0 (3.67) and K-h 1
k Rs








42

The solution for the V-problem can also be simplified to the following form (see Appendix B for derivation):


V(X, t L D(tl) ..
n=1


(3.68)


where


v Slsin(flnS1)- - e 'sn ..... j Dn(ti) = (02)2 + V)
n 2
(I v) Sicos( ,S ) - S ck-to 2.t.


[(,#2)2 + 2-)V] 2] 2 - (lv) [sin(OnsS) - e n Ati(/S)

-23l3 cos(%nSi) - e n 2COS ) }


(3.69)


3.3.4 Multiple interfaces

The method developed above can also be of problems involving multiple interfaces. In of the U-problem remains unchanged. However, problem can be obtained by using multiple equivalent governing equation for this problem


applied to the solution this case, the solution the solution of the Vsources and sinks. The is


0V 1 a _dOV+ M LdSP 6(XS(t)) 19t - d 0 1X + _ c dt


(3.70)


where M is the number of interfaces. The general solution of this equation is















V(X't1) = ZLL P r)Sd(TGXI()~T(.1
P= (C dr ) (r)G (x, t, I Sp(r), ) dr (3.71 )
0 to


Here, the interface positions are determined by solving the following nonlinear equations simultaneously




T - U(SIi ti) = L ti dST(7)Sd( r)G(Sl,ilIS (T),r)dT
Trt"


Tm -U(S2,t1) = L ti dS4r)Sd(r)G(S,, t, I SP(r),r)dr



(3.72)


M L dS P(7) d
Tm -U(SM, t1) c -s -SPr)G(S , t1 fSgr) ,r)d-r





3.3.5 Stability analysis

Stability and convergency are essential to numerical solutions. For most problems, the stability and convergency tests are limited to the analysis of linear equations, which, hopefully, provides trends for the projection of nonlinear equations. An alternative to this approach is the use of a numerical experiment. In this subsection, the stability of the linearized Stefan problem in one dimensional Cartesian system will be analyzed. Numerical tests for stability will be relegated to Chapter 5 where concrete examples will be examined.












Following the Von Neumann stability analysis, a small disturbance of T is defined at time level n so that Tn = n+ Jn (3.73) where the hat '-' denotes the exact solution and c is the disturbance. In the superposition method, the temperature at a new time level n+1 can be decomposed as


Tn+l = Un+l + Vn+l (3.74) Then, the disturbance at the new time level becomes ,n+l = n+1 +n0+l (3.75) where c comes from U and C2 comes from V. As for the present work, the solution of U is obtained by using a finite difference method. Thus,



n+1 n
E1 = 1+20(1-cosA) (3.76a) if the simple implicit method (3.45a) is used, and ,n+1 - 1-O(1-cosA) n (3.76b)
1 1+0(1-cosA)



if the Crank-Nicolson method (3.45b) is used (Anderson et at. 1984). Here, in both equations, A is the phase angle.








45

As for V and S, Eq.(3.41) must be first linearized so that the Dirac delta function b(x-S) can be treated as constant. Then, if V and S are represented as



V=V+2, S=S+. (3.77)



The disturbance C2 must satisfy the following equation

1 aC2 02 2 Ldy 6(xS) t 2 acdt (3.78)



which is of the same form as (3.41). The disturbance c2 can then be expressed in terms of 7, the deviation in the interface position, as


t0




As discussed previously, for small At, dy/dr can be treated as constant. It can thus be taken out of the integral as t1


C T G(x,tlIS(r),r)dr (3.80) 1-=t
0


It is noted that, in this work, the V problem has been solved by using a homogeneous initial condition, the disturbance in V thus will not accumulate. However, the disturbance in U can be transferred to V through S which is solved using Eq.(3.43).








46

It is clear that, at the interface,


V(S,t1) = Tm - U(, tl) (3.81) Hence,
C IX= -E (3.82) Introduction of (3.82) into (3.80) yields ti
d-1 G (9, tjIS(r), r) dr (3.83)


r~0


or

d - t1 (3 .84 ) LJ G(S,tljS(r),r)dr
i-=t0


So that Eq.(3.80) can be written as C2 -H1 (3.85) 1R is defined as t1

f G(x,tu S(r),-)dr 11= 0 (3.86) f G(S,tljS(r),r)dr
0

The disturbance at the new time level can then be expressed as ,n+1 n(3(1.8) n [1+20(1 - cosA)] (3.87a)








47

for the simple implicit method, and



cn+1 1-O(1-cosA) 1).
OEcosA) (1 (3.87b) 1+-0(1-OA




for the Crank-Nicolson method. Finally, the amplification factor, defined as the ratio of c and cn, can be related to II as



A = En+l= (1 -11) (3.88a) A = = [1+20(l-c osA)]


and

A - -(1-cosA) (1-H (3.88b) 1+O(1-cosA)


for Eqs.(3.45a) and (3.45b), respectively.

Since




1+20(1-cosA) , 1+0(1-cosA) 1



for all 0 and A, the stability condition hinges on the value of 1-i11. A II versus x plot is thus given in Figure 3.1. It is found that the integration of the Green's function in the numerator of Eq.(3.86) peaks at x=S+7. If 7=0, the denominator of Eq.(3.86) is maximized so that H is always less than or equal to one. On the other hand, if 7 is not equal to zero but remains small, H is peaked at a value that is slightly greater than one. However, it is independent of the time step. In any case, the condition of 11-H[<1 is met. It can thus be

















1.2 1.0 0.8 0.6





0.4 0.2


0.2 0.4 0.6 0.8

Position, x/H


Figure 3.1 H versus x plot








49

concluded that the amplification A can never be greater than one, providing assurance of an unconditionally stable solution.

When a numerical scheme is developed for the solution of the Stefan problem, it is necessary that the numerical solution converges to the true solution. For numerical solutions using finite difference methods, convergence analysis is usually accomplished through the application of the fundamental equivalence theorem of Lax. According to this theorem, a numerical scheme is considered to be convergent if and only if the scheme is stable and consistent. In the superposition method, the U-problem is solved by a finite difference method, while as the V-problem is solved by an analytical method. The numerical scheme will be stable and consistent If both the U- and the V-problems are stable and consistent. It has been well established that the finite-difference equations, (3.45), for the U-problem are consistent with their differential-equation counterpart in (3.33). The solution of the V-problem is exact up to Eq.(3.43). In (3.44), the velocity of the interface is approximated as a finite-difference. It is obvious that this equation is also consistent. Thus the entire solution of the T problem is consistent with the original governing equations. As it is also stable for all time steps, convergence is thus assured.

It should be noted that the above analyses have been carried out on the basis of linearized equations. The nonlinearity effects remain to be investigated. Use will thus be made of numerical experiment later to test the stability (Chapter 5).








50

3.4 Extension to Multi-Dimensional Stefan Problems

The analyses developed above for the solution of one dimensional Stefan problem will now be extended to the solution of multidimensional Stefan problems. Again, the phase change is considered to be diffusion driven in a PCM having constant and equal properties in different phase regions.


3.4.1 Solutions in two dimensional Cartesian systems

For phase change in a two dimensional Cartesian system, the governing equation, initial, and boundary conditions are, respectively,


1aT _ a2T + 2T t t
a at +O2T t>t'0 (x,y) EfQ (3.89) T(x,y,to) =To(x,y) (3.90) aTbi
kii + hi[ Tbi - T = qj, i=1,2,... (3.91)


Following Patel (1968), the interface position (Figure 3.2) can be expressed in the y direction as


S = S(:,t) (3.92) Then the interface temperature and flux conditions become


T(x,S,t) = Ts(x,S,t) = Tm (3.93)


18 +ka 1[T x ,t T,(,S ) pL Sxt (3.94)


The governing equation for the equivalent problem becomes































. . . ' , . . . . . , . . ' . . . . . . ' . . . , . . ' . . . . . . . .. . . . .l . . . : . . . . . . r . , .
. . .. . . . .. . . . . .. . - . . . . . . . . . . . . . . . . . . . . . . . . . ..: . . . . . 1 : . . : L . . .
..........T: ' ... .. :: : , , : ' . ,: , .". . v ......... ... . '. . . .
:<''::' ::I ::" :' :'I:: ::''::I : ' ': ':i''::::' .' :''::: ':.::' S o i d :, . ' . ..]., '. ...


... .. . . .. . . ....uid

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
0 . . . I- . . . . . . . . . . ~ . . . . . .H . . . .


Figure 3.2 System analyzed in a two-dimensional
Cartesian coordinate system












l0T =2T 02T L O!y
0t X2 + 6(y-S(x,t) )



Application of the superposition method



T(x,y,t) = U(x,y,t) + V(x,y,t)


the solution of two problems, for which


1IaU a2U +!21U tto (y) a at - O a2 +y2 ' t 0 (Xy) E



U(x,y,to) = To(x,y)


the U-problem


(3.97)


(3.98)


i=1,2, ...


(3.99)


For the V-problem, the governing equation is


10 V 02V + 2V L Oy a at 2'a2 + 6(y-S(z,t))a


t>to, (X,y)EQ


The initial, boundary, and interface conditions for the V-problem are


V(x,Y,to) = 0


kaVbi+hVbi=
k,-bi +h~i=0


i=1,2, ....


(3.101) (3.102) (3.103)


(3.95)


leads to satisfies


(3.96)


(3.100)


k, bi + hi b - T.o] = qi'


Vt(x, S, t) = VS (ZS, t) = T.- U (X,S, t)








53

As discussed previously in Section 3.2.2, the U-problem can be solved by FDM. For a two dimensional problem, it can be solved efficiently by an alternating direction implicit (ADI) method which has been established to be unconditionally stable for all time steps. In this method, the temperature in x- and y-directions are calculated alternately. The finite difference equations are


Sx i-l 1l_, J + (1l+20x)--6,junl2- il,jnl2= _v m -1 +(1-20Y)U!i -+O T!'~ l
-0 U!+ n1l n+l 0Un31 + (1-2 )+l/2 + T,+l/2


-- + (1+2%Y)U9+!- =0U ,'i~+ (129x)Ui2�-kTn+1,2

(3.104)



where -- and = 2Ay. These finite difference equations can be
2Az2 ad 2Ay2
solved by using two tridiagonal systems of linear algebraic equations by means of an efficient tridiagonal solver.

The solution of the V-problem can be expressed again in terms of a Green's function as


t1 H


r=t0 x-0

The interface position curve S(x,tl) is then determined implicitly by the following equation


tl H
Tm-U(X,S,tl)=i drj SC(X,S,tllx',S(x',r),r)dx' (3.106) r=t0 X=O












In this effort, since the interface S is a function of both x and t, Eq.(3.106) can not be solved directly for S. The difficulty lies in the fact that a two dimensional searching for the position is necessary if this position is to be found iteratively. Alternatively, an assumption of the S function to represent the interface position will be necessary if this position is to be found analytically. Moreover, in the analytical solution, the accuracy of the result depends greatly on the fitness of the function that is initially assumed for approximation. In any event, the interface position can not be determined easily.

To overcome these problems, a new approach is developed by splitting the two dimensional problem into an array of transversely related one dimensional problems. In this approach, the second partial differential of V with x in (3.100) is changed to a finite difference as
2V(z,y,t) V(x+Ax, y,t)-2V(x,y,t)+ V(x-Ax, y,t) (3.107) OX, Ax2


Then, by indexing x position with subscript i so that V(x,y,t) is changed to Vi(y,t) and V(x�Ax,y,t) to Viil(y,t), Eq.(3.100) can be re-written as follows


1 aV(Y, t) _ a2Vi(y, t) Vi+l(y, t) - 2Vi(y, t) + Vj-l(y, t)
a t - y2 + Ax2



6(y-Si(t)), i= 1,2,... (3.108)








55

Equation (3.108) represents a group of piecewisely related one dimensional problems which can be solved by means of the Green's function using a time marching scheme. In this effort, the last two terms on the right hand side of Eq.(3.108) can be collectively taken as source terms. Then with the help of zero initial condition and homogeneous boundary conditions imposed in the V-problem, the solution of V can be written as


t1
v (y~t C:--,)(

r=at0
ti D
+a-f drif gi(y',r)Gi(y,tl y',-r)dy' (3.109) r-tO y'=O


where


gi(y',r) = V , +i(y', 7) - 2Vi(y',7) -+ Vj-I (y', 7)(3 10
Ax (3.110) AX2


It is still not possible to solve Eq.(3.110) because of the

function gi(y,t) which depends on Vi(y,t). Some additional work is necessary, and one way to do it is to expand g about to as



Ogi(y, to)
M(Y,T)= M(Y,(to) r-to)+..., tO_ r


If higher order terms are neglected, then g can be expressed as



Og(y t o)
M3Y,T7) = M,yto)+ at (r-to) (3.112)








56

Using Eqs.(3.101) and (3.110) permits the first term on the right hand side of this equation to drop out. As for the second term, it can be shown from (3.110) that


9g1(Y't0)= - [nVi(y~to) .2V.(Y'i0)+ aV,..(y~t0) (3.113)
Ot AX2 [ t - t Ot ]1


Then with the help of the initial condition given as (3.101), Eq.(3.100) can be used to re-write (3.113) as


A 2 6y- -s.
a t C a-bt a,- '~ l - 2 tt6 ( y - ) t- ( l_1) / ( 3 . 1 1 4 )




It thus follows that

t1

V(Y,tiL)f -- G1(y,tllSi,r)dr



+ A5 (Wi+1-2wj+wj.._) (3.115)



where

t1
W f ( -O Sa,(o)l~ iIS
w = jar-to)-t- " i(y,tlS(t),r)dr, K=i,i�1 (3.116)
r~to


In practice, a first-order approximation is sufficiently accurate to solve the problem at hand. For the first-order approximation


(3.117)


gi(Y, 7) - MY(', to)













Then, from (3.101),



gM(y', r) = 0 (3.118)



and V can be obtained by solving a simple equation as t1
Vi (y, tj)=L c Gj(y,tjjSj,r-)dr (3.119)

.r=t0




3.4.2 Solutions in cylindrical coordinate systems

The governing equation for heat diffusion in a two dimensional cylindrical system can be written as follows


18T 02T 18 O{rT
aT = 'a - t > tog (z,r) E (3.120)



As is always the case, appropriate initial and boundary conditions as well as the interface conditions must be specified. For the phase change in cylindrical coordinates, two types of the interface motion will be addressed-the phase-change interface may be moving in zdirection (Figure 3.3), and the phase-change interface may be moving in r-direction (Figure 3.4).

In the former, the interface position may be expressed as



S = S(r,t) (3.121)















11~~ I


.v.' , .. . . v ......


Solid s


Figure 3.3


System analyzed in cylindrical coordinates: Interface moves in z-direction


I r


R--Liquid


-Solid


0


Figure 3.4 System analyzed in cylindrical coordinates:
Interface moves in r-direction


Liquid


--A


.J \


A








59

Then, the interface conditions are written as


T(S(r, t), r, t) = Tm


(3.122)


1(OS 2] [T aT] _piLl aS(r,t) E -a7 I z -Oz Z=s k ot In the latter, the interface position may be expressed as


S = S(z, t)


Then, the interface conditions are given by


T(z, S(z, t), t) = Tm


(OS2I [aTs aT1 _pI LI OS(z,t) [l+azj [IFar _ -=s k at


(3.123)


(3.124)


(3.125)


(3.126)


In the method of superposition, the temperature is decomposed as


T(z, r, t) = U(z, r, t) + V(z, r, t)


(3.127)


Here, the U-problem satisfies the same governing equation for T as in Eq.(3.120), and the nonhomogeneous boundary and initial conditions. This problem can be solved by a finite difference method as



Uit' - U J = 20z(Ui+j-2Ui, +Ui-lj)


+ Or'(Ui, j+1 - Ui, j1) + 20r(Ui, j+1 - 2Ui,j + Ui, j-1)


(3.128)


where 0;, aAt ,0, aAt and or= 'At
2Az 2 2rAr r 2Ar2








60

Equation (3.128) can be solved either by an explicit method by setting U=U" on the right hand side, or by an implicit method by setting U=Un+l on the same side. In any event, the solution is routine.

The V-problem can be solved by the Green's function method. Again, the interface flux condition is incorporated into the governing equation; Then, for the case of interface moving in the z-direction, the governing equation becomes



lg 02V lOI )+ L a
aPt r az rZr r(z-S(r,t)) (3.129)



The governing equation is split into two directions, where a finite difference method is used in the r-direction. This reduces the V equation to the following IaVj =02VJ L OZ (a = - a (z-Sj(t)) +gj(z,t) (3.130)



where
Vj+1- Vj_1 Vj+ -2Vj +Vj1
2rAr + Ar2 for r#O gj(z, t) = 4(V +1 - Vj) for r= (3.131) Ar2

As described in Section 3.4.1, the first-order approximation of gj leads to a very simple approximation as gj(z,r) - gj(z, to) = 0 (3.132)










The V-problem can be solved as



V =ZLTl = i Gj(z,tlISj,,r)dr (3.133) r=t0


Similarly, the second-order approximation of gj results in the following solution


t
Vjztl= -j'Gj(z,tl1Sj,r)dr+Wj (3.134) .r=t 0

where

Wj+1 - W 1 wj+1 - 2wj + Wj_f 2rAr + Ar2 for r#O Wj= w ) (3.135) Ar 2 for r =0 and

t1
w " = a s (-t � ) G j( z t llS ,,(t ),r )d -r ' c = j ,j~ l( .1 6

-r=to


For the case of the interface moving in the r-direction, the governing equation of the equivalent problem is given by



1a8 - =2 ---raVY)+k. L 6(r-S(z, t)) (3.137) aDt az2 ~r8 D



Again, the governing equation is split into two directions; this time, a finite difference method is used in the z-direction. The V equation becomes














1OVi(r, t) 1 _o_[_o t)(+,Lt)a ] (r-S(t)) + g,(r,t) (3.138)




where

gi(r,t) = V,+l(r ,t) - 2V,(r, t) + Vi-,(r, t) (3.139) Az2


The first-order accurate solution of the V-problem is readily obtained as t
Wi(r,ti) -=LC 'O-Si(,)Gj(r, tltISi, 7-)d7- (3.140) .r=t0


and the second-order accurate solution can be derived as


t1
Vi(r, tl) --Lf -S i(7)Gi (r, tjI Si, 7)d7 T~t0

+ otL (Wil2ww1) (3.141)


where
t1
= f (r-to)aSV(to)S.(to)Gi(r,tl S,(t),r)dr r=to

=i, i�1 (3.142)





3.4.3 Solutions in three dimensional Cartesian systems

The governing equation of a Stefan problem in a three dimensional Cartesian system takes the form














18T 02T + a2T + a2T
UOt -2 Oy2 -Z2, t>t0, (X,y,z) E (3.143) The interface position is expressed as



S = S(x, y, t)



Then, the interface condition can be written as (Patel,1968) +(aS\2 +(taS 21 [aT s aT] pILI aS S) Vay) j L Oaz Ozj - - k at(3.144) With the use of the superposition method, the temperature is decomposed as


T(x,y,t) = U(z,y,z,t) + V(x,y,z,t) (3.145) in which the U-problem satisfies a governing equation that is of the same form as (3.143) and relevant nonhomogeneous initial and boundary conditions. The finite difference equation for the U-problem can be derived as


U.+tI - = 20x(Ui+1,, ,j - 2Uj, jk+ U -,j, k) " 20y(Ui, j+1,k - 2Ui, j,k + U, j-1,k) + 20z(Ui, j, k+1 - 2Ui, j, k + Ui, j, k-1) (3.146)













where O, O., and 0. have been defined previously. This equation can be solved by using an iteration procedure such as Gauss-Seidel and SOR method, or by a modified ADI method.

Again, the V-problem can be solved by a Green's function method. In this method, the equivalent problem can be written as



10 2V 02V+a2V+L Lz b(zS) (3.147)
-a aX2 ay2 az2 a



Then, by using a finite difference method in both x- and y-directions, Eq.(3.147) can be re-written as


Ly,,, -2Vi, +L L t(z-S) +g, j(z,t) (3.148)
aat acz2 ct

where



gi, j(z, t) = Vj+l,j-2V,i, +Vi,I V, j+1 -2VI, j+Vi, j-1 (3.149) Ax2 Ay2


With the help of the homogeneous initial condition for the V-problem, the first-order accurate approximation of gij leads to a very simple relation as

gi, j(z,r) = gi, j(Z, to) = 0 (3.150)



and the solution of the V-problem follows as


t(
Vi jZ t)- fa ijGi, j (z, tjI ISi, j, 7) dr (3.151) T=t0












In a similar fashion (see Section 3.4.1), the second-order accurate approximation of gi,j leads to a solution as


_i ass.,
Vi, jZ' t' = L (Z ztj I Si, j, r)d-r

r=t0

+ aL 2(Wi+l,,j_-2wi, j + Wi-1, j
cAX2



+ aL 2 (3.152) "} O"-- (W i, j+ 1 -- 2 W i, j - - W i _I, j)( 3 5 2 cAz2



where

t1
W r. n , n = ( -r - t �) a t G i, j ( z t j I S , (t ) r ) d r

"-:t0 re i�1, n=j,j~l (3.153)






3.4.4 Stability analysis

For multi-dimensional problems, the explicit treatment of the second partial differential may impose additional constraints on the stability of the solution. Due to the difficulties associated with the nonlinearity of the problem, the stability can not be tested by simple analysis. However, some comments can be made as follows.

It is well known that the solution of U is unconditionally stable. As for the V, since it always starts with a homogeneous initial condition, the disturbance (error) of V will not accumulate.








66

The only source of error that may be brought into the V-problem is U(S,t) which is used to find the interface position. In this connection, because of the use of the Green's functions, disturbance in the integrand will not be magnified but will be suppressed. On the other hand, the disturbance in the V solution can still be brought to the U-problem through the superposition principle. However, this disturbance can still be suppressed by the ADI method in the solution for the temperature in the succeeding time step. It is thus expected that any disturbance will not grow in the solution, and, as a result, the method is unconditionally stable. It is noted that the discussion above is not meant to be rigorous. Numerical experiments for a twodimensional problem (See Chapter 5) will be provided later for tests.



3.5 Application to Latent Heat Thermal Energy Storage Systems


For the phase-change problems solved in the preceding sections by the source-and-sink method, the boundary conditions have been applied directly to the phase change material. In practice, the phasechange material may be heated or cooled by a fluid which is in motion itself, a conjugate heat transfer problem thus results. In this problem, the fluid temperature is only given at the inlet. Both the wall and the fluid exit temperature are unknown.

Two conjugate heat transfer problems will be solved in Sections 3.5 and 3.6. Analysis for fluid flow in a straight parallel channel of phase change material will be studied first. Results for a sample study for this problem will be given in Chapter 5.










3.5.1 Problem formulation

A thermal energy storage (TES) unit is sketched in Figure 3.5. This unit can be modeled for analysis as shown in Figure 3.6. It consists of a straight parallel channel surrounded by PCM. Heat transfer fluid (HTF) is forced to flow through the channel and it exchanges heat with the PCM. It is a time-dependent phase-change problem combined with conjugate forced convection. The governing equations for the forced convection in the HTF and the heat diffusion in the PCM must be solved simultaneously. The solution for heat transfer in the PCM has been covered in the previous sections. Only the convection problem and the matching method need be described in this section.

Analysis of forced convection for channel flow can be accomplished by solving the Navier-Stokes equations and the energy equation. However, a direct numerical simulation of the Navier-Stokes equations is time consuming. In addition, for a conjugate problem, an iteration is necessary to satisfy the boundary condition at the fluid/PCM interface, which further adds to the complexity of the solution. An alternative approach is taken that uses empirical correlations for the treatment of the heat transfer at the fluid side. This permits the use of convective heat transfer coefficient. While the accuracy of the method depends on the accuracy of the empirical correlations, the uncertainty is comparable with that introduced in the turbulence model. With this effort, the fluid temperature in the channel can be modeled as a function of position and time and analyzed by solving the following equation:



















HTF outlet Container


Figure 3.5 A schematic drawing for thermal energy storage unit







Y


H
PCM
Plane of symmetry




-- =--- 5 -z -- zz50 Plane of symmetry
HTF Wall


Figure 3.6 Simplified system for analysis














aT1 a + h(T1 Tw)
-f + i- Tf +-f- ---- =0 (3.154)




At the inlet,

T1(0, t) = Ti,(t) (3.155)



In these equations, Tf is the mixed mean temperature of the fluid, U is the mean velocity, h is the convective heat transfer coefficient, d is the half width of the channel, and p1 and cf are the density and the specific heat of the fluid, respectively. In the formulation above, the Peclet number is taken to be large, so that the heat diffusion in the x direction can be neglected.

In the conjugate problem, Eq. (3.154) is to be solved in conjunction with the phase-change problem whose solution has been given previously. If the heat conduction in the axial direction in the channel wall and the heat storage in the wall are negligibly small, heat transfer by convection from the fluid to the wall must be equal to the heat that is supplied to the PCM through its lower boundary. Hence,
h(Tf - T.) A201 q_ (3.156)




where yp is the lower boundary of the PCM (see Figure 3.6). Using an electric analogy, the temperature at this boundary can be related to the temperature at the inner surface of the channel wall by the relation














TP =TW -bwqw/kw (3.157) where bw and kw are the thickness and the conductivity of the wall, respectively.

The heat transfer coefficient h is taken from the literature. In this regard, the Nusselt number for a fully developed laminar flow varies from 7.54 for constant wall temperature to 8.24 for constant wall heat flux condition (Kays and Crawford, 1980). There is a lack of correlation for turbulent flow in a straight channel. Use is thus made of an empirical equation developed by Gnielinski (1976) as


(f /8)(ReD - 1000)Pr
NuD 1 + 12.7(f/8)"/2(Pr2/3- 1)



where f is the friction factor, which is related to the Reynolds number as

f = [0.79n(ReD)- 1.64]-2 (3.159)



In these equations, both the Nusselt number and the Reynolds number are defined as


hDh UDh (3.160) NUD --- = , R(D3



where Dh is the hydraulic diameter of the channel. It is noted that Eq. (3.158) is strictly valid for a circular pipe with any condition (constant wall temperature or flux) imposed on the boundary. The








71

Prandtl number is in the range of 0.5 to 2,000, while the Reynolds number is in the range of 2,300 to 5x106. It is adapted for use in the straight channel by redefining the dimensionless numbers in terms of Dh as shown in (3.160).



3.5.2 Solution method

In the conjugate problem, the heat flux at the channel wall and the temperature in the PCM and the fluid are all unknown. To satisfy the matching condition (3.156), it is necessary to use an iteration. However, as discussed by Barozzi and Pagliarini (1985), even for a simple steady-state conjugate problem without phase-change, instabilities may be encountered during the iteration process. Convergent results can only be obtained if the initial guesses are sufficiently close to the true values.

To alleviate these problems, the heat flux in the third term in (3.154) is re-written in terms of T1, the temperature at the first grid point above the boundary of the PCM at yp (Figure 3.6). Thus,



h(Tf - TW) T"T1 (3.161)



where R" represents the overall thermal resistance per unit area,



R" 1 bw Ay (3.162) h ks kn



It is noted that, in the numerical solution, if the interface position








72

S is less than Ay, which is the grid spacing in y direction, then the T1 in Eq. (3.161) is replaced by Tm, the phase-change temperature, and the Ay is changed to S.

The energy equation (3.154) can then be written as


OT, OUT T1-T1 5t Ox + PcjRd -


(3.163)


This equation can be solved by using a first-order accurate implicit finite difference scheme as


-n +0 + ,n+ +0T
n l = T j - "1 + f j 1 2 T 1
-f j = 1+--01+"-02


(3.164)


where


0 UAt
-'AX'


At
c= PfCR"d


The temperature imposed at the surface of the PCM follows as


TP = fRR


(3.165)


where


= + b


(3.166)


" - Ay
R2 - k-


The surface temperature given in (3.165) will then be used as the boundary condition for the solution of the phase-change problem.













3.5.3 Energy conservation test

For the solution of the conjugate problem given in this section, errors mainly come from the matching conditions between the PCM and the HTF. Since experimental data are usually unavailable for comparison, the accuracy of the numerical solution developed in this section will be tested by means of an energy conservation that is based on an open-system thermodynamic analysis.

For the heat transfer of channel flow between parallel walls of phase change material without heat and work interaction with the surroundings, the first law of thermodynamics can be written as



Ein -Eout= AEcv/At (3.167)



where the ki refers to the energy rate associated with the fluid flow entering the channel and k.,t refers to the energy rate associated with the fluid flow leaving the channel. The term on the right hand side represents the time rate of change of the internal energy inside the control volume. For the problem under investigation, there are energy changes in the PCM, the wall, and the fluid.

For a time period of t, the left hand side of Eq.(3.167) can be expressed as

t

Ein - Eout f J Qpfcf[Tf(O, r) - Tf(H,r)Idr (3.168) r=0















while the right hand side can be expressed as


M D H


H
+ f [-LpS(x,t)z
X=O

H' (3.169)
+ f Pcb4Tt(k' t) - T0Jdx
x=O




x=0


Here, the energy balance in (3.167) is formulated on the basis of discharge of heat from the PCM. The PCM, the wall, and the fluid are at the same initial temperature T., and the channel has a width W. Each PCM plate has a thickness 2D. Under these conditions, the first term on the right of (3.169) refers to the energy associated with the change of the sensible heat of the PCM. The second term refers to the energy associated with the phase change. The third term accounts for the energy change in the containing wall, and the last term refers to that energy associated with the change of the temperature of the HTF inside the channel. For an accurate numerical solution of the problem, the left hand side of (3.167) should be equal to the right hand side of the same equation. This energy conservation will be used to check the accuracy of a numerical experiment to be presented in Chapter 5.













3.6 Application to Encapsulated PCM TES Systems

In the preceding section, the phase change material analyzed is in the form of a plane slab. Since the heat transfer fluid is flowing over the surfaces of the slab, the plane of symmetry can be taken to be the half thickness of the slab as shown in Figure 3.6. An alternative form of the phase change material is one that is encapsulated in spheres as shown in Figure 3.7. This section is devoted to the analysis of the phase change in packed spheres.



3.6.1 Problem formulation

The system under investigation is sketched in Figure 3.7. A tank is packed with encapsulated PCM spheres. Heat transfer fluid enters through the bottom of the tank, exchanges heat with the PCM by forced convection, and exits from the top of the tank. It is assumed that, in the HTF, the temperature only varies in the z direction while, in the capsulated spheres, the temperature only varies in the r direction as shown in Figures 3.8 and 3.9. Other assumptions that are related to the phase-change analysis earlier still apply.

For the analysis of the heat transfer in the HTF, the control volume is taken to have a thickness that is equal to the diameter of the packed sphere (Figure 3.8). The energy equation for the fluid can be written in terms of the mean temperature Tf and mean velocity U as


V T+T1 + z h(T- Ts)=O, T= T (z't) (3.170)









































lquid SONd



Spherical capsule


Figure 3.7


A schematic drawing for the encapsulated latent heat TES unit


Number Of Waen.


tIMU


LI


.777











T f U


2Rs Tf, U





Figure 3.8 Control volume for fluid analysis





TS T
ST,




h . ..... . . .
T, Tp Ts~ Tf


Figure 3.9 Boundary conditions for PCM/HTF analysis












where, in addition to those that have been defined previously, T. is the surface temperature of the sphere, V! is the volume of the free space in a packed layer, and A. is the total surface area of the spheres that are in contact with the fluid. For a single sphere of radius Rs, the surface-to-volume ratio is given by the simple relation As = 3 (3.171) Vs R$

which also holds for multiple spheres in a packed layer. Thus, by using a void fraction p, defined as the ratio of the free volume V1 to the total volume Vc, Eq. (3.170) can be recast as


OT +OT +3(1 - )h (T (3.172) at OU"-z ppf c Rs ( -Ts) = 0(


This equation is equivalent to that derived by Schumann (1929) and Riaz (1978) for packed rock beds. For the present problem, the initial and boundary conditions for the fluid are given as



T1(z,O) = T,, T!(0, t) = Ti.(t) (3.173)



As for the spheres, the governing equation, initial, boundary, and interface conditions have been given in Section 3.3.3. These equations must be solved in conjunction with (3.172) for heat transfer analysis. As before, the solutions of the fluid and the PCM are to be matched at the surface of the spheres by the relation










h(T8-T=-k aT| = (3.174) h(s-T)= Pa-r =p=

In this effort, the convective heat transfer coefficient h for the packed spheres is obtained by curve-fitting the experimental data of Chou et al. (1992) as


Nu=11 + 0.06Pe (3.175) where, Pe is the Peclet number (Re.Pr), The Reynolds number is defined on the basis of the diameter of the sphere and the average velocity of the fluid as
U-(2Rs)
Re U (3.176)




The heat transfer coefficient can then be evaluated by means of the Nusselt number as


Nu.k (3.177) 2R$



where k! is the conductivity of the fluid.



3.6.2 Solution method

Solution of the conjugate problem for fluid flow in a packed sphere bed follows the same line described earlier for channel flow. In this case, it is assumed that the heat storage in the sphere wall is negligibly small, so that the heat transfer from the fluid to the sphere can be expressed in terms of the temperature difference between the fluid and the first nodal point that is located inside the surface










of the PCM as



h(T-Ts) = TT (3.178)



For the spherical PCM, the total thermal resistance based on the surface area of the sphere is



R(3.179)



where k. is the conductivity of the sphere wall, RP is the radius of the PCM, and R. is the radius of the sphere (see Figure 3.9). Equation (3.172) can then be written as OTf + u(TT 3(1 - ) T)= 0 (3.180)




As before, in the numerical solution, if the interface position S is greater than rl, the first nodal point position that is inside the surface of the PCM, but is less than Rp, the radius of the PCM, then T1 in Eqs. (3.178)-(3.180) is replaced by Tm and r1 is changed to S.

Equation (3.180) can be solved by using an implicit finite difference scheme as



T+ f Th f id 1 (3.a181) f 1 +01+ 02



where Tfi is the temperature of the fluid at the inlet of the control










volume as shown in Figure 3.8, and



= Atu - 3(1 - p)At
0- s 82- PPc " (3.182)



The temperature at the surface of the PCM can then be evaluated as


R'T1 +R2Tf (3.183)




where

R" R"-i (3.184) !R/s p s




This temperature is to be used as the boundary condition for the solution of the phase change in the PCM as described in Section 3.3.3.



3.6.3 Energy conservation test

As mentioned in Section 3.5.3, errors for the conjugate problem mainly come from the matching conditions between the PCM and the HTF. Since experimental data are usually unavailable for comparison, the accuracy of the numerical solution developed in this section will be tested by means of an energy conservation that is based on an opensystem thermodynamic analysis.

For the heat transfer of fluid flow in packed spheres inside a tank without heat and work interaction with the surroundings, the first law of thermodynamics can be written as














kin - F- out = AEcv/At (3. 185)



where the E, refers to the energy rate associated with the fluid flow entering the tank and Eout refers to the energy rate associated with the fluid flow leaving the tank. The term on the right hand side represents the time rate of change of the internal energy inside the control volume. For the problem under investigation, there are energy changes in the PCM, the shell of the PCM sphere, and the fluid.

For a time period of t, the left hand side of Eq.(3.185) can be expressed as


t
Ein - Eut = J QPc[ Tf (O,'r) - Tf (H, )]dr (3.186) T=O



while the right hand side can be expressed as

rRp

NZ{47rJf ptTkr ) -T ,,] rdr k=1 r=-O


+4[S3(t)R]



(3.187)
47r 3 3 +' Pwcw(Rs- R3Twk(t) - T.]


H
+irR2 PfC cT(z, t) -T z Z=O













Here, the energy balance is formulated on the basis of discharge of heat from the PCM. The PCM, the shell, and the fluid are at the same initial temperature T0, and the tank has a height H and a radius Rt. The number of spheres in the tank is N which is related to the volume fraction of the PCM sphere as

volume of the tank

N=(volume fraction of sphere) volume of a sphere



Under these conditions, the first term on the right of (3.187) refers to the energy associated with the change of the sensible heat of the PCM. The second term refers to the energy associated with the phase change. The third term accounts for the energy change in the spherical shell of the PCM spheres, and the last term refers to that energy associated with the change of the temperature of the HTF inside the tank. It is expected that, for an accurate numerical solution of the problem, the left hand side of (3.185) should be equal to the right hand side of the same equation. This energy conservation will be used to check the accuracy of a numerical experiment to be presented in Chapter 5.














CHAPTER 4
TWO DIMENSIONAL PHASE-CHANGE EXPERIMENT


In this chapter, an experiment is described for the purpose of validating the numerical solutions described in the previous chapter. The experimental setup, method of measurements, and data reduction are discussed in detail.



4.1 Objective and Basic Considerations

The objective of the experiment is to validate the numerical solutions developed in the previous chapter. Since there is lack of data for two dimensional phase change, experimental measurements become necessary. In this effort, the interface motion is visualized and measured. The measured interface position is then compared with the predicted position for validation of the numerical solution.

It has been well established in phase-change studies that the interface position provides the best indication for the solution error. Comparison of the temperature field is often inconclusive even if the interface position is conspicuously in error. This is due to the fact that, for a phase-change problem, the boundary and interface conditions can be satisfied exactly. This leaves little error for the temperature. As a result, the temperature provides little resolution for discerning the numerical solution error.










To facilitate visualization of the interface position, a test cell of the configuration of a rectangular parallelepiped is constructed. Paraffin wax is chosen as the phase change material because of its very narrow melting range or nearly distinct melting point (Lane, 1986). It also possesses the desired characteristic of becoming translucent to visible light upon phase change from solid to liquid state. The interface position can thus be measured by illuminating the wax with an incandescent light source and recording the interface position on photographic films by means of a camera. The thermophysical properties of the paraffin wax are summarized in Table 4.1 (source: Lane, 1986), in which, the melting point listed is the result of experimental measurements that are specifically performed for the study described in this work.



Table 4.1 Thermophysical properties of paraffin wax Properties Solid Liquid Melting point (K) 327 Heat of fusion (kJ/kg) 266
Specific heat (kJ/kg K) 2.51 2.95 Conductivity (W/m K) 0.24 0.24 Density (kg/m3) 760 818


In the experiment, the wax is heated electrically. The apparatus is designed so that the boundary condition can be controlled precisely. Provisions are thus made to simulate the prescribed temperature at the heated side of the wax while keeping the other boundaries of the wax insulated. Since no convection has been




Full Text

PAGE 1

SOURCE-AND-SINK METHOD OF SOLUTION OF TWO AND THREE DIMENSIONAL STEFAN PROBLEMS By HONGJUN LI A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1993

PAGE 2

TO MY WIFE YU LUO, MY DAUGHTER ANGELA, AND MY PARENTS.

PAGE 3

ACKNOWLEDGEMENTS I would like to express my sincere gratitude and appreciation to Dr. C. K. Hsieh , chairman of the supervisory committee, for his guidance and encouragement throughout my entire graduate study, especially for his patience in revising this manuscript. His enthusiasm and great scholarly attainments have contributed substantially to the quality of this work. I also want to extend my appreciation to Dr. Y. D. Goswami , cochairman of the supervisory committee, for his continuous support, encouragement and guidance during the research. I would like to thank Dr. W. Shyy, Dr. J. F. Klausner, and Dr. R. Mei for serving on my supervisory committee and for their fruitful comments and suggestions. I acknowledge my colleagues and fellow graduate students for their friendship and valuable help. Finally, my deepest gratitude goes to my wife, Yu Luo, for her unselfish love, patience, and understanding, and to my parents for their love, encouragement and support. Ill

PAGE 4

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS Hi LIST OF TABLES vii LIST OF FIGURES NOMENCLATURE xi ABSTRACT xiv CHAPTERS 1 INTRODUCTION 1 2 REVIEW OF LITERATURE 6 2.1 Analytical Solution of the Stefan Problems 8 2.1.1 Exact solutions 8 2.1.2 Approximate analytical solutions 9 2.1.3 Other approximate solutions 12 2.1.4 Analytical solutions of two dimensional Stefan problems 14 2.2 Numerical Solutions of the Stefan Problems 15 2.2.1 Classical formulation 15 2.2.2 Weak formulation 17 2.2.3 The boundary element method 18 3 NUMERICAL ANALYSIS 20 3.1 Problem Formulation 20 3.1.1 Governing equations 20 3.1.2 Brief review of the SSM 22 3.1.3 Time marching method 25 3.2 SSM With Superposition Techniques 26 3.2.1 The validity of the superposition method 26 3.2.2 Solution method 30 3.3 Solution of One dimensional Stefan Problems 35 3.3.1 Solutions in one dimensional Cartesian system . 35 3.3.2 Solutions in cylindrical coordinate system .... 39 3.3.3 Solutions in spherical coordinate system 40 3.3.4 Multiple interfaces 42 3.3.5 Stability analysis 43 i v

PAGE 5

3.4 Extension to Multi-Dimensional Stefan Problems 49 3.4.1 Solutions in two dimensional Cartesian system . 50 3.4.2 Solutions in cylindrical coordinate system .... 57 3.4.3 Solutions in three dimensional Cartesian system 62 3.4.4 Stability analysis 65 3.5 Application to Latent Heat Thermal Energy Storage Systems 66 3.5.1 Problem formulation 67 3.5.2 Solution method 71 3.5.3 Energy conservation test 73 3.6 Application to Encapsulated PCM TES Systems 75 3.6.1 Problem formulation 75 3.6.2 Solution method 79 3.6.3 Energy conservation test 81 4 TWO DIMENSIONAL PHASE-CHANGE EXPERIMENT 84 4.1 Objective and Basic Considerations 84 4.2 Experimental Setup and Measurement 86 4.2.1 The setup g 2 4.2.2 Test procedure 92 4.3 Sample Results and Experimental Uncertainty 95 5 RESULTS AND DISCUSSION 95 5.1 Numerical Results of One Dimensional Problems 95 5.1.1 Test of accuracy and stability 95 5.1.2 Subcooling effect on the interface position ... 105 5.1.3 Solution with oscillating interfaces 107 5.1.4 Solution with multiple interfaces 109 5.1.5 Melting of a sphere 114 5.2 Numerical Results of Two Dimensional Problems 116 5.2.1 Test of accuracy and stability 118 5.2.2 Long-time solutions 121 5.2.3 Comparison with experimental results 124 5.2.4 Test for coalescing fronts 126 5.3 Application to Two Dimensional Latent Heat Thermal Energy Storage Systems 128 5.4 Application to Encapsulated PCM TES Systems 134 5.5 Visualization of Numerical Results 141 6 CONCLUSIONS AND RECOMMENDATIONS 144 APPENDICES A DERIVATION OF EQUATIONS (3.52) AND (3.53) 148 v

PAGE 6

B DERIVATION OF EQUATION (3.68) 152 C FORTRAN PROGRAM FOR TWO DIMENSIONAL STEFAN PROBLEMS 153 REFERENCES 174 BIOGRAPHICAL SKETCH 18n VI

PAGE 7

LIST OF TABLES 2^±e Page 3.1 Expression to account for the effects of boundary conditions 23 4.1 Thermophysical properties of paraffin wax 85 5.1 Cases tested for one dimensional Stefan problems 96 5.2 Properties of PCMs analyzed 97 5.3 Cases tested for two dimensional Stefan problems 117 3.4 Input data for TES analysis 129 5.5 Specifications of encapsulated TES unit 134 vi 1

PAGE 8

LIST OF FIGURES Eig ure Page 2.1 Solution methods of Stefan problems 7 3.1 II versus x plot 4 g 3.2 System analyzed in a two dimensional Cartesian coordinate system 51 3.3 System analyzed in cylindrical coordinates: Interface moves in ^-direction 58 3.4 System analyzed in cylindrical coordinates: Interface moves in r-direction 58 3.5 A schematic drawing for the thermal energy storage unit . . 68 3.6 Simplified system for analysis 68 3-7 A schematic drawing for the encapsulated latent heat TES unit 7 g 3.8 Control volume for fluid analysis 77 3.9 Matching conditions for PCM/HTF analysis 77 4.1 Test setup §7 4.2 Assembly drawing of the test cell 88 4.3 Exploded view of the test cell 89 4.4 Exploded view of the heater 91 4.5 Sample pictures for interface position measurements 93 5.1 Temperature profile of the superposition method 98 5.2 Test of accuracy and stability using different At and Ax . 99 5.3 Errors in the interface positions for different At 101 5.4 Interface positions for time-variant temperature condition 102 vi i i

PAGE 9

5.5 5.6 5.7 5.8 5.9 5.10 5.11 5.12 5.13 5.14 5.15 5.16 5.17 5.18 5.19 5.20 5.21 5.22 5.23 5.24 5.25 5.26 5.27 Interface positions for time-variant flux condition 102 Interface positions in finite domain 104 Temperature history at the insulated wall ( x-H ) 104 Interface positions for different subcooled conditions.... 106 Test of oscillating interfaces 108 Sketch of double interfaces generated by conditions imposed on two boundaries HO Temperature profiles of a double interface problem Ill Interface positions of a double interface problem 112 Sketch of double interfaces generated by conditions impose on one boundary 113 Interface positions of a melting-freezing problem 113 Interface positions of sphere melting 115 System used in two dimensional stability test 119 Interface positions at different times at x=0 and x-H /2 .. 120 Interface position curves at <=3000s for different At .... 120 System used for long time solution 122 Comparison of numerical and analytical solutions for a two dimensional problem at long time 123 Comparison between numerical results and experimental data for two dimensional phase change 125 Interface position curves in a two dimensional coalescing-f ront problem 127 Test of energy conservation 131 Water exit temperature at different flow rates 132 Isotherms in a half PCM plate 133 Test of energy conservation 136 Water exit temperature at different flow rates 137 IX

PAGE 10

5.28 Water exit temperature for different ball diameters 138 5.29 Water exit temperature for different shell thickness 139 5.30 Temperature profiles in the PCM and the shell 140 5.31 Post-processor for a two dimensional TES unit 142 5.32 Post-processor for an encapsulated TES unit 143 x

PAGE 11

NOMENCLATURE A Amplification factor, area b Thickness of TES containing wall B,C Constants D Constant, dimension of PCM in y-direction Dj j Hydraulic diameter c Specific heat d Weighting function, half width of channel E Energy / Friction factor G GreenÂ’s function g Function defined in Eq.(3.110) etc. H Dimension of PCM in ^-direction h Convection heat transfer coefficient k Thermal conductivity L Latent heat Nu Nusselt number n Normal direction J Bessel function Pe Peclet number Pr Prandtl number Q Volume flow rate xi

PAGE 12

9 Heat flux R Radius r Coordinate in cylindrical and spherical systems R” Thermal resistance Re Reynolds number S Interface position Ste Stefan number T Temperature t Time U Decomposed temperature, mean velocity V Decomposed temperature, volume W PCM width x,y,z Coordinates in Cartesian system Greek a Thermal diffusivity p Eigenvalues in the Green’s function 6 Dirac delta function e Small constant e Parameter in finite difference method n Void fraction n Function defined in Eg. (3.86) p Density T Dummy variable of time n Solution domain xi 1

PAGE 13

Subscripts and Superscripts 6 Boundary value cv Control volume / Fluid i Boundary index, spatial index j Spatial index l Liquid phase m Melting point n Time level index p Phase regions, PCM s Solid phase, interface position, sphere w Wall, shell xi 1 1

PAGE 14

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SOURCE-AND-SINK METHOD OF SOLUTION OF TWO AND THREE DIMENSIONAL STEFAN PROBLEMS By Hongjun Li December 1993 Chairman: Dr. Chung K. Hsieh Cochairman: Dr. Yogi D. Goswami Major Department: Mechanical Engineering A source-and-sink method has been developed for the solution of conduction heat transfer with phase change in multiple dimensions. In this method, the heat transfer in one direction is decoupled from the other directions by changing the second partial differential in these directions where the phase change is less dominant to a finite di-ff erence form and solving the problem in the other direction with an analytical solution that accounts for the motion of the interface in a phase change problem. The solution developed in this work is thus independent of the equations used to represent the interface as well as the conditions imposed on the boundaries. The method has been applied to the tracking of single, double, and coalescing interface fronts formed by different phases assuming xiv

PAGE 15

equal properties. The method has been demonstrated to be accurate, convergent, and stable by numerical computation as well as experimental measurements. The method has also been applied to the solution of conjugate heat transfer including convection in a fluid medium and conduction in the phase change material, a problem applicable to the analysis of latent heat thermal energy storage systems . xv

PAGE 16

CHAPTER 1 INTRODUCTION Phase-change problems involving melting and freezing have found numerous applications in science and engineering. Phase-change processes are encountered in manufacturing and material processing, and in annealing and continuous casting. In the important area of food processing, phase change has been an integral part of a freezing and thawing process in which one substance can be separated from another by means of a phase change. The latent heat absorbed and released by a substance during a phase change also works favorably in energy storage. Since the latent heat is usually large for some substances, a small volumeto-energy ratio can be achieved. This makes the process particularly attractive for solar energy applications. Recently, with the increased use of lasers in industrial processes such as cutting, drilling, surface hardening and densif ication, the phase change has added another dimension to contemporary research. Phase-change problems have been commonly referred to as a class of moving boundary problems. In contrast to many conventional physical problems in which the solutions are sought in a fixed domain, the phase-change problems are solved in a variable domain. 1

PAGE 17

2 There are two phases in the medium which are separated by a phasechange interface. Under the continuous heat transfer with the surroundings through the boundary, this interface position can change with time. The interface position is unknown a priori and must be determined by solving nonlinear equations. The phase-change problems thus pose a real challenge for analysis. Phase-change problems that have been solved in the literature have often been treated with various degrees of simplification. Earliest phase-change analyses have been confined strictly to a heat diffusion study in a domain imposed with boundary conditions that are invariant with position and time. Later models account for the property variations in different phase regions and for a medium imposed with time-variant conditions. While these are more realistic physically, they also add complexity to the analysis to the extent that these problems can no longer be solved exactly. Particularly noteworthy is the recent development of methods for the solution of phase change coupled with convection and radiation. These problems can only be solved numerically. These numerical methods have recently been further extended for the solution of phase change due to a moving heat front, a problem of great importance in welding and annealing. A detailed literature review is provided in Chapter 2. Phase-change heat transfer problems can be classified according to the configuration and the formation of the interface in the medium. Thus, in some analyses, the phase change is assumed to occur at a distinct temperature which does not change with time. In such cases,

PAGE 18

3 the solid and liquid regions are separated by a sharp front whose position can be determined rigorously. This permits the interface position to be verified experimentally. In some materials, however, the phase change occurs over a narrow temperature range. Then, the phase change takes place in a mushy zone where the phase change is smoothly dispersed. Probably most intriguing of all is the phase change in alloys. In the alloys, the phase-change temperature is not distinct. There is a crystalline-like structure at the interface which is not necessarily smooth or even continuous. In practice, materials with crystalline structure and mushy zone can be analyzed more conveniently with computational approximations. This is not so, however, for a material with a distinct phase-change temperature. The interface position in such a material must be determined more rigorously. This interface position can also serve as a basis for evaluation of the accuracy of the analysis. Phase-change problems that have been solved in the literature are mostly for domains in one dimension. Solutions of two dimensional problems have rarely been attempted. Although an enthalpy method can be used to solve such problems conveniently, that method only works well if the phase change occurs in a mushy zone. It is thus the purpose of this work to develop a hybrid scheme that employs a finite-difference method in one direction and a sourceand-sink method in another direction. The former will be used in the direction where the phase change is less dominant, while the latter will be used in the other direction where the phase-change interface

PAGE 19

4 moves rapidly. The source-and-sink method is selected primarily for its convenience and accuracy as demonstrated in eight previous studies (see Hsieh, 1993; Hsieh and Choi, 1992a, 19926; Choi and Hsieh, 1992; Choi, 1991). As will be shown, the source-and-sink method only requires one set of equations to solve for temperatures in ‘all’ phase regions. By contrast, it is necessary to use a set of equations for ‘each’ phase region in the traditional method of solution by the Neumann-Stefan formulation. The work involved in the solution can thus be reduced considerably for the source-and-sink method. In addition, in the traditional method of solution, the interface Stefan condition must be given as a separate condition, whereas in the source-and-sink method, this condition has been incorporated into the governing heat diffusion equation. The source-and-sink method is thus highly efficient in the solution of the phase-change problems. In the present work, the source-and-sink method is subjected to a rigorous test for its feasibility for superposition. The method is then used to develop strategies for the solution of the problems in Cartesian, cylindrical, and spherical coordinates. For problems in Cartesian coordinates, the method has been tested analytically for its accuracy for long-time solutions. An apparatus has also been built for experimentally testing its accuracy for short-time solutions. In both efforts, the interface positions are used for comparison. Finally, the method has been applied for the analysis of two latent heat, thermal energy storage systems.

PAGE 20

5 This dissertation is divided into six chapters. A brief literature review is provided in Chapter 2. Analytical methods that are related to the solution of phase change in multiple dimensions are compiled in Chapter 3. An experimental setup is described in Chapter 4. The analysis developed in this dissertation is used in a wide range of case studies and their results are presented in Chapter 5. Finally, conclusions and recommendations are given in Chapter 6 to close this dissertation. To facilitate use of the analysis developed in this work, computer programs written in FORTRAN language are compiled in the appendices. It is hoped that the insight gained from this study will contribute to the solution of phase change in multiple dimensions in the future.

PAGE 21

CHAPTER 2 REVIEW OF LITERATURE Study of phase-change problems dates back to 1831 when Lame and Clapeyron (see Rubinstein, 1971) first studied the thickness of the solid crust formed by cooling a liquid under a constant temperature condition. Ever since, the problems have drawn attention from various quarters, and today the phase-change principles have been greatly extended and applied to study various scientific phenomena. The phasechange problems have been in studied for over one hundred and sixty years, and the phase-change literature is spread over several disciplines. It is thus a daunting task even for a brief review in this chapter. Since the present work strictly deals with the development of a two dimensional, diffusion-driven, phase-change solution scheme that combines an analytical solution in one direction and a numerical solution in another direction, only mathematical methods for diffusion-driven phase change are relevant here. It is thus the purpose of this chapter to provide a brief review of the mathematical techniques popular for the solution of phase-change problems by mechanical engineers; see a summary of methods diagramed in Figure 2.1. Interested readers are referred to books and monographs (e.g. , Crank, 1984; Yao and Prusa, 1989) for a detailed treatise on the subject. 6

PAGE 22

Exact Solutions Neumann Soluti 7 c _o h-» M T> C/3 cO 0 X -a 0 3 «T V X w TJ cr H — C/3 Js 2 0) V 2 0 3 «r S3 w 3 u X Cw T5 .2 S3 0) t . c •4-9 C 0 crj a> c 0) X c$ U H4 O X 0) 2 w S3 > £ -o 0 PQ a £ o Q >> a T8 X 4S c 3 w c o o cn T3 CA T3 IA u • «-4 3 O u *b 3 .2 >> "33 *-w J3 0) j3 c •3 P o < cn z CD c 0 T3 O X -w 2 V bC X U J = &< ca E V 3 o ICU CO s Figure 2.1 Solution methods of Stefan probl

PAGE 23

8 2.1 Analytical Solutions of the Stefan Problems 2.1.1 Exact solutions In the heat transfer literature, phase-change problems have often been referred to as Stefan problems in recognition of Stefan who first reported the solution of such problems in the open literature in 1889 (Carslaw and Jaeger, 1959). However, it has now been fully settled that the first exact solution of the problems was attributed to Neumann some thirty years earlier. NeumannÂ’s work was concealed in a series of lecture notes but was complete in the sense that it provided derivation not only of the square root time law that governs the motion of the phase-change interface as predicted earlier by Lame and Clapeyron but also of the value of the proportionality constant which was missing in Lame and ClapeyronÂ’s work. Thus, for the first time in history, Neumann solved the Stefan problems completely. NeumannÂ’s work is important for the fact that, despite innumerable efforts made by others in search of additional exact solutions after Neumann, Stefan problems that can be solved exactly today remain those that were originally solved by Neumann some one hundred and thirty years ago. Indeed, it has now been firmly established that, for a Stefan problem that can be solved exactly, the problem must be amenable to a similarity transformation. This limits the problems to be solved in an infinite domain, composed of a material of constant phase properties, and imposed with a constant temperature boundary condition.

PAGE 24

9 2.1.2 Approximate analytical solutions With the great limitation imposed by the exact solution, it is not surprising to see a wide range of approximate solutions developed and documented in the literature. One of the early approximate methods expands the interface position and the temperature in the medium in terms of a power series or error and complementary error functions. These series are then substituted into the governing equations and boundary conditions to determine the coefficients in the series expansion. Tao used this method in a series of studies. He also extended the method for study of phase change imposed with arbitrary initial and boundary conditions (Tao, 1978, 1980, 1981). Power series methods find great usefulness in the determination of short-time solutions. At short time, the series converge rapidly; only a few terms are thus sufficient for an accurate solution. However, at long time, more terms are necessary. The methods become tedious and the determination of the coefficients becomes time consuming. They lose their appeal. Phase-change problems can also be solved by using integrodif f erential equations. Lightfoot (1929) developed a moving heat source method for the solution of solidification imposed with a constant temperature boundary condition. In his method, liberation of latent heat due to solidification is taken as a moving heat source. The temperature in the medium is expressed in terms of a GreenÂ’s function. Since Lightfoot used this method for solving a problem that is possible for an exact solution, the integrodif f erential equations

PAGE 25

10 obtained in his work can be solved by a similarity transformation. This leads to results expressed in error and complementary error functions (Crank, 1984). Thus, Lightfoot’s integrodifferential solution effectively retrieves Neumann’s exact solution. Using a moving heat source front for solidification naturally leads to the use of a moving heat sink for melting. Lightfoot’s method can thus be taken as a source-and-sink method. Chuang and Szekely (1971) have used this method in the development of a numerical scheme for the solution of Stefan problems. Recently, the source-and-sink method has been greatly extended for the solution of Stefan problems imposed with temperature and flux conditions that are either constant or variable with time (Hsieh and Choi, 1992a, 1992b; Choi and Hsieh, 1992). In these efforts, the method has been shown to be highly accurate particularly for solution at long time. It has thus been used to assess the errors in some other solution methods frequently used in the literature. The method has also been applied to the solution of phase change subjected to cyclic temperature and flux conditions (Hsieh and Choi, 1992a). In these problems, both melting and freezing fronts appear simultaneously in the medium. Recently, the source-andsink method has been shown to be equivalent to the boundary element method (Hsieh, 1993). This discovery carries promise of broadening the utility of both methods for the solution of heat transfer problems. An important feature of the source-and-sink method is that the method uses only one set of equations to solve temperatures in ‘all’ phase regions. By contrast, in the conventional methods, one set of

PAGE 26

11 equations is required to solve temperatures in ‘each’ phase region. Furthermore, in the source-and-sink method, the interface flux condition has been incorporated into the governing equation. The number of equations solved in this method is thus far fewer than with the conventional methods. However, the method works best for a medium of constant and equal properties in different phase regions. For a medium with different property values, double sources and sinks or a combination of the source-and-sink and the boundary element can be used as suggested by Kolodner (1956) and Hsieh (1993). However, the equations are lengthened considerably. The method loses its attractiveness . Integrodif f erential equations have also been used by Boley (1974) in the solution of Stefan problems. Following a different approach, Boley embedded the time-variant phase regions into a large domain enclosed by a fixed boundary. The problem can then be solved in a fixed domain. However, since the conditions imposed on the fixed boundary are intended to reproduce the conditions existing on the real physical boundaries including the interface, the conditions on the fixed boundary become unknown. In one case study, Boley (1974) introduced a fictitious flux condition at the fixed boundary. This results in the solution of a system of integrodiff erential equations for the unknown flux and the interface position. Stefan problems can also be solved by a heat balance integral method (Goodman, 1958). Essentially a branch of a large class that is commonly known as the method of weighted residuals, the heat integral

PAGE 27

12 method is convenient to use particularly in the solution of phase change in a medium that is initially at the phase-change temperature (i.e., one-phase problem) and imposed with a timeinvar i ant boundary condition. In this method, the governing equation is changed to an integral form and a trial function (usually polynomials) is used to represent the temperature of the medium. This temperature function is then substituted into the integral equation to solve for the relation of the interface position with time. Since the temperature function must be assumed a priori, the accuracy of the solution is related to the appropriateness of this function in the solution. Goodman and Shea (1960) and Yuen (1981) have also applied this method to solve phasechange probems with two phase regions. It has been widely accepted that, in the heat balance integral method, using a higher-order polynomial for temperature in the trial function provides no assurance for improvement of accuracy in the solution (Goodman, 1964; Langford, 1973). However, a spatial subdivision has been found to be useful as reported by Noble (1975) and Bell (1978). Bell (1979) has further applied this method to the analysis of solidification around a cylindrical pipe by a Runge-Kutta numerical integration algorithm. 2.1.3 Other approximate solutions Stefan problems can also be solved by using quasi -steady and quasi-stationary approximations. In the quasi-steady approximation, the unsteady terms in the governing equations are dropped. This

PAGE 28

13 results in a temperature that is good for long-time solution. As for the quasi -stationary approximation, the interface velocity term in the governing equations and the interface condition are all discarded. As such, the interface is rendered stationary. Since the initial condition must still be satisfied, the temperature is good for short-time solutions . In the asymptotic method of solution, these quasi-steady and quasi -stationary limits are used for inner and outer expansions in the construction of temperature and interface position. Duda and Vrentas (1969) have demonstrated that the quasi -stationary approximation is actually the zeroth order solution of a regular perturbation, which is, indeed, good for a short-time solution. On the other hand, Weinbaum and Jiji (1977) have shown that the lowest order term of a long-time perturbation actually leads to a quasi-steady limit and is thus good for a long-time solution. These perturbation methods work well for small Stefan numbers. Since the mathematics involved in the asymptotic solution is prohibitively tedious, the method is practically useful only for first-order approximations. Solomon (1978) developed an approximate solution for a phase change imposed with a constant heat flux condition. His method follows that of the integral method; however, the energy balance is not preserved in the solution. Phase-change problems have also been solved by a Fourier transform (Ockendon, 1974). In this method, the Stefan problems are reduced to an initial-value problem. In a similar fashion, a Laplace transform can be used to change the problems into boundary-value

PAGE 29

14 problems. While the transforms are relatively easy to apply, the inversion poses a problem. The inversion process often leads to the solution of integrodif f erential equations. 2.1.4 Analytical solutions of two dimensional Stefan problems Stefan problems that have been solved in the literature are mostly one dimensional. There has been a lack of solutions for two dimensional problems. Poots (1962) extended the heat balance integral method for the solution of phase change in a square prism imposed with a uniform temperature condition. He was able to show that his analytical results were in reasonable agreement with the numerical results obtained earlier by Allen and Severn (1962). Two dimensional Stefan problems can also be solved by using an embedding method as shown by Sikarski and Boley (1965). Following a totally different approach, Rathjen and Jiji (1971) were able to combine LightfootÂ’s moving heat source method and PatelÂ’s treatment of the interface condition (Patel, 1968). They developed an analytical solution for solidification in a two dimensional corner. Budhia and Kreith (1973) followed essentially the same approach in the analysis of phase change in a planar wedge. In these efforts, the interface position must be assumed a function a priori, and the constants in this function are determined by substituting the function into integrodif f erential equations. The success of the method thus hinges on the use of the proper function in the early stage of solution.

PAGE 30

15 2.2 Numerical Solutions of the Stefan Problems It is expected that, in the numerical solution, the phasechange problems can be solved by using a finite difference method (FDM), a finite element method (FEM), and a boundary element method (BEM). It is noted that, in the study of heat diffusion without phase change, the problems are usually analyzed in a temperature domain. This is not so, however, for diffusion with phase change as in a Stefan problem. In the Stefan problem, because the phase-change interface is unknown a priori and must be determined as a part of the solution, there is an alternative to the solution in the temperature domain. In fact, the problem can be solved more conveniently in an enthalpy domain. The former (solution in temperature domain) requires a classical formulation, while the latter (solution in enthalpy domain) requires a weak formulation. They are briefly reviewed in the sections that follow. 2.2.1 Classical formulation In the classical formulation, temperature is targeted for solution. The governing equations, the initial and boundary conditions are all expressed in terms of temperature and they are used together in the solution of the interface position as well as the temperature in the medium. This method works well if the interface position is the focus of investigation. As such, the method is useful if the material has a distinct phase-change temperature.

PAGE 31

16 In the solution of the phase change by a classical formulation, it is difficult to track the interface position. Since the predicted interface position may not fall on the grid lines initially set up in the solution, iteration is usually necessary which leads to complexity in the solution. A remedy to this is to use adaptive grids or transformation of variables. In these methods, the interface position is always tied to a grid line or is immobilized in the transformed domain. Furzeland (1980) has made an extensive study for comparison of numerical methods. Some popular methods have been provided in works by Wilson et al . (1978) and Crank (1981). In the numerical solution by the FEM and FDM, the computational grids can be either fixed or time-variant. In the fixed grid method, an interpolation technique is usually necessary at each time step so that the interface positions can be determined accurately (Murray and Landis, 1959; Hastaoglu, 1986; Hastaoglu and Baah, 1991). By contrast, in the variable grid method, the grid must be modified at each time step (Gupta and Kumar, 1980,1981). Bonnerot and Jamet (1981) solved a one dimensional Stefan problem by using a third-order-accurate finite element method with an adaptive grid. In their work, the computational grid is re-constructed at the end of each time step and the results of the previous time step are mapped onto the new grid by interpolation. These methods are often referred to as front-tracking methods and have been used extensively in the numerical solution of phase-change problems .

PAGE 32

17 An alternative to the fixed grid method is to introduce an appropriate set of space coordinates so that the governing equation and boundary conditions can be transformed from the physical domain to a computational domain. In the computational domain, the interface position rests on a grid line. This method fits into the category of front-fixing or boundary-fixing. Furzeland (1980) and Lacroix and Voller (1990) have made extensive studies for comparison of the fixedgrid method and the transf ormed-gr id method. As expected, when the grid transformation is extended to the solution of phase change in two and three dimensions, the choice of the transformation coordinates, known as body-fitted coordinates, becomes particularly important. It is thus necessary that the interface front curve be smooth and continuous and the number of grid lines between the interface and the boundaries be fixed. Hsu et al . (1981) combined the transformation method with an energy-conserving finite-difference scheme in the analysis of freezing outside a cooled tube. Shyy et al . (1992) used body-fitted grids together with an adaptive grid technique in the solution of solidification in continuous casting. Additional numerical solutions by FDM and FEM have been documented in Oberkampf (1976); Saitoh (1978); Pham (1985) and Liang (1993), among others 2.2.2 Weak formulation In the weak formulation, the problems are solved in an enthalpy domain. Analogous to the Rankine-Hugoniot equation in the shock-wave analysis, the weak formulation can be used to solve phase-change problems without searching for the exact position of the interface. It

PAGE 33

18 is thus superior if the interface position is not targeted for investigation (Elliot and Ockendon 1982; Boucheron and Smith 1985). In the enthalpy method, the phase-change interface is determined a posteriori from the numerical solution carried out in a fixed domain. A survey of the recent developments of the enthalpy method has been given in Tamma and Namburu (1990) and a detailed formulation of the method can be found in Alexiades and Solomon (1993). Although the enthalpy methods have been used in the solution of phase change in many engineering problems, it is found that those methods work well if the phase change occurs over a mushy zone. For problems with a distinct phase-change temperature, the nonlinearity behavior of the enthalpy near the phase-change temperature presents great numerical difficulties which necessitate the application of an arbitrary smoothing function of the enthalpy (Goodrich, 1978) and suppression of the oscillation in the solution (Voller and Cross, 1983; Pham, 1985, 1986). While these difficult ies can be resolved by tracking the interface, it also complicates the solution. Most of all, it defeats the purpose of using the enthalpy method which is intended for bypassing the search for the interface in the first place (Yao and Prusa, 1989). 2.2.3 The boundary element method Phase-change problems have recently been solved by a boundary element method (BEM) and a combination of BEM and FDM. In the BEM, with the use of a time-dependent fundamental solution, only boundary

PAGE 34

19 grids are required in the solution. OÂ’Neill (1983) developed a noniterative boundary element analysis, which permitted the temperature gradients normal to the elements at the phase-change interface to be determined directly in the solution. Zabaras and Mukherjee (1987) developed an iterative implicit algorithm using BEM, which was shown to be free from major time-step limitations. Sadegh et al. (1987) applied BEM to the solution of solidification around a buried pipe in a semi-infinite medium with initial superheating. The planar boundary was thus initiated at the surface and propagated inward with time. The BEM has also been applied to the solution of solidification of metal casts in metal or sand molds (Hong et al. 1984). This problem has also been solved by coupling the boundary element method and the finite difference method in a conjugate heat transfer analysis (Hong et al. 1984). More recently, two and even three dimensional phase-change problems have been analyzed for laser heating, drilling, and cutting (see, for example, Basu and Date, 1990; Sarler et al. 1991; Hsieh et al. 1992). Hsieh (1993) showed that the BEM is closely related to the source-and-s i nk method and their unification is discussed in a great detai 1 .

PAGE 35

CHAPTER 3 NUMERICAL ANALYSIS In this chapter, the solution of Stefan problems by a sourceand-sink method (SSM) is presented. To ensure both accuracy and efficiency, a superposition technique with a new time marching scheme is introduced. The validity of the superposition technique for the solution of nonlinear Stefan problems is verified. Finally, the SSM with superposition technique is extended to the solution of multidimensional Stefan problems. 3.1 Problem Formulation 3.1.1 Governing equation s Consider a phase-change process that is diffusion driven in a phase change material (PCM) having constant and equal properties in different phase regions. For such a process, the governing equation can be written as ldT p (r,t) a dt V 2 T p (r,<), p = s, l , t > t 0 , r 6(1 (3.1) where the subscripts / and s denote the liquid-phase region and solidphase region, respectively. 20

PAGE 36

21 The initial and boundary conditions are given as follows: T(?,t 0 )=T 0 (r) (3.2) (3.3) where n bi denotes the outward-drawn normal of the ith boundary. For the sake of generality, the initial condition, Eq.(3.2), starts at t Q . The boundary condition, Eq.(3.3), is also general in that it reduces to the boundary condition of the first kind when kand q • are taken to be zero. It reduces to the boundary condition of the second kind when his taken to be zero. In addition to the initial and boundary conditions, the interface between the liquid and solid regions is subjected to the temperature and flux conditions. Known as interface conditions, they are given by the following relations: where n is in the normal direction of the interface. L is the latent heat of fusion (L>0 for freezing, L<0 for melting). Equation (3.5) is also known as the Stefan condition. It is this flux condition that makes the phase-change problem nonlinear. T l (S,t)=T s (S,t)=T m (3.4) dT s (S,t ) dT,(S,t) _p\L\ dS(t) (3.5) dn dn k dt

PAGE 37

22 3.1.2 Brief review of the SSM To describe the source-and-sink method, a one dimensional Stefan problem is considered. In this method, the interface condition, Eq.(3.5), is incorporated into the governing equation to form an equivalent problem (Hsieh and Choi, 1992a , 1992b) , which, in one dimensional Cartesian system, takes the following form * W = 0 + *0~ 5(<) )’ <>*0’ (3.6) The corresponding initial , boundary, and interface conditions are given as follows: t 0 ) ~ T 0 ( x ) dT(x:,t) k r~dn ’ + h il T (*i> t )T oc]=g, bl T(S,t) = T r (3.7) (3.8) (3.9) The general solution of this problem can be expressed in terms of a Green’s function as t T(x,t)= T 0 (x')G(x,t\x',t 0 )dx' £[BC ,]dr i=l + 1 } ^G(x,t\S(r),r)d7 + a T=t o t T—t (3.10)

PAGE 38

23 In Eq.(3.10), the first term on the right hand side accounts for the contribution of the initial condition, the second term accounts for the boundary conditions, and the third term accounts for the motion of the interface. The expressions for [BC] in the second term vary with the specific condition imposed at the boundary and are listed in Table 3.1 . Table 3.1 Expressions to account for the effects of boundary conditions Boundary conditions Expressions for [BC] in Eq.(3.10) T(x b ,t) = T b (t) dT(x b ,t) q b (t) dn b k b dT(x h ,t) h h i dn b + fc 6 T(Xb,< ) k h Hb ^ dG(x,t | x b ,r) T ‘M -T-G( x ’t I x b ,T) K b i G(x, 1 1 z b , t) Hsieh and Choi (1992a, 19926) and Choi and Hsieh (1992) used the SSM to solve one dimensional Stefan problems in a semi-infinite domain. In the problems they solved, the medium was initially at a uniform temperature. The boundaries were imposed with constant or time-variant temperature and flux conditions. To account for possible initial subcooling or superheating of the medium, the solutions were divided into two stages: a pre-melt stage and a phase-change stage. The solutions for these stages were then combined into a single equation as

PAGE 39

24 T(z, t) = T 0 («, t) + H{t-t 0 )^ d -^G{x, 1 1 5(r), r)dr (3.11) r=t o where t Q is the time when phase change took place, and H(t-t 0 ) is the Heaviside function given by 0 1 for t < t Q for t >t 0 Since the initial temperature was uniform, the first term on the right hand side of (3.11) could be expressed in a compact form as r 0 (x, t ) t [af e { t ) V t T ~ l o e 4 “(‘T W (3.12) where xT b (r) 2a t-T ^( r ) for the first-kind boundary condition for the second-kind boundary condition In practice, since the interface position is unknown a priori, Eq. (3.11) must be solved numerically. In this effort, the entire time range is divided into small time intervals, Ai, over which numerical integration is performed as N T-t F(r)d T = n— 1 .n-1 F(t)cIt (3.13)

PAGE 40

25 where = t and t n = t 0 + (n — l)At . One of the advantages of this numerical method is that its accuracy does not degrade with time. Large errors only occur over the first few time steps and they diminish rapidly with time. Accurate results can thus be obtained with this method. However, since the time integration always traces back to the initial condition at t Q , the effort of computation may be considerable, particularly at long time. 3.1.3 Time marching metho d Another approach to the solution of Eq.(3.11) is the use of a time marching scheme. In this approach, the solution is expressed as T(x,t n+1 ) = | T(x',t n )G(x,t n+1 1 x',t n )dx' ,n+l + # E ( T ) 4 a(t n+1 -r) „ (t n+1 ~r) 1/2 + t\ ^G(x,t n+ 1 \S(rlr)dr dr T-t" n+1 (3.14) The time integration is thus evaluated one step at a time, and the temperature evaluated at each time step serves as the initial condition for the succeeding time step and it moves on with time. For this method, a domain integral emerges as the first term on the right. In this term, the temperature at time t n is T(x,f n ) which is nonuniform. A numerical integration must therefore be performed, giving rise to additional computational effort. It also degrades the accuracy, which

PAGE 41

26 is particularly pronounced when the domain for integration extends to infinity . 3.2 SSM With Superposition Te chniques To alleviate the problems stated in the preceding section, a superposition method with a new time marching scheme is developed. Before addressing this method, the validity of the superposition method for the solution of nonlinear Stefan problems is critically analyzed . 3.2.1 The validity of the superposition method The superposition method has been widely used in the solution of linear boundary value problems. Although this method has also been applied to the solution of nonlinear Stefan problems, its validity has never been established. For the ease of derivation in what follows, it is assumed that the domain is one dimensional. The PCM has constant and equal property values in different phase regions. Thus, the governing equations are all linear. The initial and boundary conditions are also linear if the radiation effect is incorporated into the convection effect by means of an overall heat transfer coefficient. The only nonlinear relation comes from the interface flux condition. It is thus necessary to decompose this condition if the superposition method is to be applied to the solution of the Stefan problem. To start the derivation, a small constant e is introduced such that, for a melting problem,

PAGE 42

27 T(S-e, t)=T,(S, t), T(S+€,t)=T s (S,t) As e approaches zero, lim T(S±e,t) = T m (3.15) By definition, the differential of T at x equal to S can be recast as dT dt x-S dT(S ± e, t) dx dx dT(S±e,t) ^ x-S dt (3.16) where, on the right hand side, the partial differentials are to be evaluated in the same phase region. That is to say, a S+e in dT/dx is to be followed by a S+e in dT/dt. Then dS/dt can be expressed as dS_ dt dT(S±e,t)/dt dT(S ± e ,t)/dx e->o (3.17) Substituting (3.17) into (3.5) changes the interface condition to the f ol lowing dT s dT, p\L\ dT(S ±e,t)/dt _ dx dx x=s~ k dT(S±e,t)/dx (3.18) which is clearly nonlinear. An alternative form of Eq.(3.17) can be derived by writing it in both solid and liquid regions and equating their results as dS _ dT(S+c,t)/dt dt dT(S+e,t)/dx dT(S-e,t)/dt dT(S-€,t)/dx e-^o e— >o (3.19)

PAGE 43

28 which can be recast as dT{S-e,t) fdT(S-e,t) '<9T(5+£, t)/0x ) dx 5 1 o T dT(S+e,t)/dt 1 Since for the Stefan problem (3.20) dT(S+e,t) dT(S-e,t ) ] dx dx -l-o ^ (3.21) because of the flux jump at the interface, one can multiply the numerator and denominator of Eq.(3.17) by the partial differential difference appearing on the left hand side of (3.21) as dS _ dT(S+e,t)/dt dT(S+e,t)/dx-dT(S-e,t)/dx dt [, dT(S+e,t)/8x ' dTls+e,t)/dx-dT(S-e,t)/dx (3.22) Then, expanding the product in the numerator of the fractions in this equation and introducing (3.20) and canceling out terms enable the following relation to be derived dS_ dt dT(S+e,t)/dt-dT(S-e,t)/dt dT(S+c,t)/dx-dT(S-e,t)/dx e-»o (3.23) Equation (3.18) can then be expressed as ’dT s dT, P\L | dT(S+€,t)/dt-dT(S-€,t)/dt ' dx dx 1 CO II H dT(S+e,t)/dx-dTis-e,t)/dx These equations will now be used to validate the (3.24) use of the superposition principle in the solution. In this effort, the

PAGE 44

29 temperature is decomposed into a linear part U and a nonlinear part V such that T(x, t) — U(x, t) + V{x, t) . For the diffusi on problem solved, U is a smooth function of x and t. It thus has continuous first derivatives even at the interface, given by the relation dU(S +e,t ) 8U(S-c,t) _5U(S+e,t) dU(S-c,t ) dx dx dt dt ~ as c — >0 (3.25) Equations (3.23) and (3.24) can then be reduced to dS_ dt dV(S+e,t)/dt-dV(S-e,t)/dt dV(S+e,t)/dx-dV(S-e,t)/dx e— o (3.26) dV, _P\L\ dS(t) dx dx x=S k dt (3.27) Equations (3.26) and (3.27) show clearly that U can be totally removed from the interface flux condition; only the V problem has a nonlinear condition at the interface. The key point in the superposition method is that the [/-problem is linear and smooth at the interface, while the V problem is nonlinear which contains all the features characteristic of the Stefan problem. It is noted that the derivation above can be easily extended to the solution of multi-dimensional Stefan problems. The derivation also illustrates well that the superposition principle is only valid for situations when the PCM has constant and equal properties in different regions .

PAGE 45

30 3.2.2 Solution method To describe the solution procedure, a general one dimensional Stefan problem is taken for analysis. For this problem, the governing equation is given by the relation LdT = ±dJ ddT ) t t o n oon « dt x d tik dx)’ t>l ° ’ ( 3 . 28 ) where x d is the Sturm-Liouvi 1 le weighting function, given as follows: 0 slab d = < 1 cylinder 2 sphere The initial, boundary, and interface temperature and flux conditions are, respectively, 1? o hi II o hT (3.29) d T (x-.t) K Qn! + h il T (*i,t) T 00 ]_ 9i bl (3.30) T(S,t) = T m (3.31) dT s dT, 1 P \L\ dS(t) dx dx k dt x—S (3.32) Application of the superposition as described in the preceding section leads to the solution of two problems in which the linear problem U(x,t) satisfies the following equations

PAGE 46

31 I QU_± dJ r ddU\ a dt~ x d M dxr t> t 0 , (3.33) U(*,t 0 ) = T 0 (x) (3.34) k qfc + h iWiT oo) = 9i (3.35) In this problem, Eq.(3.25) is satisfied implicitly because of the smoothness of U. Correspondingly, the V(x,t) problem satisfies the following relations 1 dV _ J_ djdd•.\ a dt~ x d M dxr t> t 0 , x e n V(x,t 0 ) = 0 + h,V, = 0 Vi(S,t) = V.(S,t)=T m -U(S,t) _ p \L\ dS(x ) dx dx ~ k dt -‘X-S (3.36) (3.37) (3.38) (3.39) (3.40) The {/-problem, being linear, can be solved exactly by established methods, such as separation of variables or integral transformations. However, these analytical solutions often lead to

PAGE 47

32 infinite series of trigonometric or Bessel functions, which are slow to converge. In addition, the eigenvalues must often be found numerically. For these reasons, the analytical solutions may be impractical. Use will therefore be made of a numerical solution by a finite difference method (FDM) for the solution of the {/-problem. To solve the U-problem, the source-and-sink method is employed. Following Hsieh and Choi (1992a), the interface condition (3.40) is incorporated into the governing equation (3.36) to form an equivalent problem as follows I dV_ 1 d(ddV\ , LdS a dt~ x d dA x dx) +ac dt &( x ~ S(t)) , t>t 0 , x e Q (3.41) This equation can be solved in terms of the GreenÂ’s function. In this effort, since the initial and boundary conditions are all homogeneous, (see Eqs.(3.37) and (3.38)), the V problem can be solved as t V(x,t)=^ ^ p-S d ( T )G(x,t\S(T),T)dT (3.42) r=t o Application of the interface condition, Eq.(3.39), gives a relation T m ~ U(S,t) = 4 dS(r) dr S d (r)G(S(t ) , { | S(t) ,r)dr (3.43) which is to be solved implicitly for the interface position. In this effort, it is noted that, if S is located at a grid point, {7(5, < ) can be d irectly obtained from the finite difference solution. Otherwise, a cubic spline will be used to determine U at point 5.

PAGE 48

33 The solution procedure can now be stated as follows. For each time step after phase change, the T-problem is decomposed into a Uproblem and a V'-problem. The {/-problem satisfies the nonhomogeneous initial and boundary conditions (see Eqs.(3.33) through (3.35)), while the ^-problem satisfies homogeneous initial and boundary conditions together with interface conditions (see Eqs.(3.36) through (3.40)). In the solution, it is convenient to use (3.41) as the governing equation for the F-problem, which incorporates (3.36) and (3.40), a reduction of one equation. Then the {/-problem is solved first by FDM . The Vproblem is solved next by SSM. These solutions are then superimposed and this is done for each time step. In this method, the temperature evaluated at the new time level will be used as the initial condition for the (/-problem at the next time level and it moves on incrementally with time until the desired time is reached. A summary of the procedure can be given as follows: (1) Start from an initial condition T Q (x) and set n = l, (2) Set t 0 =t n ,t = ( n+1 and solve Eq.(3.33) for U(x,t n+1 ), (3) Solve Eq . (3 . 43) for S(t n+1 ), (4) Calculate V(x,t n+l ) from Eq.(3.42), (5) Find T{x,t n+1 ) = U(x,t n+1 ) + V{x,t n+1 ), (6) Set T 0 (x) = T(x, < n+1 ) , n = n + 1 , and repeat (2) to (6). An important feature in this solution method is that the integral in V (x,t) (see Eq.(3.42)) is evaluated one step at a time and

PAGE 49

34 the F-problem always starts with homogeneous initial and boundary conditions. No domain integration is thus involved. Since a time marching scheme is employed, the method is equally efficient and accurate for short time as well as for long time. In the numerical evaluation of the integrals in Eqs.(3.42) and (3.43), a multi-point Gauss-Quadrature formula is used. With sufficiently small time steps, the interface velocity can be approximated as dS(r) _ S(t n+1 )-S{t n ) dr ~ At (3.44) where t n < r < < n+1 , and A t = t n+1 — t n . Then, only one unknown 5(<” +1 ) appears in Eq.(3.43). This unknown can be found by using a rootfinding method. To ensure a fast and convergent solution of this method, a conjugate Newton-Raphson and bisection method is employed. It is expected that, for problems imposed with a sudden change of temperature at the boundary, the interface velocity ^ tends to be infinity at time equal to t Q . A large error may thus result if a finite and constant velocity is assumed in the first time step (see also Hsieh and Choi, 1992a). To reduce this error, a slow-start procedure is used in which the time step is treated as a variable. It starts with a very small value and increases gradually with time until it reaches a preset upper bound. In this way, the errors can be reduced considerably .

PAGE 50

35 3.3 Solutions of One dimensional Stefan Problems The analysis given in the preceding section will now be used to solve one dimensional Stefan-problem examples in Cartesian, cylindrical, and spherical systems. For the examples in Cartesian system, both semi-infinite and finite domain (plan slab) will be solved, of the cylinders and spheres will be considered. For these examples, the governing equation, initial, boundary, and interface conditions have been given as Eqs.(3.28) through (3.32). The U and V problems are formulated as Eqs.(3.33) through (3.40). The V problem has been solved using the equivalent problem (3.41) with the result given as (3.42). The interface position is to be solved using (3.43). Solutions for these examples are now compiled in the subsections that follow. 3.3.1 Solutions in one dimensional Cartesian s ystem U -problem: The governing equation for the (/-problem in a Cartesian system can be obtained from Eq.(3.33) by setting d equal to zero. To ensure numerical stability for all time steps, an implicit finite difference method is used with the results given as follows: whereas in the cylindrical and spherical systems, only finite domain (3.45a) or (3.45b) where 9 x — aAt/Ax 2 . Here, the superscript n refers to the time level and

PAGE 51

36 the subscript j refers to the nodal point location. Thus, U'} = U(x,t n ), —U(x± Ax, < n+1 ) Notice that Eq. (3.45a) is second-order accurate in space and firstorder accurate in time and is often referred to as the simple implicit method. Equation (3.45b), however, is second-order accurate in both space and time and is commonly known as the Crank-Nicolson method. Both methods have been well established and are unconditionally stable for all time steps. These equations can be solved by using a tridiagonal system of linear algebraic equations with a simple forwardand backward-substitution, an efficient numerical procedure. V'-problem: The solution of the V'-problem can be directly obtained from Eq . (3 .42) by setting d to zero. For simplicity sake, let t n be t 0 and t be . Then, at the end of the first time step, V(x,t n+1 )=V(z,t 1 )=i -^ 7 ^
PAGE 52

37 A close examination of the Green’s function in (3.47) reveals that a singularity is located at r = < 1 . This singularity may degrade the accuracy of the numerical integration and can be resolved by introducing a transformation r/ = ^j4a(< 1 -r) (3.48) so that Eq.(3.46) is changed to >j4aAt -(x-sy (x+sy L dS f i 2 2 e n ± e 77 c dt 2Wa *7—0 dr] (3.49) For the PCM in a finite domain, such as a plan slab (0
PAGE 53

38 oo 2 / G( x ,t\x',r) = C 0 + e~ af3 '' {t T) cos(p n x)cos((3 n x l ) (3.51) n = 1 where , for a temperature condition at x = H, for a flux condition at x = H. for a temperature condition at x = H , for a flux condition at x — H. In practice, when the Green’s functions are substituted into Eq. (3 .46) , a tedious and time consuming algorithm will be necessary to evaluate the integral because of the repetitive computation of the infinite series. However, by assuming a constant interface velocity from t 0 to ij , Eq.(3.46) can be simplified to the following forms (see Appendix A for derivation): v (^ t i)=j^f^B n (t 1 )sin(/3 n x) (3.52) n—l (2n-l)7r H ’ and / °’ H , H ’ for a temperature condition imposed at x = 0, and V(*,t 1 ) = C 0 £ ^At + j^f^CMcosi^x) n= 1 (3.53)

PAGE 54

39 for a flux condition imposed at r = 0. In these equations, S n(*l) = t) ) 2 |( a ^n) s »' n (/ ? n 5 l) “ (^„w)cos(/?„5i) -e a0nAt [(or Pl)sin(/3 n S 0 )-(p n v)cos((l n S 0 j\ | (3.54) and C n(
PAGE 55

40 V'-problem: The solution of the V'-problem can be obtained by setting d equal to unity in Eq.(3.42). Then, *1 V ^h)=k\ ^r-S(T)G(x,t 1 \S(r),T)dr (3.58) T-t O The Green’s function for a solid cylinder (0 < x < R c ) with a temperature condition imposed on the boundary at R c is given by the relation C(*,l I |*',r) = i (3.59) H c n=l J \\Pn H c) where J Q and J x are Bessel functions of the first kind of order zero and one, respectively. The eignvalues (3 n are roots of J o ((1 n R) = 0, n= 1,2,3... (3.60) 3.3.3 Solutions in spherical coordinate system U -problem: The governing equation for the {/-problem in a spherical system can be obtained from Eq.(3.33) by setting d equal to two. An implicit finite difference solution for such a problem can be derived as -(0 x -29' x )U^l + (l + 2e x )U^ +1 -(9 x+ 20' x )U^l = U t ;, if x^O (3.61) and ( 1 + 60 x )t/" +1 — — U" , if x = 0 (3.62)

PAGE 56

41 where O x and 0' x have been defined previously. These equations can also be solved by using the tridiagonal solver. K-problem: For a problem in spherical system, the solution of V can be obtained directly from Eq.(3.42) by setting d equal to two. This gives V(xJ i)=-g(r)G(x,t l \S(T),T)dr (3.63) O The GreenÂ’s function for a solid sphere (0
PAGE 57

42 The solution for the ^-problem can also be simplified to the foil owing form (see Appendix B for derivation): v(^i)=^E^w^ sin ^) n=l (3.68) where ^n(*l) — /. o2\2 {<*&)* + {fi n v) ^j“4 5 1 Sm( / 3 " 5 1 ) S o e~ a0nAt sin(P n S o ) ~ (Pn v ) S,co,^S t )V'" 3 '^ »»(£),, Sj a/3^At 2aP^v cos{P n S 1 )-e n cos(P n S 0 ) (3.69) 3.3.4 Multiple interfaces The method developed above can also be applied to the solution of problems involving multiple interfaces. In this case, the solution of the {/-problem remains unchanged. However, the solution of the Vproblem can be obtained by using multiple sources and sinks. The equivalent governing equation for this problem is 1 dV__ 1 df r ddV\, V' L dS P a dt~ x d JA dxr 2^ ac dt x P = l «(*-5 p (0) (3.70) where M is the number of interfaces. The general solution of this equation is

PAGE 58

43 = Y 1 ds »(r) li\ = 1 1 , T—t dr Si(r)G(x,t 1 \S Jr),r)dT (3.71) Here, the interface positions are determined by solving the following nonlinear equations simultaneously T m ~U(S v t 1 )= ££[ d ^ls d p (r)G(S 1 ,t 1 \S p (T),T)dr p=1 A M r dS (t) T n -U(Svt l )= Y,7\ Sp(r)G(S 2 , t x | 5 p (r),r)dr P = 1 T=t (3.72) Tm-U{S M ,t l )=Y j ± dSJr) t P=1 r=t dr S d
PAGE 59

44 Following the Von Neumann stability analysis, a small disturbance of T is defined at time level n so that T n = T n + e n (3.73) where the hat denotes the exact solution and e is the disturbance. In the superposition method, the temperature at a new time level n+1 can be decomposed as T n+1 -u n+ 1 +V n+l (3.74) Then, the disturbance at the new time level becomes e n+l =( n+l +e n+l (3.75) where £j comes from U and e 2 comes from V. As for the present work, the solution of U is obtained by using a finite difference method. Thus , n+1 _ e T Cl — 1 l+20( 1-cosA) (3.76a) if the simple implicit method (3.45a) is used, and n+1 l-g(l-cosA) „ 1 l+0(l-cosA) (3.76b) if the Crank-Nicolson method (3.45b) is used (Anderson et al. 1984). Here, in both equations, A is the phase angle.

PAGE 60

45 As for V and 5, Eq.(3.41) must be first linearized so that the Dirac delta function 6(x S) can be treated as constant. Then, if V and 5 are represented as V = V + e 2 , 5 = 5 + 7 (3.77) The disturbance e 2 must satisfy the following equation i ^2 “ dt (3.78) which is of the same form as (3.41). The disturbance e 2 can then be expressed in terms of 7 , the deviation in the interface position, as 69 — ' ^G(x,t 1 \S(T),T)dT T=t „ (3.79) As discussed previously, for small At, dj/dr can be treated as constant. It can thus be taken out of the integral as 1 727 j G{x,t x \S(T),T)dl T—t (3.80) It is noted that, in this work, the V problem has been solved by using a homogeneous initial condition, the disturbance in V thus will not accumulate. However, the disturbance in U can be transferred to V through 5 which is solved using Eq.(3.43).

PAGE 61

46 It is clear that, at the interface, V(S,t l ) = T m -U(S,t 1 ) Hence , “ £ l Introduction of (3.82) into (3.80) yields ,fi |5(r),r)dr or dy dt L c G(S,t 1 \S(T),r)dT So that Eq.(3.80) can be written as f 2 = n fl II is defined as [ G(x,t 1 \S(T),r)dT G(S ,t 1 \S(r),T)dr The disturbance at the new time level can then be expressed as c n +1 _ (1-n) e n [1+20(1 — cosA)] (3.81) (3.82) (3.83) (3.84) (3.85) (3.86) (3.87a)

PAGE 62

47 for the simple implicit method, and € n+1 _ l+0( 1 -cosA )' 1 (3.87b) for the Crank-Nicolson method. Finally, the amplification factor, defined as the ratio of t n+ l and e n , can be related to II as and f n + 1 (i-n) e n [1+20(1 cosA)] (3.88a) A l-0(l-cosA) l+0( 1-cosA) ( i — n) (3.88b) for Eqs. (3.45a) and (3.45b), respectively. Since 1
PAGE 63

48 Position, x/H Figure 3.1 II versus x plot

PAGE 64

49 concluded that the amplification A can never be greater than one, providing assurance of an unconditionally stable solution. When a numerical scheme is developed for the solution of the Stefan problem, it is necessary that the numerical solution converges to the true Solution. For numerical solutions using finite difference methods, convergence analysis is usually accomplished through the application of the fundamental equivalence theorem of Lax. According to this theorem, a numerical scheme is considered to be convergent if and only if the scheme is stable and consistent. In the superposition method, the {/-problem is solved by a finite difference method, while as the ^-problem is solved by an analytical method. The numerical scheme will be stable and consistent If both the {/and the K-problems are stable and consistent. It has been well established that the finite-difference equations, (3.45), for the {/-problem are consistent with their differential-equation counterpart in (3.33). The solution of the F-problem is exact up to Eq.(3.43). In (3.44), the velocity of the interface is approximated as a finite-difference. It is obvious that this equation is also consistent. Thus the entire solution of the T problem is consistent with the original governing equations. As it is also stable for all time steps, convergence is thus assured. It should be noted that the above analyses have been carried out on the basis of linearized equations. The nonlinearity effects remain to be investigated. Use will thus be made of numerical experiment later to test the stability (Chapter 5).

PAGE 65

50 3.4 Extension to Multi-Dimensional Stefan Problems The analyses developed above Tor the solution oT one dimensional Stefan problem will now be extended to the solution of multidimensional Stefan problems. Again, the phase change is considered to be diffusion driven in a PCM having constant and equal properties in different phase regions. 3.4.1 Solutions in two dimensional Cartesian systems For phase change in a two dimensional Cartesian system, the governing equation, initial, and boundary conditions are, respectively , ldT_d^T,d 2 T ... . a dt ~ dx 2 + dy 2 ' t>l o' ( x >2/)£ft (3.89) T(x,y,to) = T (3.90) dT ^^+\-[T 6t .-T 0O ] =9f. * = 1 >2, . . • (3.91) Following Patel (1968), the interface position (Figure 3.2) can be expressed in the y direction as S = S(x, t) Then the interface temperature and flux conditions become (3.92) T l (x,S,t) = T s (x,S,t)=T„ [i.fasyi dT s (x, S ,t) dT,(x,S,t) L + \dx) \ dy dy p | L | dS(x , t) ~k Ft (3.93) (3.94) The governing equation for the equivalent problem becomes

PAGE 66

51 Figure 3.2 System analyzed in a two-dimensional Cartesian coordinate system

PAGE 67

52 1 dT _ d 2 T , d 2 T , L dy a dt ~ dx 2 dy 2 ac dt 6(y-S(x,t)) Application of the superposition method T(x,y,t ) = U(x,y,t ) +V(x,y,t) (3.95) (3.96) leads to the solution of two problems, for which the {/-problem satisfies , Q T, *2t T q2t, (3.97) W.#u + £u t t>toi {x>y)en a dt Q X 1 dyl U(x,y,t 0 ) = T 0 (x,y) (3.98) k i dnl] + T °°~^ ~ q i ’ i=l ,2, (3.99) For the ^-problem, the governing equation is IdV _ d 2 V , d 2 V . L dy c ,_. , x c0 a dt ~ dx 2 + dy 2 + ac dt ° ^ d(x,t)), t>t 0 , (i,y)GS2 (3.100) The initial, boundary, and interface conditions for the ^-problem are V(x,y,t 0 ) = 0 (3.101) hs£ + *. y « = 0 ’ ‘ =1 > 2 — (3.102) V l (x,S,t)=V s (x,S,t)=T m -U(x,S,t) (3.103)

PAGE 68

53 As discussed previously in Section 3.2.2, the £/-problem can be solved by FDM. For a two dimensional problem, it can be solved efficiently by an alternating direction implicit (ADI) method which has been established to be unconditionally stable for all time steps. In this method, the temperature in xand y-directions are calculated alternately. The finite difference equations are vr+ ,Y?+ (l-wpti* 3 + 1 VW-i + V.'+.V, = (3.104) where 0,,. ~ 2Ax* an< ^ ^ y ~ 2A ^ 2 " These finite difference equations can be solved by using two tridiagonal systems of linear algebraic equations by means of an efficient tridiagonal solver. The solution of the ^-problem can be expressed again in terms of a Green’s function as V(*,y,* 1 ) = £[ dr f &G(x,y,t 1 \x',S(x',T),r)dx' (3.105) r=t 0 x '—0 The interface position curve 5(x,ij) is then determined implicitly by the foil owing equation H -U{x,S,t 1 ) = ± T — t dr dS | x',5(x',r),r )dx' x —0 (3.106)

PAGE 69

54 In this effort, since the interface S is a function of both x and t, Eq. (3.106) can not be solved directly for S. The difficulty lies in the fact that a two dimensional searching for the position is necessary if this position is to be found iteratively. Alternatively, an assumption of the S function to represent the interface position will be necessary if this position is to be found analytically. Moreover, in the analytical solution, the accuracy of the result depends greatly on the fitness of the function that is initially assumed for approximation. In any event, the interface position can not be determined easily. To overcome these problems, a new approach is developed by splitting the two dimensional problem into an array of transversely related one dimensional problems. In this approach, the second partial differential of V with x in (3.100) is changed to a finite difference as d 2 V(x,y,t) _ K(z+Az, y, t) 2V(x, y, t) + V(x-Ax, y, t) dx* Ax* (3.107) Then, by indexing x position with subscript i so that V(x,y,t) is changed to Vi(y,t) and V(x±Ax,y,t) to V,± j ( y, t ) , Eq.(3.100) can be re-written as follows ldV,(y,t) _d 2 V,(y,t) , V i+1 (y,t)-2V i (y,t) + V i _ 1 (y,t) a dt ~ Qyt + A -v.2 Ax* +4 % Ky-Si(t)), i = 1,2,... (3.108)

PAGE 70

55 Equation (3.108) represents a group of piecewisely related one dimensional problems which can be solved by means of the Green’s function using a time marching scheme. In this effort, the last two terms on the right hand side of Eq.(3.108) can be collectively taken as source terms. Then with the help of zero initial condition and homogeneous boundary conditions imposed in the V'-problem, the solution of V can be written as where ^t(y^l) — c -foT G i(y’ t il S i’ T ) dr r=t. l l D + a dr 9i(y' ,r)G t (y,t 1 \y l ,r)dy' r =*o y'=0 (3.109) 9i(y',r) y, + l(y',r)-2V i (y',r)+V^ 1 (y',T ) Ax 2 (3.110) It is still not possible to solve Eq.(3.110) because of the function g^(y,t) which depends on l^ ( (j/,i). Some additional work is necessary, and one way to do it is to expand g about t 0 as 9i(y, t) = <7,(2/ , t 0 ) + g,( ^ o) (r t 0 ) + . . . , (3.111) If higher order terms are neglected, then g can be expressed as 9,{y, r) = g,(y , t 0 ) + <0 \ r 1 0 ) (3.112)

PAGE 71

56 Using Eqs. (3.101) and (3.110) permits the first term on the right hand side of this equation to drop out. As for the second term, it can be shown from (3.110) that dg { (y,t 0 ) 2 dt Ax 2 dV i+l (y,t Q ) o 0V,.(y,i o ) dV { _,(y, t Q )] dt z Ft + Ft (3.113) Then with the help of the initial condition given as (3.101), Eq.(3.100) can be used to re-write (3.113) as d 9,(y,t 0 ) dt 1 2®^5 .+i) +§f%-S i 1 ) cAx (3.114) 1 t=t„ It thus follows that V i(y> t G i(y^i\ S i’ T ) dr + ^2 ( u ’i+i 2u '< + «’i1 ) (3.115) where ( T ~ t o)— §f^ G i(y> t il s K ( t o)’ T ) dT ’ k = i, i ± 1 r=t_ (3.116) In practice, a first-order approximation is sufficiently accurate to solve the problem at hand. For the first-order approx i mat i on g i(y',r)~ gi (y', t 0 ) (3.117)

PAGE 72

57 Then, from (3.101), r) = 0 (3.118) and V can be obtained by solving a simple equation as g s v i(v,t i)=|J -q? G,(y,< 1 |5,,r)dr (3.119) r=t o 3.4.2 Solutions in cylindrical coordinate systems The governing equation for heat diffusion in a two dimensional cylindrical system can be written as follows = 0 + l>t o ’ 0 4)en (3.120) As is always the case, appropriate initial and boundary conditions as well as the interface conditions must be specified. For the phase change in cylindrical coordinates, two types of the interface motion will be addressed — the phase-change interface may be moving in zdirection (Figure 3.3), and the phase-change interface may be moving in r-direction (Figure 3.4). In the former, the interface position may be expressed as S = S(r,t) (3.121)

PAGE 73

58 Figure 3.3 System analyzed in cylindrical coordinates Interface moves in ^-direction Figure 3.4 System analyzed in cylindrical coordinates Interface moves in r-direction

PAGE 74

59 Then, the interface conditions are written T(S(r,t),r,t) = T„ (3.122) Mf) 2 ] dT s dT, dz dz z-S _ P\L | dS(r,t) k dt (3.123) In the latter , the interface position may be expressed as S = S{z,t) (3.124) Then, the interface conditions are given by T(z,S(z,t),t) = T Mi ) 2 dT s dT z dr dr J r=S _ P\L\ dS(z,t) k dt (3.125) (3.126) In the method of superposition, the temperature is decomposed as T(z, r,t) = U (z, r,t) + V (z, r, t) (3.127) Here, the [/-problem satisfies the same governing equation for T as in Eq.(3.120), and the nonhomogeneous boundary and initial conditions. This problem can be solved by a finite difference method as Kl U li = 2 W + 1, i 2 U if j + Ui_ hj ) + W lJ+ ,-P,, r i) + 2« r (f |)J+1 -2f I . J + f | . rl ) (3.128) where 0 = -9l4^ 0 ’ = a nd a-aAt 2 2A z 2 r 2rAr r 2 A r 2

PAGE 75

60 Equation (3.128) can be solved either by an explicit method by setting U = U n on the right hand side, or by an implicit method by setting U = U n+1 on the same side. In any event, the solution is routine. The ^-problem can be solved by the Green’s function method. Again, the interface flux condition is incorporated into the governing equation; Then, for the case of interface moving in the 2 -direction, the governing equation becomes IdV _ d 2 V . ld( r dV\ , L dz 01 dt d z 2 r dr\ dr)' ac dt 6(z-S(r,t)) (3.129) The governing equation is split into two directions, where a finite difference method is used in the r-direction. This reduces the V equation to the following a dt d 2 V , dz T + h §f 6(*-SjV)) + gj ( Z,t) (3.130) where » y ( 2 . 0 = 2rAr ^ r 2 4 (^ + i Ar 2 ’ for r ^ 0 (3.131) for r = 0 As described in Section 3.4.1, the first-order approximation of g leads to a very simple approximation as g j (z,T)~g j (z,t 0 ) = 0 (3.132)

PAGE 76

61 The F-problem can be solved as t f dSV j( z >h)=?\ G J (z,t 1 |5 J ,r)d7 r=t(3.133) Similarly, the second-order approximation of gresults in the following solution where 1 ) — c Gj( z , < 1 1 5^ ,r)dr + Wj 55 , (3.134) T — t .. W i = w j + 1 *Vi , ^j+1 2u >j + Ar 2 2rAr 4 K+l «>>) Ar 2 ’ for r / 0 for r = 0 (3.135) and tn. i = | ( r < 0 ) -Qt° ) G j (z,h\S K (t 0 ),T)dT, K = j,j± 1 (3.136) T = f _ For the case of the interface moving in the r-direction, the governing equation of the equivalent problem is given by IdV _ d 2 V , 1 d (_dV\ , f dr C/ a 5< 5 Z 2 + 5r J + QC 5< (3.137) Again, the governing equation is split into two directions; this time, a finite difference method is used in the z-direction. The V equation becomes

PAGE 77

62 1^,M) id a dt ~ ^ dV t (r,t) dr + fj K r ~ S >,(<)) +ff,(M) (3.138) where j( /) _ y .+ l(»-.0-2y,-(r > 0+V 1 -1 (r > 0 ' ’ A z 2 (3.139) The first-order accurate solution of the V-problem is readily obtained as l l ^(^i) = £ ^5.(r)G t (r,t 1 |5,.,r)dr T—t r ds and the second-order accurate solution can be derived (3.140) ^>(Mi)=£j ^5,(r)G | .(r,l 1 |5 | .,r)d7 + ~£? (“i+l 2w i + ”i-l) where t W K = T — t (3.141) * ^ ^ ^ \ k = i, z'±l (3.142) 3.4.3 Solutions in three dimensional Cartesian systems The governing equation of a Stefan problem in a three dimensional Cartesian system takes the form

PAGE 78

63 IdT _ d 2 T ,d 2 T , <9 2 T a dt ~ dx 2 + dy 2 dz 2 ’ t > t 0 , (x,y,z) £ ft (3.143) The interface position is expressed as S S(x,y,t) Then, the interface condition can be written as (Patel, 1968) With the use of the superposition method, the temperature is decomposed as _P\L | dS dt (3.144) T{x,y,t) = U(x,y,z,t) + V(x,y,z,t) (3.145) in which the [/-problem satisfies a governing equation that is of the same form as (3.143) and relevant nonhomogeneous initial and boundary conditions. The finite difference equation for the {/-problem can be derived as ^ i,j,k ^x^i+l,j,k 2Ui,j,k + Ui-l,j,k) + i, k + u it (3.146)

PAGE 79

64 where 6 X , 6 y , and O z have been defined previously. This equation can be solved by using an iteration procedure such as Gauss-Seidel and SOR method, or by a modified ADI method. Again, the F-problem can be solved by a Green’s function method. In this method, the equivalent problem can be written as IdV _ d 2 V d 2 V a Q x 2 dy 2 i d V . L dz c / + dz 2 +ac dt ^ S ' (3.147) Then, by using a finite difference method in both xand y-directions , Eq. (3.147) can be re-written as l^.i a dt d 2 V. dz T^ + 53 |f S(z-S) +g ifj (z,t) (3.148) where 9i, ]{*>*) = Az 2 + Ay 2 (3.149) With the help of the homogeneous initial condition for the F-problem, the first-order accurate approximation of g { leads to a very simple relation as 9i, j(*> t) ~ g fa t 0 ) = 0 (3.150) and the solution of the ^-problem follows as T f dS i i V i,i(*>h)=T] ST dr T — t^ (3.151)

PAGE 80

65 In a similar fashion (see Section 3.4.1), the second-order accurate approximation of g { leads to a solution as r 95 = ~c q t G i 'j( z ,t 1 \S^j, T )dT T—t~ + 7a? where +JK? + (3.152) w. i,n | ( r t o ) ~qI G i ^j( z ,t l 1 5 m n (< 0 ),r)dr m = i, i± 1, n = j,j± 1 (3.153) T — t „ 3.4.4 Stability analysis For multi-dimensional problems, the explicit treatment of the second partial differential may impose additional constraints on the stability of the solution. Due to the difficulties associated with the nonlinearity of the problem, the stability can not be tested by simple analysis. However, some comments can be made as follows. It is well known that the solution of U is unconditionally stable. As for th e V , since it always starts with a homogeneous initial condition, the disturbance (error) of V will not accumulate .

PAGE 81

66 The only source of error that may be brought into the ^-problem is U(S,t) which is used to find the interface position. In this connection, because of the use of the GreenÂ’s functions, disturbance in the integrand will not be magnified but will be suppressed. On the other hand, the disturbance in the V solution can still be brought to the t/-problem through the superposition principle. However, this disturbance can still be suppressed by the ADI method in the solution for the temperature in the succeeding time step. It is thus expected that any disturbance will not grow in the solution, and, as a result, the method is unconditionally stable. It is noted that the discussion above is not meant to be rigorous. Numerical experiments for a twodimensional problem (See Chapter 5) will be provided later for tests. 3.5 Application to Latent Heat Thermal Energy Storage Systems For the phase-change problems solved in the preceding sections by the source-and-sink method, the boundary conditions have been applied directly to the phase change material. In practice, the phasechange material may be heated or cooled by a fluid which is in motion itself, a conjugate heat transfer problem thus results. In this problem, the fluid temperature is only given at the inlet. Both the wall and the fluid exit temperature are unknown. Two conjugate heat transfer problems will be solved in Sections 3.5 and 3.6. Analysis for fluid flow in a straight parallel channel of phase change material will be studied first. Results for a sample study for this problem will be given in Chapter 5.

PAGE 82

67 3.5.1 Problem formulation A thermal energy storage (TES) unit is sketched in Figure 3.5. This unit can be modeled for analysis as shown in Figure 3.6. It consists of a straight parallel channel surrounded by PCM. Heat transfer fluid (HTF) is forced to flow through the channel and it exchanges heat with the PCM. It is a time-dependent phase-change problem combined with conjugate forced convection. The governing equations for the forced convection in the HTF and the heat diffusion in the PCM must be solved simultaneously. The solution for heat transfer in the PCM has been covered in the previous sections. Only the convection problem and the matching method need be described in this section. Analysis of forced convection for channel flow can be accomplished by solving the Navier-Stokes equations and the energy equation. However, a direct numerical simulation of the Navier-Stokes equations is time consuming. In addition, for a conjugate problem, an i teration is necessary to satisfy the boundary condition at the fluid/PCM interface, which further adds to the complexity of the solution. An alternative approach is taken that uses empirical correlations for the treatment of the heat transfer at the fluid side. This permits the use of convective heat transfer coefficient. While the accuracy of the method depends on the accuracy of the empirical correlations, the uncertainty is comparable with that introduced in the turbulence model. With this effort, the fluid temperature in the channel can be modeled as a function of position and time and analyzed by solving the following equation:

PAGE 83

68 Figure 3.5 A schematic drawing for thermal energy storage unit Figure 3.6 Simplified system for analysis

PAGE 84

69 (3.154) At the inlet, T / (0,<) = T m (<) (3.155) In these equations, Tj is the mixed mean temperature of the fluid, U is the mean velocity, h is the convective heat transfer coefficient, d specific heat of the fluid, respectively. In the formulation above, the Peclet number is taken to be large, so that the heat diffusion in the x direction can be neglected. In the conjugate problem, Eq. (3.154) is to be solved in conjunction with the phase-change problem whose solution has been given previously. If the heat conduction in the axial direction in the channel wall and the heat storage in the wall are negligibly small, heat transfer by convection from the fluid to the wall must be equal to the heat that is supplied to the PCM through its lower boundary. Hence , where y p is the lower boundary of the PCM (see Figure 3.6). Using an electric analogy, the temperature at this boundary can be related to the temperature at the inner surface of the channel wall by the is the half width of the channel, and pj and Cj are the density and the (3.156) relation

PAGE 85

70 T P = T W ~KUK (3.157) where b w and k w are the thickness and the conductivity of the wall, respectively . The heat transfer coefficient h is taken from the literature. In this regard, the Nusselt number for a fully developed laminar flow varies from 7.54 for constant wall temperature to 8.24 for constant wall heat flux condition (Kays and Crawford, 1980). There is a lack of correlation for turbulent flow in a straight channel. Use is thus made of an empirical equation developed by Gniel inski (1976) as (//8)(iZe g -1000)Pr l + 12.7(//8) 1/2 (Pr 2 / 3 -l) (3.158) where / is the friction factor, which is related to the Reynolds number as / = [0.79ln(Re D ) — 1.64] -2 (3.159) In these equations, both the Nusselt number and the Reynolds number are defined as hD. Nu d = U ’ *7 „ UD. Re D ~ u (3.160) where D h is the hydraulic diameter of the channel. It is noted that Eq. (3.158) is strictly valid for a circular pipe with any condition (constant wall temperature or flux) imposed on the boundary. The

PAGE 86

71 Prandtl number is in the range of 0.5 to 2,000, while the Reynolds number is in the range of 2,300 to 5xl0 6 . It is adapted for use in the straight channel by redefining the dimensionless numbers in terms of D h as shown in (3.160). 3.5.2 Solution method In the conjugate problem, the heat flux at the channel wall and the temperature in the PCM and the fluid are all unknown. To satisfy the matching condition (3.156), it is necessary to use an iteration. However, as discussed by Barozzi and Pagliarini (1985), even for a simple steady-state conjugate problem without phase-change, instabilities may be encountered during the iteration process. Convergent results can only be obtained if the initial guesses are sufficiently close to the true values. To alleviate these problems, the heat flux in the third term in (3.154) is re-written in terms of 7j , the temperature at the first grid point above the boundary of the PCM at y p (Figure 3.6). Thus, h(T j T w ) = (3.161) where R" represents the overall thermal resistance per unit area, R" i_ . K ~h + k w (3.162) It is noted that, in the numerical solution, if the interface position

PAGE 87

72 5 is less than Ay, which is the grid spacing in y direction, then the T i i n Eq . (3.161) is replaced by T m , the phase-change temperature, and the Ay is changed to 5. The energy equation (3.154) can then be written as 8T f dT f Tf-T , ~df + U 1h + p fCf R"d = 0 (3.163) This equation can be solved by using a first-order accurate implicit finite difference scheme as rpn + 1 1 ti i + g. 1 0 ^ + ^2 2^1 (3.164) where 9 _U At Ax' a At 2 pfCjR"d The temperature imposed at the surface of the PCM follows as R \ T \ + R 2 : T f R’{ + R‘1 (3.165) where p// 1 i Rl ~h + k~: (3.166) The surface temperature given in (3.165) will then be used as the boundary condition for the solution of the phase-change problem.

PAGE 88

73 3.5.3 Energy conservation te st For the solution of the conjugate problem given in this section, errors mainly come from the matching conditions between the PCM and the HTF. Since experimental data are usually unavailable for comparison, the accuracy of the numerical solution developed in this section will be tested by means of an energy conservation that is based on an open-system thermodynamic analysis. For the heat transfer of channel flow between parallel walls of phase change material without heat and work interaction with the surroundings, the first law of thermodynamics can be written as ^in~E out = AE cv /At (3.167) where the E in refers to the energy rate associated with the fluid flow entering the channel and E out refers to the energy rate associated with the fluid flow leaving the channel. The term on the right hand side represents the time rate of change of the internal energy inside the control volume. For the problem under investigation, there are energy changes in the PCM, the wall, and the fluid. For a time period of t, the left hand side of Eq.(3.167) can be expressed as t E ,U E out = J Qp f c f [T f (0,T)-T f (H,r)]dT T — 0 (3.168)

PAGE 89

74 while the right hand side can be expressed as m f D c n AE cv=T,wl J j pc{T(x,y,t)-T 0 }lxdy y— 0 x=0 H + J [~ LpS(x,tj\h x—0 H x—0 H x + | p f c f ^T f (x,t)-T 0 ]ix (3.169) Here, the energy balance in (3.167) is formulated on the basis of discharge of heat from the PCM. The PCM, the wall, and the fluid are at the same initial temperature T Q , and the channel has a width W. Each PCM plate has a thickness 2D. Under these conditions, the first term on the right of (3.169) refers to the energy associated with the change of the sensible heat of the PCM. The second term refers to the energy associated with the phase change. The third term accounts for the energy change in the containing wall, and the last term refers to that energy associated with the change of the temperature of the HTF inside the channel. For an accurate numerical solution of the problem, the left hand side of (3.167) should be equal to the right hand side of the same equation. This energy conservation will be used to check the accuracy of a numerical experiment to be presented in Chapter 5.

PAGE 90

75 3.6 Application to Encapsulated PCM TES Systems In the preceding section, the phase change material analyzed is in the form of a plane slab. Since the heat transfer fluid is flowing over the surfaces of the slab, the plane of symmetry can be taken to be the half thickness of the slab as shown in Figure 3.6. An alternative form of the phase change material is one that is encapsulated in spheres as shown in Figure 3.7. This section is devoted to the analysis of the phase change in packed spheres. 3.6.1 Problem formulation The system under investigation is sketched in Figure 3.7. A tank is packed with encapsulated PCM spheres. Heat transfer fluid enters through the bottom of the tank, exchanges heat with the PCM by forced convection, and exits from the top of the tank. It is assumed that, in the HTF, the temperature only varies in the z direction while, in the capsulated spheres, the temperature only varies in the r direction as shown in Figures 3.8 and 3.9. Other assumptions that are related to the phase-change analysis earlier still apply. For the analysis of the heat transfer in the HTF, the control volume is taken to have a thickness that is equal to the diameter of the packed sphere (Figure 3.8). The energy equation for the fluid can be written in terms of the mean temperature Tj and mean velocity U as / dT f dT f s pff v K -dr +u -df)+ hA s(T f -T s ) = v, T f = T f (z, t) (3.170)

PAGE 91

76 A schematic drawing for the encapsulated latent heat TES unit Figure 3.7

PAGE 92

77 T f . U 2Rc Figure 3.8 Control volume for fluid analysis Figure 3.9 Boundary conditions for PCM/HTF analysis

PAGE 93

78 where, in addition to those that have been defined previously, T g is the surface temperature of the sphere, V f is the volume of the free space in a packed layer, and is the total surface area of the spheres that are in contact with the fluid. For a single sphere of radius R s , the surface-to-volume ratio is given by the simple relation which also holds for multiple using a void fraction p, defined the total volume V c , Eq. (3.170) 9T f dT f w +u sr + spheres in a packed as the ratio of the can be recast as 3(1 — fi)h VPfCfR-s (T f — T s ) = 0 (3.171) layer. Thus, by free volume Vj to (3.172) This equation is equivalent to that derived by Schumann (1929) and Riaz (1978) for packed rock beds. For the present problem, the initial and boundary conditions for the fluid are given as T f (z,0) = T o , T / (0,t) = T,„(t) (3.173) As for the spheres, the governing equation, initial, boundary, and interface conditions have been given in Section 3.3.3. These equations must be solved in conjunction with (3.172) for heat transfer analysis. As before, the solutions of the fluid and the PCM are to be matched at the surface of the spheres by the relation

PAGE 94

79 In this effort, packed spheres is Chou et al. (1992) h(T s -T f )= -k dT = 9 , p dr\r-R p the convective heat transfer coefficient h obtained by curve-fitting the experimental as (3.174) for the data of Nu =11 + 0.0 6Pe (3.175) where, Pe is the Peclet number ( Re-Pr ), The Reynolds number is defined on the basis of the diameter of the sphere and the average velocity of the fluid as Re = U-(2R S ) V (3.176) The heat transfer coefficient can then be evaluated by means of the Nusselt number as Nu-k f h = ~2R^~ (3 . 177) where kj is the conductivity of the fluid. 3.6.2 Solution method Solution of the conjugate problem for fluid flow in a packed sphere bed follows the same line described earlier for channel flow. In this case, it is assumed that the heat storage in the sphere wall is negligibly small, so that the heat transfer from the fluid to the sphere can be expressed in terms of the temperature difference between the fluid and the first nodal point that is located inside the surface

PAGE 95

80 of the PCM as K T fT s) = — (3.178) For the spherical PCM, the total thermal resistance based on the surface area of the sphere is R" I +. _ _L) + _ J_\ h KSR P R.J TVi RJ (3.179) where k w is the conductivity of the sphere wall, R p is the radius of the PCM, and R s is the radius of the sphere (see Figure 3.9). Equation (3.172) can then be written as dT , dT f 3(1 -u) + + ( r /T .) = » (3-180) As before, in the numerical solution, if the interface position S is greater than , the first nodal point position that is inside the surface of the PCM, but is less than R p , the radius of the PCM, then T l i n Eqs . (3. 178)-(3. 180) is replaced by T m and r l is changed to 5. Equation (3.180) can be solved by using an implicit finite difference scheme as nnn+1 1 f ~ T f + + e 2 r l 1 + Q-y + $2 (3.181) where T fi is the temperature of the fluid at the inlet of the control

PAGE 96

81 volume as shown in Figure 3.8, and n _AtU 0 3(1 fi)At 1_ 2 Jl,» °*-]Tp f c f R s R" (3.182) The temperature at the surface of the PCM can then be evaluated as R'{T l + R^T f R'{ + R'l (3.183) where *i'4+ ^s( 1 V
PAGE 97

82 E,n-Eout = AE cv/ At (3.185) where the E in refers to the energy rate associated with the fluid flow entering the tank and E out refers to the energy rate associated with the fluid flow leaving the tank. The term on the right hand side represents the time rate of change of the internal energy inside the control volume. For the problem under investigation, there are energy changes in the PCM, the shell of the PCM sphere, and the fluid. For a time period of t, the left hand side of Eq.(3.185) can be expressed as t E in ~ E out ~ J" QPfC f [T f (0,T)-T f (H,T)]dT (3.186) T — 0 while the right hand side can be expressed as R n f r AE cv= Epj pi T k(r,t)-T 0 1 r 2 dr ^ r=0 + + ^!> w cJK 3 ,-RliT wk (t)-T 0 ] (3.187) H + rttR 2 Pj c^T f {z,t)-T o y z j=0

PAGE 98

83 Here, the energy balance is formulated on the basis of discharge of heat from the PCM. The PCM, the shell, and the fluid are at the same initial temperature T 0 , and the tank has a height H and a radius R t . The number of spheres in the tank is N which is related to the volume fraction of the PCM sphere as N = (volume fraction of sphere) . vo 1 ume of the tank volume of a sphere Under these conditions, the first term on the right of (3.187) refers to the energy associated with the change of the sensible heat of the PCM. The second term refers to the energy associated with the phase change. The third term accounts for the energy change in the spherical shell of the PCM spheres, and the last term refers to that energy associated with the change of the temperature of the HTF inside the tank. It is expected that, for an accurate numerical solution of the problem, the left hand side of (3.185) should be equal to the right hand side of the same equation. This energy conservation will be used to check the accuracy of a numerical experiment to be presented in Chapter 5.

PAGE 99

CHAPTER 4 TWO DIMENSIONAL PHASE-CHANGE EXPERIMENT In this chapter, an experiment is described for the purpose of validating the numerical solutions described in the previous chapter. The experimental setup, method of measurements, and data reduction are discussed in detail. 4.1 Objective and Basic Considerations The objective of the experiment is to validate the numerical solutions developed in the previous chapter. Since there is lack of data for two dimensional phase change, experimental measurements become necessary. In this effort, the interface motion is visualized and measured. The measured interface position is then compared with the predicted position for validation of the numerical solution. It has been well established in phase-change studies that the interface position provides the best indication for the solution error. Comparison of the temperature field is often inconclusive even if the interface position is conspicuously in error. This is due to the fact that, for a phase-change problem, the boundary and interface conditions can be satisfied exactly. This leaves little error for the temperature. As a result, the temperature provides little resolution for discerning the numerical solution error. 84

PAGE 100

85 To facilitate visualization of the interface position, a test cell of the configuration of a rectangular parallelepiped is constructed. Paraffin wax is chosen as the phase change material because of its very narrow melting range or nearly distinct melting point (Lane, 1986). It also possesses the desired characteristic of becoming translucent to visible light upon phase change from solid to liquid state. The interface position can thus be measured by illuminating the wax with an incandescent light source and recording the interface position on photographic films by means of a camera. The thermophysical properties of the paraffin wax are summarized in Table 4.1 (source: Lane, 1986), in which, the melting point listed is the result of experimental measurements that are specifically performed for the study described in this work. Table 4.1 Thermophysical properties of paraffin wax Properties Solid Liquid Melting point (A) 327 — Heat of fusion (kJ /kg) 266 Specific heat ( kj/kg K ) 2.51 2.95 Conductivity (W/m K) 0.24 0.24 Density (kg/m 3 ) 760 818 In the experiment, the wax is heated electrically. The apparatus is designed so that the boundary condition can be controlled precisely. Provisions are thus made to simulate the prescribed temperature at the heated side of the wax while keeping the other boundaries of the wax insulated. Since no convection has been

PAGE 101

86 accounted for in the previous analysis, convection in the melted wax must be suppressed in the experiment. This motivates the heating of the wax from above, effectively minimizing the free convection. The large density difference between the solid and the liquid wax (see Table 4.1) presents another concern that is related to the expansion of the wax upon phase change. Another provision is thus made to permit the wax to overflow due to volume change. The experiment is now discussed in detail in the section that follows. 4.2 Experimental Setup and Measurement 4.2.1 The setup The test setup is provided in a sketch given in Figure 4.1. It consists of three major components: the test cell (e), the light source (/), and a camera (a), as elaborated below. The test cell is shown in an assembly drawing given in Figure 4.2 and an exploded view shown in Figure 4.3. The test cell measures 12(V)x6(H)x2(D) cm. It consists of a wooden frame (e, Figure 4.3) covered by double-walled transparent lexan windows (c,d) on its front and back sides. To minimize the convection heat loss from these windows, outer windows (c) are separated from the the inner windows ( d ) by a small gap (0.003m). A small Rayleigh number of 91 is found for the gap, providing justification of use of a two dimensional heat transfer analysis for the cell (more justification to follow). To heat the phase change material, a heater (/) is provided at the top of the

PAGE 102

87 e a. Camera b. Temperature control console c. Timer d. Digital thermometer e. Test cell f. Light source g. Power supply Figure 4.1 Test setup

PAGE 103

88 Figure 4.2 Assembly drawing of the test cell

PAGE 104

89 F!<^* e. Frame a. Nuts f • Electric heater assembly b. Gaskets 9Thermocouple wires c. Outer windows h. Bolts d. Inner windows i . Levelling screws Figure 4.3 Exploded view of the test cell

PAGE 105

90 cell. Driven by a DC power supply (g in Figure 4.1), it is specially designed to incorporate several features as detailed below. As shown in the drawing given in Figure 4.4, the heater is composed of three heating elements (6) which are individually wired to their connectors (e). The heater is affixed to the top of a copper plate (c), which keeps the temperature distribution uniform in the transverse direction while maintaining a gradual temperature variation in the other direction. Copper-constantan thermocouples ( d ) are attached to the lower side of the copper plate so that they measure the temperature imposed on the phase change material when the heater is in place in the cell. The four holes (/) located at the two ends of the plate are lined up with identical holes provided at the flanges of the lid (a). These holes are for vent purposes so that, when the phase change material expands, the material can be vented through these holes to the small recess formed between the flange and the wooden frame wall. Paraffin wax has been used as the phase change material. It is poured into the cell through the hole at the base of the frame (e in Figure 4.3). The cell is placed in an upright position with the heater at the top during tests. This causes the wax to be heated from above, a position which effectively suppresses free convection circulation in the cell. As a further precaution for minimizing convection, screws (i in Figure 4.3) are used to adjust the interface position so that the interface is leveled at all times in the course of the experiment. A semi-transparent grid paper is attached to the inner window facing the light source (/ in Figure 4.1). The grid paper serves dual

PAGE 106

91 a e a • Lid b* Heating elements c . Copper plate d • Thermocouples e • Connectors to power supply f . Vents Figure 4.4 Exploded view of the heater

PAGE 107

92 purposes of diffusing the light while also providing grid lines to measure the phase-change interface position. Finally, a camera (a in Figure 4.1) is used to record the images of the wax during phase change . 4.2.2 Test procedure Prior to the experiment, the melting point of the wax is measured by using a thermometer. During the test, the digital thermometer readings are recorded once every 5 minutes and the interface positions are photographed once every 15 minutes. Because of the time delay that is unavoidable in recording temperatures from one thermocouple station to another, provision is also made to correct for the temperature variation for synchronization. During the time when the camera is not being used to record the interface position, the test cell is covered with an insulation sleeve jacket to eliminate the convection heat loss from the cell. 4.3 Sample Results and Experimental Uncertainty Sample pictures taken for the interface positions are shown in Figure 4.5. The interface positions are clearly visible. Also, because of the temperature imposed on the wax being nonuniform, there are slight curvatures of the interface close to two ends. In this effort, the grid paper used has a scale of five divisions per centimeter. Repeated tests of reading the interface position over the background grids provides an uncertainty of position reading to be about 0.5mm.

PAGE 108

93 Figure 4.5 Sample picture for interface position measurements

PAGE 109

94 Since the thermocouples are highly accurate (better than ±0.5°C), errors of the experiment come primarily from the (interface) position reading .

PAGE 110

CHAPTER 5 RESULTS AND DISCUSSION In this chapter, numerical experiments will be provided for the validation of the solution methods described in Chapter 3. Numerical solutions of one dimensional Stefan problems will be discussed first. Results for such problems will be compared with exact solutions and the results published in the literature. Solutions of two dimensional problems follow next and they will be validated by comparing them with some limiting cases and the experimental measurements described in Chapter 4. Finally, numerical results of the latent heat TES systems will be discussed. Numerical Results of One Dimensional Problems Eight test cases (Table 5.1) are presented for verifying the methods developed for the solution of the Stefan problems. Results are compared with analytical solutions as well as the solution methods documented in the literature. Both aluminum and paraffin wax have been tested whose property values are listed in Table 5.2. It is noted that the properties of the wax are taken to be the average of the values of the solid and the liquid state as reported in Hsieh and Choi (1992a). 5.1.1 Test of accuracy and stabilit y Prior to the test of the stability and accuracy for the numerical solutions developed in this work, the superposition concept 95

PAGE 111

96 Table 5.1 Cases tested for one-dimensional Stefan problems Test Case Test Performed Boundary and Initial Conditions (T in K, q in W/m 2 , t in s) PCM and Processes 1 Validation of the superposition method and test of stability and accuracy T U) =383 , T,=T m =327 wax, one-phase melting in a semi-infinite domain 2 Further validation of the superposition method T |u =732+8< 2 and 9 u =6.39x10 6 +3.87x 10 5 < 2 , T,=300 < T m aluminum, two-phase melting in a semi-infinite domain 3 Further validation of the superposition method T UI =336.7 at x=0, and q«, = 0 at x=H, T, = 317.3 < T m wax, two-phase melting in a finite domain 4 Stefan-number effect on interface position in finite domain T UJ = 336.7 at x=0, and q w = 0 at x=Il, T,=327, 278.5, and 230 wax, two-phase melting in a finite domain 5 Test of oscillating interface for long-time solution T u ,=354+20«n(ir
PAGE 112

97 Table 5.2 Properties of PCMs analyzed Aluminum Paraffin Wax T m [K] 932 327 L [kj /kg] 376.56 266 k [W/mK] 200 0.24 p [kg/m 3 ] 2,710 789 c [kJ/kgK] 1.2 2.73 is first tested by applying it to the solution of a Stef an-Neumann problem described in Test Case 1 (5
PAGE 113

98 Position, x (cm) Figure 5.1 Temperature profiles of the superposition method

PAGE 114

Interface position, S (cm) Interface position, S (cm) 99 Figure 5.2 Test of accuracy and stability using different At and Ax

PAGE 115

100 The bottom figure gives the interface positions for At equal 60 s. Here, the number of grid points varies from 10 to 90. Again, the results are good; in fact, they are stable in both tests. It is noted that the implicit method of solution employed here permits use of large time steps together with a coarse grid, which is not possible for explicit method of solution. Efforts have also been made for identifying the time when the largest numerical error occurs. As shown in Figure 5.3, with Nx taken to be 90 and At varying from 60 to 600 s, the largest error (107.) is always found at the first time step, which is consistent with the findings for the SSM as reported by Choi and Hsieh (1992). This percentage error is diminished rapidly (less than 17.) toward large time, a distinct feature for this method. Case 2 (Table 5.1) provides a more thorough test of the superposition method based on the development of the numerical solutions given in this work. In this case, time-variant temperature and flux conditions are imposed on a semi-infinite medium (aluminum), which is initially subcooled (see Table 5.1). A two-phase melting thus occurs and the numerical results are compared with the SSM results reported by Hsieh and Choi (1992a) and Choi (1992). For Case 2, 60 nonuniform grids are used for the solution of U, while the time step is taken to be 1 s. The interface-position-versustime results are shown in Figures 5.4 and 5.5. In Figure 5.4, the medium is imposed with a time-variant temperature condition, while in Figure 5.5, the medium is imposed with a time-variant flux condition. Numerical results are in good agreement with the SSM results obtained

PAGE 116

Errors in interface position, 101 Figure 5.3 Errors in the interface positions for different At

PAGE 117

Interface position, S (cm) f Interface position, Figure 5.5 Interface positions for time-variant flux condition

PAGE 118

103 by Hsieh and Choi (1992a), providing further assurance of the success of the superposition method. The third test case (Case 3 in Table 5.1) deals with a problem in a finite domain. In this case, the slab is initially subcooled in a solid state with one side (x = H) insulated and the other side (x = 0) imposed with a constant temperature condition. To compare the results with the published data given in the literature (Dursunkaya and Nair, 1990), the following dimensionless groups are introduced Sb = Tm-Tj T u,~Tm x-2H' — T-T 7 1 m ~T:-T ’ i m i = ^4 H 2 (5.1) where Ste is the Stefan number and Sb is a subcooling number. Paraffin wax is again used for tests, which is taken to be 1 cm thick divided into 10 grid points. The initial and boundary conditions listed in Table 5.1 are specially selected to yield a Stefan number of 0.1 and a subcooling number of 1, which have been used by Dursunkaya and Nair in their original paper. Results are shown in Figures 5.6 and 5.7. In Figure 5.6, the dimensionless interface position is plotted as a function of the dimensionless time. The present results are in good agreement with Dursunkaya and Nair’s data. As shown in this figure, it takes a value of t of 5.3, which corresponds to 4600 s, for the slab to be melted completely. Figure 5.7 gives the dimensionless temperature history at the insulated side (x = H) of the wall. The wall temperature

PAGE 119

Dimensionless wall temperature Dimensionless interface position 104 Figure 5.6 Interface positions in a finite domain Figure 5.7 Temperature history at the insulated wall (x=H)

PAGE 120

105 remains constant for a short while after time zero, signifying the time period when the heat has not penetrated through the depth of the wall. Then, the dimensionless temperature drops, indicating a temperature rise at the insulated side. Again, the results are in good agreement with the published data in the literature. 5.1.2 Subcooling effect on the interface position Case 3 provides a good test for the assessment of the effect of subcooling on interface motion. The initial temperature of the wax in Case 3 is varied (230, 278.5, and 327 K) as listed in Case 4 (Table 5.1) and the results are plotted in Figure 5.8. Here, Neumann solutions in a semi-infinite domain are also attempted and their results plotted in curves for comparison. As shown in this figure, for a given subcooled condition, the interface front in a finite domain (symbols) first follows that in a semi-infinite domain. Later, they depart from each other. Eventually, all curves for the finite domain (symbols) move at same speed, which is evidenced by their identical slope at large time. These phenomena can be attributed to the fact that, since the medium is subcooled initially, it takes time for the heat supplied at x = 0 to penetrate through the depth of the wall. Consequently, during the initial period of time to heat the wall, the heat transfer in a finite domain is equivalent to that in a semi-infinite domain with the same degree of subcooling. Once the heat penetrates through the wall and reaches the insulated boundary, the temperature at this boundary starts to rise. Then, the interface in a finite domain moves in a heated medium at a higher speed. From another

PAGE 121

Dimensionless interface position Figure 5.8 Interface positions for different subcooled condit

PAGE 122

107 perspective, the medium in a finite domain with back side (x = H) insulated may be considered to be pre-heated prior to melting. Thus, for the final stage of melting, the medium is fully pre-heated to its melting point; the problem reduces to a one-phase melting in a StefanNeumann problem. 5.1.3 Solution with oscillating interfaces An interesting case arises when the wall is initially at its melting point but in solid state. For time greater than zero, the boundary at x = 0 is imposed with a sinusoidal temperature condition which is at a level higher than the melting point of the medium. The boundary at x = H is imposed with a constant temperature condition which is lower than the melting point of the medium. Case 5 represents such conditions and the results are plotted in Figure 5.9, where two curves are plotted for comparison. The dashed line represents the case when the boundary at x = 0 is raised to a constant temperature of 354 K. The solid curve represents the boundary imposed with a sinusoidal temperature, which differs from the previous condition by that sine term (see legend). For both curves, the conditions at the other side (x-H) are identical. Clearly, because of the boundary conditions imposed, the interface position rapidly reaches a quasi-steady state. The solid curve fluctuates about the dashed curve, which eventually reaches the mid-plane of the slab. Using this dashed curve as a basis, there is approximately a 10 degree angle of phase delay. The interface oscillates at the same frequency as the imposed temperature cycle. The long-time cyclic solution presented here is believed to be the first

PAGE 123

Interface position, S (cm) 108 Figure 5.9 Test of oscillating interface

PAGE 124

109 that addresses quasi -steady temperature variation accompanied by a phase change. A more detailed follow-up study is being made by Nyros and Hsieh (1993). 5.1.4 Solution with multiple interfaces Case 5 deals with a situation where the phase change is initiated from one side of the wall. A variation of that is for interfaces initiated from two sides of the wall. Two cases (Cases 6 and 7) are included to address such problems as listed in Table 5.1. Case 6 is for solidification of a 5 cm aluminum slab initially in liquid state at its freezing point. Both sides of this slab are imposed with constant temperature conditions which are at the levels lower than the freezing point of aluminum. Phase change is thus initiated from both of the boundaries and the interfaces move inward with time as shown in Figure 5.10. Based on the Stefan number defined in (5.1), the boundary at x = 0 is imposed with a Stefan number of 0.637, while the boundary at x = H is imposed with a Stefan number of 0.319. For this problem, the interface motion can be analyzed by means of a one-phase Stefan-Neumann solution. The superposition method is tested for this case as shown in Figure 5.11. The interface positions are plotted versus time as shown in Figure 5.12, where the numerical results are compared with the Stefan-Neumann solutions. The two interface fronts move independently of each other until they merge at 2.88 cm. The results are in good agreement with the analytical solutions.

PAGE 125

110 Figure 5.10 Sketch of double interfaces generated by conditions imposed on two boundaries

PAGE 126

Ill Figure 5.11 Temperature profiles for a doubl e-interface problem

PAGE 127

Interface positions, S (cm) 112 Figure 5.12 Interface positions of a doubleinterface problem

PAGE 128

Interface position, S (cm) 113 T(0, t) 0 S 2 Si H X Figure 5.13 Sketch of double interfaces generated by conditions imposed on one boundary Time, t (s) Figure 5.14 Interface positions of a melting-freezing problem

PAGE 129

114 Case 7 deals with a different situation in that the phase changes are initiated from the same side of the boundaries in a plane wall. A thick aluminum wall is considered, which is initially at its melting point but in solid state. At time greater than zero, the boundary at z = 0 is imposed with a time-variant temperature condition while the other at x = H is insulated (see Table 5.1). Because of the conditions imposed, the medium first melts and then re-freezes as shown in Figure 5.13. The interface histories are plotted in Figure 5.14. Here, because of the rapid temperature drop at the boundary for time greater than 100 s, the freezing front rapidly overtakes the melting front. They eventually merge at 153 s, when the medium return to its solid state. 5.1.5 Melting of a sphere In all the tests performed above, the phase change has been considered in a Cartesian system. Case 8 deals with a one dimensional phase change in a sphere which is initially at its melting point but in solid state. The surface is suddenly heated to a constant temperature greater than the melting point as described for this case in Table 5.1. There is no exact solution for this problem. Hill (1987) solved this problem by an integral method, and derived the lower and upper bounds for the interface position as given by the following formulae : **36 ( 5*02 ( 1 -^) 2 [( 1 + 2 / ? )( 1 +^) + 2 (/?l ) 5 2 l 2 60(S
PAGE 130

Dimensionless interface position 115 Figure 5.15 Interface positions of sphere melting

PAGE 131

116 and 1 6(i
PAGE 132

117 Table 5.3 Cases tested for two-dimensional Stefan problems Test Cas< Tests Performed Boundary and Initial Conditions (T in K, q in W/m 2 , t in s) PCM and Processes (x and y in cm ) 1 Test of stability and accuracy ^=397-100 |§| at t/=0, other boundaries insulated, T, =T m = 327 wax, one-phase melting in a finite domain of -5 1800 where /(x, 1) = 8.9X10' 5 (1.5 --§) tl (18001 J 2 ) T,=327 K < T m wax, two-phase melting and freezing in a finite domain of 0
PAGE 133

118 5.2.1 Test of accuracy and stability Case 1 (Table 5.3) tests the stability and accuracy by solving the phase change in a two dimensional domain sketched in Figure 5.16. The bottom boundary is imposed with a space-variant temperature condition, while the other three sides are all insulated. Wax initially at melting point in a solid state is used for experiments. To test the stability of solution, a 40x30 grid is used to solve the problem in x and y direction. The time step varies from 30 to 600 s, which corresponds to dimensionless time steps (Fourier modulus) of 0.5 to 10. Figure 5.17 gives the numerical results for interface position histories at x = 0 and x = H/ 2. Two groups of curves are plotted. In each group, the solid line is an analytical solution of a Stefan-Neumann problem imposed with a uniform temperature condition. This provides for the upper bound for the interface position at x = 0 and the lower bound for the interface position at x = H/ 2. The dashed curves give the numerical results for the two dimensional problem. These dashed curves stay within the bounds. Most importantly, all dashed curves are smooth and converge nicely irrespective of the values of the time step tested. Figure 5.18 gives the interface position at a time equal to 3000 s. Four different time steps are tested and the curves are smooth and convergent. No instability is encountered even with the use of a time step as large as 600 s. It is thus concluded that the explicit treatment of the partial differentiation in the x-direction of the V problem, Eq.(3.107), does not introduce any instability in the solution. In two-dimensional problems, the time steps are determined

PAGE 134

119 T w =T w( X ' t ) Figure 5.16 System used in two-dimensional stability test

PAGE 135

Interface position, S (cm) £ Interface position. S (cm) 120 re 5.17 Interface positions at different times at x=0 and x=H /2 Figure 5.18 Interfac e position curves at <=3000s for different At

PAGE 136

121 strictly on the basis of the accuracy rather than the stability. This also serves to verify numerically the comments made earlier at the end of Chapter 3. It should be noted that the stability of the numerical solution demonstrated in this work is independent of the number of the grid points used in the x direction. A reduction of the number of the grids by a factor of four (Nx=10) and application of a time step of 300 s yield data (circles) which fall right on the curves in Figure 5.18. This also serves to demonstrate the efficiency of the method. 5.2.2 Long-time solutions Case 2 is intended for testing the convergence of the results to a long-time solution. As summarized in Table 5.3, a column of medium of infinite height initially subcooled is used for tests. A uniform temperature is imposed at the bottom boundary while the lateral boundaries dissipate heat to ambient by convection (Figure 5.19). Heat transfer in the column is thus two dimensional. There is no exact solution of this Stefan problem. However, a steady-state solution can be derived as follows T T oo _ V' 2 (Pm +Bi2 ) sin P m T ” ~ T °° £l( P 2 m +Bi 2 ) +Bi P m 2V0JH COS (2 xfim/H) (5.4) where Bi — hH /2k is the Biot number, and /? m ’s are roots of Pm tan P m = Bi ’ m= 12 .. (5.5)

PAGE 137

122 Figure 5.19 System used for long-time solution

PAGE 138

Interface position, S (cm) 123 Figure 5.20 Comparison of numerical and analytical solutions for a two-dimensional problem at long time

PAGE 139

124 This steady-state solution can be used for checking the interface position at long time by setting y to S(x) and T to T m in (5.4) as T m T oo _ 2 (/ ? m + B * 2 ) sin P m T tv~ T oc £?' 1 (P 2 m +Bi 2 )+Bi P m cos(2x(3jH)e 2S0 ™/ H (5.6) The interface position results are shown in Figure 5.20. Here, because of the temperature symmetry, only one half of the width of the medium is used for presentation. As shown in this figure, the interface position curves determined by the source-and-sink method (Section 3.4) move upward with time. The topmost curve, representing the steady state limit, is in close agreement with the analytical results (circles) predicted by using (5.6). In this analytical efforts, sufficient number of terms have been taken in the partial sum in the infinite series to ensure convergence (10 -6 ). The two dimensional test is thus successful at long time as shown clearly in this example. 5.2.3 Comparison with experimental results What remains to be done is to test the two dimensional solution at short time. This is accomplished by comparing the short-time experimental results with those obtained using the numerical methods given in this work. A composite diagram is constructed as shown in Figure 5.21. Here, the data at the top refer to the imposed temperature conditions at the surface of the paraffin wax as measured

PAGE 140

Interface position, S (cm) Boundary temperature, TV (C) 125 Figure 5.21 Comparison between numerical results and experimental data for two-dimensional phase change

PAGE 141

126 from the experiment described in Chapter 4, while the bottom figure shows the positions of the phase-change interface. The y axis for this lower figure is intentionally drawn up-side-down so that it confirms with the analysis (Figure 3.2). It is noted that, in the computation, the ^-direction is discretized into 24 elements and, once again, sufficient number of terms are taken in the Green’s function to ensure convergence (1(T 6 ). As shown in this lower figure, the agreement between the experiment (symbols) and the analysis (curves) is good. There are small discrepancies only in the middle section of the curve which can be attributed to the curvature of the phase-change interface position. In the experiment, the leveling screws described in Section 4.2 are used to adjust the slope of the middle section of the interface curves. A leveling of this section of the curves upsets the slope at two ends. As a result, the convection can not be totally eliminated from the test cell, causing the slight disagreement between the prediction and the experiment. 5.2.4 Test for coalescing fr onts A general case of multiple coalescing fronts is tested as Case 4 (Table 5.3). Wax at melting temperature in solid state is tested which is imposed with a positionand time-variant condition at t/ = 0 (Figure 5.16) as follows: t < 1800 s 327 + 8. 9xl0 _5 (1.5 — j|)(1800<-* 2 ), rjx,o = | 327 8 . 9xl0 _5 ( 1 .5 -jj) (1800* t 2 ) , t > 1800 s (5.13)

PAGE 142

127 a o CO a o CO a o CO a o CO x (cm) Figure 5.22 Interface position curves in a two-dimensional coalescingfront problem

PAGE 143

128 Here, the sign for the second term is changed so that, for time less than 900 s, the boundary temperature decreases in the * direction yet the local boundary temperature increases with time. However, for time greater than 900 s, this local temperature decreases with time. Since the first term in both conditions refers to the melting point of the wax ’ T w ( x >0 ls greater than the melting point for time less than 1800 s while it is less than the melting point for time greater than 1800 s. The wax thus melts and re-freezes unevenly in the x direction, giving rise to a two dimensional multiple-coalescing-front problem. The wax melts and re-sol idifies as shown in Figure 5.22. No literature data have been available for comparison. Yet, the two dimensional interface motions are clearly visible due to the positionvariant temperature conditions imposed. 5.3 Application to Two Dimensional Latent Heat Thermal Energy Storage Systems In this section, numerical results of a sample latent heat TES unit will be presented. The unit is composed of stacked rectangular PCM plates with fluid passages in between (see sketch in Figure 3.5). The input data for analysis has been compiled in Table 5.4. The unit is precharged to 80°C . Since the melting point of the PCM is 72°C, the PCM is initially in a liquid state. Cold water at 23°C is passed through the unit at three different flow rates: 1 gpm, 2gpm, and 3 gpm. The total amount of energy pre-stored in the unit is 63.5 MJ. In the numerical tests, calculation ceases as soon as the exit water temperature drops down to 40°C.

PAGE 144

129 Table 5.4 Input data for TES analysis PCM Na 3 P0 4 • 12 H 2 0 Conductivity, k ( W/m-K ) 0.6 Heat capacity, c ( kJ/kg-K ) 1.7 Density, p (kg/m 3 ) 1,630 Latent heat, L ( kJ/m 3 ) 110 Melting point, T m (°C) 72 HTF — Water Conductivity, k f (W/m-K) 0.58 Heat capacity, c f (kJ/kg-K) 4.19 Density, p f (kg/m 3 ) 1,000 Overall initial temperature, T Q (°C) 80 HTF inlet temperature, T f (0,t) (°C) 23 Volume flow rate, Q (gpm) 1,2,3 Other specifications — PCM length, H (m) 2 PCM width, W (m) 0.5 PCM half thickness, D (m) 0.02 Channel width, 2 d (m) 0.01 Number of PCM plates, M 4 Containing wall thickness, b w (m) 0.0

PAGE 145

130 The accuracy of the numerical simulation is tested by means of the energy conservation described in Section 3.5.3. A comparison of the left and right hand side of Eq.(3.167) is shown in Figure 5.23. For all the flow rates from 1 gpm to 3 gpm, the maximum error is found to be about 37. at the end of the process, which provides a good indication of the accuracy of the present model. The water exit temperatures during the discharge process are plotted in Figure 5.24. For a flow rate of 1 gpm, the unit runs for 170 minutes before the exit temperature drops down to 40°C. For a flow rate of 2 gpm, it works for 95 minutes; and for 3 gpm, it works for 70 minutes. Out of the 63.5 MJ energy initially stored in the unit, the energy extracted from it is 55.5, 54.3, and 52.9 MJ for the flow rate of 1, 2 and 3 gpm, respectively. It is noted that, in the temperature curves shown in Figure 5.24, the slope of the curves changes erratically. This is believed to be the result of the release of the latent heat from the PCM. The temperature distribution in the PCM at 2-gpm flow rate is shown in Figure 5.25. It should be noted that the plots are re-scaled to make them easy to read. The topmost figure shows the position of the isotherms in the PCM after 20 minutes of energy discharging. At that time, the lower region of the PCM is solidified where there is a large temperature gradient (see the dense temperature contours). The upper region remains in the liquid state. The solid/liquid interface can be identified by the sharp change of the density of the isotherms. At 40 minutes (second figure), the left part of the PCM is completely solidified and is uniform at phase-change temperature. The liquid

PAGE 146

131 Figure 5.23 Test of energy conservation

PAGE 147

132 Time, t (minute) Figure 5.24 Water exit temperature at different flow rates

PAGE 148

133 Figure 5.25 Isotherms in a half PCM plate

PAGE 149

134 region diminishes in size until the PCM is completely solidified (bottom figure). After that, only sensible heat remains in the unit. 5.4 Application to Encapsulated PCM TES Systems Numerical results for the encapsulated PCM TES system will be discussed in this section. The input data for the PCM, the HTF, and the working conditions are the same as those listed in Table 5.4. Only the specifications are different, which are listed in Table 5.5. The parametric study will include the effects of the HTF flow rate, and the diameter and the thickness of the PCM ball on the thermal performance of the TES. Table 5.5 Specifications of encapsulated TES unit Height of tank, H (m) 1.1 Radius of tank, R t (m) 0.21 Overall void fraction, p 0.4 Diameter of sphere, 2 R s (cm) 2, 3, 4 Shell thickness, R s -R p (cm) 0.0, 0.2, 0.5 The accuracy of the numerical results is evaluated first by means of the energy conservation as described in Section 3.6.3. In this test, the diameter of the balls is taken as 3 cm, and the thickness of the shell is 0.2 cm. A comparison of the two calculations

PAGE 150

135 of the energy change from Eqs.(3.186) and (3.187) is shown in Figure 5.26. It is found that, for all the flow rates tested (1 gpm to 3 gpm), the maximum error is approximately 3%. The numerical results are thus accurate. Parametric study results will be discussed next. In Figure 5.27, the water exit temperature is plotted versus time for three flow rates. As expected, for a lower flow rate, the unit operates longer; for a higher flow rate, the unit operates shorter. However, the total amount of energy extracted from the unit remains almost the same. The effect of the ball diameter on the water exit temperature is shown in Figure 5.28. For the three diameters tested (2 cm to 4 cm), the thickness of the shell is taken to be constant at 0.2 cm. As expected, for a ball of small diameter, because of the presence of the shell, the volume fraction of the PCM is smaller. As a result, less energy can be stored in and extracted from the ball. As the ball diameter is increased, the volume fraction of the PCM increases correspondingly. However, a large ball is penalized for the smaller surface-to-volume ratio. The heat transfer rate is thus lower, as shown in Figure 5.28. The effect of the shell thickness is shown in Figure 5.29. It may be said that the main drawback of using a thick shell is the reduction of the volume fraction of the PCM. For example, for a ball of 3 cm diameter, when the shell thickness is increased from 0.2 cm to 0.5 cm, the volume occupied by the PCM is reduced by 54.57„. This greatly reduces the energy that can be stored in the unit, as clearly shown in Figure 5.29.

PAGE 151

136 Figure 5.26 Test of energy conservation

PAGE 152

Water exit temperature, T (C) 137 Time, t (minute) Figure 5.27 Water exit temperature at different flow rates

PAGE 153

138 Figure 5.28 Water exit temperature for different ball diameters

PAGE 154

139 Figure 5.29 Water exit temperature for different shell thickness

PAGE 155

140 Radius, r (cm) Figure 5.30 Temperature profile in the PCM and the shell

PAGE 156

141 It is also informative to compare the temperature distribution in the PCM ball (see Figure 5.30). Here, the topmost curve refers to a situation when the PCM is still in the liquid state since there is no abrupt change of the temperature slope in the ball. Phase change takes place at 3 though 9 minutes. The interface moves inward with time. At 11 minutes, the PCM has completely solidified. Only sensible heat is being extracted from the ball. 5.5 Visu alization of Numerical Results To graphically visualize the operation of the latent heat TES system, two animating post-processors are developed using a QuickBASIC language. Sample outputs of the post-processors from a computer screen are shown in Figures 5.31 and 5.32. From the graphics, the interface position, the flow rate, the exit water temperature, the available energy left in the unit, and other system information can be viewed simultaneously.

PAGE 157

142 Figure 5.31 Post-processor for a two dimensional TES unit

PAGE 158

143 Figure 5.32 Post-processor for an encapsulated TES unit

PAGE 159

CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS Based on the analysis and results of the present work, the following conclusions can be drawn: 1. A source-and-sink method has been developed for the solution of multi -dimensional phase-change problems. The method has been shown to be accurate by comparison with exact analysis and experimental measurements. In the source-and-sink method, only one set of equations is required to solve for the temperatures in both liquid and solid phase regions. 2. The success of the present method relies on the use of superposition technique. The validity of the superposition method for the solution of Stefan problems with constant and equal properties in different phases has been rigorously analyzed. In the superposition method, the solution of the Stefan problem was superimposed by a finite difference method, which was used to calculate the temperature distribution due to the initial and boundary conditions, and an analytical solution, which was used to track the interface position and to update the temperature due to the motion of the interface. For multi-dimensional problems, a hybrid scheme was developed which combines an analytical solution in one direction and a finite difference method in another direction. With these efforts, a time 144

PAGE 160

145 marching scheme was developed with which both short-time and long-time solutions can be obtained. The stability of the solution methods for one and two dimensional Stefan problems was analyzed and numerically tested. No instabilities were encountered in the tests even for very large time steps. This provides assurance of stable and converged solutions. 3. In the hybrid method of solution, the subcooling and superheating conditions can be handled conveniently. It is not necessary to divide the solution into pre-melting, melting, prefreezing, and freezing stages. Furthermore, the initial condition can be uniform and nonuniform, and the domain of solution can be finite and semi-infinite. In this work, the effects of subcooling on the interface positions in a finite medium have been studied. Problems with oscillating interfaces resulting from time-variant boundary conditions have also been solved. 4. With the use of multiple source and sink method, multiple interfaces resulting from periodic melting and freezing as encountered in the discharging/recharging processes of a latent heat TES unit were successfully calculated. Coalescing fronts in a two dimensional melting and freezing problem with multiple phase regions were also investigated. 5. A two dimensional phase-change experiment was conducted for testing the accuracy of the short-time solutions of the numerical method. An analytical solution was used for the test of accuracy for long-time solutions. In both cases, the interface positions were used

PAGE 161

146 for comparison. Numerical results were in good agreement with the experimental measurements as well as the analytical solutions. 6. Applications of the solution methods to the simulation of latent heat thermal energy storage units were discussed. Two conjugate heat transfer problems were solved for the prediction of thermal performance of two TES units. The first one was composed of parallel channels formed by stacked PCM plates, and was modeled as two dimensional conjugate problems. The second was an encapsulated TES unit consisting of packed PCM spheres. In both cases, the accuracy of the numerical method was tested by means of energy conservation based on an open-system thermodynamic analysis. It has been shown from the numerical results presented in Chapter 5 that the present methods were accurate and efficient. Parametric studies of the encapsulated TES unit were also attempted. The effects of the flow rate, the diameter and the shell thickness of the PCM sphere on the thermal performance were analyzed. Two post-processors for graphical visualization of the TES system performance were also developed. The following are recommended for further work: 1. The source-and-sink method with superposition technique will be extended to the solution of three dimensional phase-change problems such as laser cutting, drilling, and surface hardening. 2. The properties of the PCM in the present analysis have been considered as constant and equal in different phase regions. In the

PAGE 162

147 future work, unequal properties for different phase regions will be solved by using double sources and sinks as analyzed in Kolodner (1956). The variation of the thermal conductivity with temperature will also be accounted for by using Kirchhoff transformation if convection does not appear at the boundary. 3. In the present analysis, the source-and-sink method requires the formulation of the GreenÂ’s functions. In the case when the GreenÂ’s function can not be analytically formulated, a unification scheme that combines the source-and-sink method and the boundary element method reported recently (Hsieh, 1993) can be used. In the boundary element methods, the GreenÂ’s free space solutions can be used to replace the GreenÂ’s functions which strongly rely on the geometry of the domain as well as the boundary conditions.

PAGE 163

APPENDIX A DERIVATION OF EQUATIONS (3.52) AND (3.53) In Cartesian coordinate systems, the decomposed temperature V is given by (see Eq.(3.46)) ‘l v ( x i t l)=rj } T ^ G{x,t J I 5(r),r)dr (A.l) T=t o In a finite domain, the Green’s function for a temperature condition imposed at x = 0 is given by (see Eq.(3.50)) OQ 2 / G(x,t\ x',t) =^2,jje 0nt T ^sin(P n x)sin(P n x') n=l (A. 2) Here , 7p for a temperature condition at x = H (2n-l)ir 2 jJ > f° r n flux condition at x = H /?„=' Substituting (A. 2) into (A.l) yields the following for soluti Lon V ( X ' *l) ~ ThY^ B ri( t l) sin (l 3 n X ) where n— 1 (A. 3) T — t (A. 4) 148

PAGE 164

149 For a small At, it is assumed that the interface velocity is constant from t 0 to t x . Thus, S ( T ) =S(t 0 ) + v(r-t 0 ) (A. 5) where v = It follows that dS = vdr, and Eq.(A.4) can be transformed as *»('.) = j {e< S /\i„^ nS )ys S(>„) This equation can be integrated and the result is given by (A. 6) n( l) ~ ' {T<*P 2 nM 2 +/3 2 n [\ a ^ /v r n ^nSlh))-0nCos(,fi^ l» e a P n S ( t o)/ v (“ P 2 Jv) 2 +P 2 n (a P 2 Jv)sin{p n S(t 0 ) ) (3 n cos(f3 n S(t 0 ) ) (A. 7) After rearrangement, Eq.(A.7) can be written as S " (
PAGE 165

150 Here , G{x,t\x\T)-C Q +^2^e a/3 " (t T) cos((3 n x)cos(f3 n x') (A. 9) n=l C 0 = ' H ( 2n-l)7r 2 H , for a temperature condition at x = H, for a flux condition at x = H. H Substituting (A. 9) into (A.l) yields the following for solution y(*. *i) = C o•ft At + jft'52 C n( t l) co *(P n x ) n— 1 (A. 10) where c n (< 1)= ^jr{ e a/3 " (tl r)cos (/ ? n 5 ( r ))} dT (A. 11) T — t By following the same line as derived in (A.5)-(A.6), Eq.(A.ll) can be written as S'(tj) C n (t l) = ea ^ )/v \ {e a ^/\ os ^ nS) } dS S(t 0 ) (A. 12) This equation can be integrated and the result is given by G n (h) ~ e -"tiSM/v f e a ^*i)/ v (<*Pn / v )°os ( Pn S Vl ) ) +P n sin ( /? n 5(
PAGE 166

151 After rearrangement, Eq.(A.13) can be written as Cn(<1 ) _ (a/?£) 2 I(/? nV ) 2 |[°^n co *(^„5(*i)) +vp n sin(p n S(t 1 ))] -e a/3 n(h * 0^/32 \<*Plcos(0 n S(t o )) + ^„ S m(/? n 5(< o) )]| This completes the derivation of Eq.(3.53). (A. 14)

PAGE 167

APPENDIX B DERIVATION OF EQUATION (3.68) In a spherical coordinate system, the decomposed temperature V is given by (see Eq.(3.63)) V ( x ’ t i)=%^ ^ S 2 (t )G(x,t 1 | S(t),t)<1t (B. 1) T — t O For a solid sphere, the Green’s function for a temperature condition imposed at the sphere surface is given by (see Eq.(3.64)) Here , oo 2 1 1 x', t) = -i-j ^2 -fi-e r ^m(/? n x)sm(/? n x') n=l s 3 Pn R e (B.2) Substituting (B.2) into (B.l) yields where _ 2 L v' n ,, >n(3 n x) v l‘' t *> = tt E».W * n=l f _ 2 = J 5(r)e' Q/3 n (t r) s in(/?„5(r))dr T = £ (B.3) (B.4) Following the same line as derived in APPENDIX A, for a small time step At, the interface position can be written as 152

PAGE 168

153 5(r) _ 5(t 0 ) + v(r — < Q ) , dS = vdr Then, Eq.(B.4) can be transformed as 5 ( tl ) «<„> From an integration table, one can get J xe ax sin(bx)dx = f^ 2 [ a • sin(bx) 6 • cos(6x)] 2 [(a 2 — b 2 )sin(bx) — 2ab • cos(6a;)J e
PAGE 169

n o o APPENDIX C FORTRAN PROGRAM FOR TWO DIMENSIONAL PHASE CHANGE PROBLEMS CC C C C C C C C C C C C c= c* Numerical Solution of Phase Change Problems One/Two-dimensional phase change problems, Combined FDM and SSM method, Time marching scheme, Combined Newton-Raphson/Bi section for searching for S, One/Two interfaces, With/without free convection on the vertical boundaries, Finite domain. 1st. kind B.C. imposed at x=0, insulated at x=b Time/space variant boundary condition ============—— PROGRAM MAIN ========== = C c c c c c c c c c c c ==c PARAMETER (MAXX=100 , MAXY=50 , MAXX2=MAXX*2 , MAXT=MAXX*MAXY , MM2=MAXY*10 ) DIMENSION Y(MAXT) ,T(MAXT) ,TO(MAXT) ,TV0(MAXX) , TW(MAXX) ,S(MAXX2) ,S0(MAXX2) , DUM2(MM2) , ITFS(MAXX2) COMMON /ALL/XDIM , YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , II FG COMMON /IO/IREAD, IOUT, IREST, IDUMP CHARACTER*40 FREAD, FOUT, FREST, FDUMP , TITLE IREAD = 2 IOUT = 10 IREST = 98 IDUMP = 99 =READ INFO. 900 WRITE(*,*) 5 Input date file ?’ READ( * , 1000) FREAD 0PEN( IREAD, FILE=FREAD,ERR=910) GOTO 920 910 PRINT*,’ Open error, No file found' ’ GOTO 900 920 CONTINUE 154

PAGE 170

155 READ( I READ , 1000 ) TITLE RE AD (I READ,*) READ (I READ, 1000) FOUT READ (I READ,*) READ( I READ , 1000) FREST READ (I READ,*) READ ( I READ ,1000) FDUMP READ (IREAD,*) READ( IREAD , * ) CK , CP , RHO , HFG , TINIT , TMELT READ( IREAD,*) READ( IREAD,*) HCOV, TAIR READ( IREAD,*) READ( IREAD,*) XDIM, YDIM , NX, NY RE AD (IREAD,*) READ ( IREAD , * ) DTT, TOUT , TSTOP , KREST C C== Initial ization . NXY = MAX(NX,NY) ITOUT = TOUT TIK = TINIT + 273. TMK = TMELT + 273. TINF= TAIR + 273. PRINT* , Â’ ====== INITIAL I ZATION========== Â’ CALL INIT(NX, NY, KBC,IBC, FREST, KREST, TSTART, TIK, Y,T, SO, S,TW) C C==Main solution procedure. CALL PROCES ( NX , NY , NXY , KBC , I BC , ITFS , TSTART , ITOUT , TSTOP , 1 DTT , HCOV , TINF , TMK , TIK , Y , T , TO , S , SO , TV , TWO , 2 DUM2( 1 ) , DUM2(NXY+1 ) ) C C==Output OPEN (I OUT, FILE=FOUT) CALL OUTPUT(NX, NY, TSTOP, Y,T) CLOSE (I OUT) CALL DUMP(NX, NY, FDUMP, TSTOP, SO, S,Y,T, TV) 1000 FORMAT (A40) STOP END C C========== SUBROUTINE INIT. ============= C Initialization. SUBROUTINE INIT(NX, NY, KBC, IBC, FREST, KREST, TS, TIK, Y,T,SO,S,TV) PARAMETER (MAXBC=50) DIMENSION SO (NX, 2) ,S(NX,2) ,TV(NX) ,Y(NX,NY) ,T(NX,NY) COMMON /ALL/XDIM , YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , HFG COMMON /BC/XBC(MAXBC) ,TBC(MAXBC,2) ,TIB1 ,TIB2 COMMON / I 0 /IREAD , I OUT , IREST , I DUMP CHARACTER*40 FREST CSCAL =0.01 XDIM = XDIM*SCAL YDIM = YDIM*SCAL

PAGE 171

O Cl o o 156 IF(NX.EQ.l) THEN DX = 1. ELSE DX = XDIM/FL0AT(NX-1 ) END IF DY = YDIM/FL0AT(NY-1 ) ALPHA = CK/ (RH0*CP) IF(KREST.EQ.l) THEN OPEN(IREST,FILE=FREST,FORM=Â’ UNFORMATTEDÂ’ ,STATUS=Â’OLDÂ’ ) CALL RESTART ( NX , NY , TS , SO , S , Y , T , TW ) CLOSE(IREST) PRINT* , Â’=== RESTART OK ! ======Â’ ELSE TS = 0. DO 100 J=1 ,NY DO 100 1=1, NX Y(I,J) = (J-1)*DY T( I , J) = TIK 100 CONTINUE DO 120 1=1, NX S0(I,1) = 0. S (1,1) = 0. S0(I ,2) = 0. S (1,2) = 0. 120 CONTINUE END IF CINPUT BOUNDARY CONDITIONS. READ ( 2 , * ) READ(2 , *) KBC , IBC -KBC>0: Uniform T. B.C. KBC<0: X-variant T. B.C. -IBC>0: Constant T. B.C. IBC<0: time-variant T. B.C. NBC = IABS(KBC) IF(NBC.GT.MAXBC) STOP Â’KBC is too largeÂ’ READ(2 , *) READ(2,*) (XBC( I ) , 1=1 ,NBC) READ(2,*) READ(2 , *) TIB2 READ(2,*) (TBC(I ,2) , 1=1 ,NBC) TIBI = TIB2 DO 300 1=1, NBC TBC( 1,1) = TBC(I ,2) 300 CONTINUE 310 CONTINUE IF(TS.LE.TIB2) GOTO 500 TIBI = TIB2 DO 320 1=1, NBC 320 TBC( 1 , 1 ) = TBC(1 ,2) READ(2 , * ,END=500) TIB2 READ(2 , *) (TBC(I ,2) , 1=1 ,NBC) GOTO 310

PAGE 172

o o 157 500 CONTINUE CALL I NTBC ( NX , NBC , DX , TS , TW ) RETURN END C10 20 ===== SUBROUTINE I NTBC ============ SUBROUTINE INTBC ( NX , NBC , DX , TIME , TW) PARAMETER (MAXBC=50) COMMON /BC/XBC(MAXBC) ,TBC(MAXBC,2) ,TIB1 ,TIB2 DIMENSION U(MAXBC) ,Y2(MAXBC) ,TB(MAXBC) ,TW(NX) DTI ME = TIB2-TIB1 DO 9 1=1, NBC IF(ABS(DTIME) . LT.0.01) THEN TB( I ) = 0 . 5*(TBC( I , 1 )+TBC( I ,2) ) ELSE TB( I ) = TBC(I,1) + (TBC(I ,2)-TBC( I , 1 ) )*(TIME-TIB1 )/DTIME END I F CONTINUE YP1 = 0. YPN = 0. Y2( 1 ) = -0.5 U(l) = (3 . / (XBC(2)-XBC( 1 ) ) )*( (TB(2) -TB( 1 ) )/ (XBC(2)-XBC(1))-YP1) I =2, NBC1 (XBC( I)-XBC( 1-1 ) )/ (XBC( 1+1 )-XBC( 1-1 ) ) SIG*Y2(I-l)+2. = (SIG-1 . ) /P = (6.*((TB(I+1)-TB(I))/(XBC(I+1)-XBC(I))(TB( I )-TB( 1-1 ) )/(XBC( I )-XBC( 1-1 ) ) )/ (XBC(I+1)-XBC(I-1))-SIG*U(I-1))/P CONTINUE QN = 0.5 UN = (3. / ( XBC ( NBC ) -XBC( NBC-1 )) )*(YPN-(TB( NBC) -TB( NBC-1 ) ) / ( XBC (NBC) -XBC ( NBC1 ) ) ) Y2(NBC) = (UN-qN*U(NBC-l))/(QN*Y 2 (NBC-l)+l. ) DO 20 I=NBC-1 ,1,-1 Y2( I ) = Y2( I )*Y2( 1+1 ) + U(I) CONTINUE DO 50 1=1, NX XX = (I-1)*DX IF(XX.LE.XBC(1)) THEN TWC = TB( 1 ) + 273.0 ELSE DO 30 J=2,NBC IF(XX.GE.XBC(J-1). AND . XX . LE . XBC ( J ) ) THEN H = XBC( J)-XBC( J-l ) A = (XBC( J)-XX)/H B = (XX-XBC( J-l ) ) /H TERM1 = A*TB( J-l )+B*TB( J) TERM2 = ( ( A**3-A)*Y2( J-l ) + (B*=t<3-B)*Y2( J) )=*(H*H) /6 TWC = TERM1+TERM2 + 273.0 DO 10 SIG = P Y2( I ) U(I)

PAGE 173

158 30 40 50 CC== GOTO 40 END IF CONTINUE TVC = TB(NBC) + 273.0 END IF TW( I ) = TVC CONTINUE RETURN END :=== SUBROUTINE RESTART ============ SUBROUTINE RESTART( NX , NY , TS , SO , S , Y , T , TV ) DIMENSION Y(NX ,NY) ,T(NX ,NY) , S(NX,2) ,S0(NX,2) ,TV(NX) COMMON /IO/I READ , I OUT , I REST , I DUMP READ(IREST) TS READ( I REST) ( (Y( I , J) ,T( I , J) , 1=1 ,NX) , J=1 ,NY) READ( I REST) (S( I , 1 ) ,S( I , 2) ,S0( I , 1 ) ,S0( I ,2) ,TV( I ) , 1=1 ,NX) C C: C C: C C RETURN END ===SUBROUT I NE DUMP===============— = SUBROUTINE DUMP (NX , NY , FDUMP , TS , SO , S , Y , T , TV) DIMENSION Y(NX,NY) ,T(NX,NY) , S(NX,2) ,S0(NX,2) ,TV(NX) COMMON / I 0/ IREAD , I OUT , IREST , I DUMP CHARACTERS FDUMP OPEN ( IDUMP , FILE=FDUMP , FORM= ’ UNFORMATTED ’ ) VRITE(IDUMP) TS VRITE(IDUMP) ((Y(I , J) ,T(I , J) ,1=1 ,NX) , J=1 ,NY) CLOSE(IDUMP) (S(I,1) ’ S(I ’ 2 >’ S0 ( I ’ 1 )’ S °( 1 ’ 2 )’ TV(I),I=1,NX) RETURN END ===== SUBROUTINE PROCES =============== Main solution procedure. FDM and SSM methods. CSUBROUTI NE PROCES ( NX , NY , NXY , KBC , IBC , ITFS , TSTART , ITOUT , TSTOP 1 DTT,HCOV,TINF,TMK,TIK,Y,T ,T0,S,S0,TV,TV0,TV1 DUM2) PARAMETER (MAXBC=50) ’ ' COMMON /BC/XBC(MAXBC) ,TBC(MAXBC,2) ,TIB1 ,TIB2 DIMENSION S(NX,2) ,S0(NX,2) ,TV(NX) ,TVO(NX) ,Y(NX,NY) ,T(NX,NY) ,TO(NX,NY) DIMENSION TV1(1), DUM2(1), ITFS(NX,2), TO(NX,NY) COMMON / ALL/XDIM , YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , HFG COMMON /SS/TSll ,TS12,TS21 ,TS22 NCY = 0 DX2 = DX*DX HCP = HFG/CP TIME = TSTART

PAGE 174

o o 159 AMP = 50. PRID = 30. NBC = IABS(KBC) DO 2 1=1, NX TV1(I) = TV(I) ITFS( I , 1 ) = 0 ITFS(I ,2) = 0 2 CONTINUE DO 900 WHILE (TIME. LT.TSTOP) NCY = NCY + 1 C-Slow start IF(DTT . LT. 0 . AND. NCY.LE.10) THEN DT = ABS(DTT)*(0. 1*NCY)**2 ELSE DT = ABS(DTT) END IF CTIME = TIME + DT IF(KBC.LT.O) CALL GETTV(NX, NBC, DX, TIME, TV) DO 100 1=1, NX IF(IBC.EQ.-l) TW( I ) = TV1(I) + 2 AMP*SIN(3 . 14159* (TIME-0 . 5*DT) /PRID) TWO(I) = T( I , 1 ) ’’ ’ DO 100 J=1 ,NY T0( I , J) = T( I , J) 100 CONTINUE C C-Solve for U-problem CALL FINDT1(NX,NY,NCY,HCOV,TINF,TIK,TO,T,TW, 1 T p, Anc , UI?r , D ! J ! 2(1) ’ DUM2 ( NXY+1 )’ DUM2 ( 2 * NXY+ l)»DUM2(3*NXY + l)) IF(ABS(HFG) . LT. 1 .E-8) GOTO 600 C C-Check for interface status CALL CHKS (NX , NY , ITFS , S , SO , TV , TwO , TMK , YDIM) C PRINT* PRINT* , ’Find interface at t=’,TIME PRINT*,’ S1(I) S2( I ) Twl(I)’ DO 500 1=1, NX K ’ IF( ITFS ( I > 1 ) • EQ . 0 . AND. ITFS( I , 2) .EQ . 0) GOTO 500 DO 200 J=1 ,NY J1 = NXY+J DUM2( J) = T(I , J) DUM2( Jl) = Y( I , J) 200 CONTINUE IF( ITFS( I , 1 ) .EQ . 1 . AND . ITFS( I , 2) . EQ . 0) THEN -SI. (Melting) HFGCP = HCP SS = S(I,1) SSO= S0(I,1) c

PAGE 175

o o o o C Search for single melting interface 160 C cc Cc c c cCALL FIND1S(NY,SS0,SS,HFGCP,TMK,DUM2(1) ,DUM2(NXY+1) 1 DUM2(2*NXY+1 ) ,DUM2(3*NXY+1) ,DUM2(4*NXY+1 ) , DUM2(5*NXY+1 ) ) -Solve for V-problem, update S and T CALL FINDT2(NY , TMK , HFGCP , DUM2( 1 ) , DUM2(NXY+1 ) , SSO , SS ) S (1,1) = SS S0(I,1) = SSO TS = TS11 + TS12 ELSE IF( ITFS( I , 1 ) . EQ . 0 . AND. ITFS( I , 2) . EQ . 1 ) THEN -S2. (Freezing) HFGCP = -HCP SS = S(I ,2) 550 =S0(I ,2) -Search for single freezing interface CALL FIND1S(NY, SSO, SS, HFGCP, TMK, DUM2(1) ,DUM2(NXY+1 ) 1 DUM2(2*NXY+1) , DUM2(3*NXY+1 ) , DUM2(4*NXY+1 ) , DUM2(5*NXY+1 ) ) -Solve for V-problem, update S and T CALL FI NDT2 ( NY , TMK , HFGCP , DUM2 ( 1 ) , DUM2 ( NXY+ 1 ) , SSO , SS ) S (1,2) = SS S0( I ,2) = SSO TS = TS11 + TS12 ELSE IF(ITFS(I,l).EQ.l.AND.ITFS(I,2).EQ.l) THEN SI L S2 . (Si-Melting, S2-Freezing) 551 = s ( 1 , 1 ) 5501 = S0(I,1) 552 = S (1,2) 5502 = S0( I ,2) HFGCP = HCP -Search for double interfaces CALL FIND2S (NY , SS01 , SS02 , SSI , SS2 , HFGCP , TMK , 1 DUM2( 1 ) , DUM2(NXY+1 ) , DUM2(2*NXY+1 ) , 2 DUM2(3*NXY+1 ) , DUM2(4*NXY+1 ) , DUM2(5*NXY+1 ) ) -Solve for V-problem, update S and T CALL FINDTT(NY , tmk , HFGCP , SS01 , SS02 , SSI , SS2 , 1 DUM2(1) ,dum2(nxy+l) ) s (1,1) = SSI S 0 (I, 1 ) = ssoi s (1,2) = SS2 S0( I , 2) = SS02 END IF C DO 300 J=1 ,NY T(I , J) = DUM2( J) 300 CONTINUE

PAGE 176

o n 500 161 WRITER, 1000) S(I,l),S(I,2),T(I,l),(i-i)*dx CONTINUE ' 600 700 900 1000 CONTINUE IP = MOD(NCY, ITOUT) DO 700 1=1, NX IF(IP.EQ.O) WRITE(30, 1000) CONTINUE (I-1)*DX,S(I,1),S(I,2),TIME IC = NX/2 WRITE(50 , 1000) TIME, WRITE(51 , 1000) TIME, CONTINUE FORMAT (5E1 5 .5) RETURN END S(1,1),S(IC,1),S(NX,1) S(1,2),S(IC,2),S(NX,2) C C===== SUBROUTINE GETTW =============== SUBROUTINE GETTW ( NX , NBC , DX , T I ME , TV ) PARAMETER (MAXBC=50) COMMON /BC/XBC(MAXBC) ,TBC(MAXBC , 2) ,TIB1 ,TIB2 DIMENSION TV (NX) C— 310 CONTINUE IF(TIME. LE.TIB2) GOTO 500 TIBI = TIB2 DO 320 1=1, NBC 320 TBC( I , 1 ) = TBC( 1,2) READ(2,*,END=400) TIB2 READ(2,*) (TBC( I ,2) , 1=1 ,NBC) GOTO 310 400 TIB2 = 1.E08 500 CONTINUE CALL INTBC(NX, NBC, DX, TIME, TW) RETURN END C= === SUBROUTINE CHKS============ SUBROUT I NE CHKS ( NX , NY , I TFS , S , SO , TW , TWO , TMK , YD I M ) DIMENSION I TFS (NX, 2) ,S(NX,2) ,S0(NX,2) , TWO (NX) ,tw(nx) DO 10 1=1, NX IF(Tw(I)-TMK.GT.0.5. a nd.TW0(I)-TMK. le.0.5) THEN C-Tw > Tm : Melting. ITFS(I , 1) = 1 r T . ELSE ! f ( tw ( i )-™ K lt -0.5.and.TW0(I)-TMK.ge.0.5) then t-lw < Im : Freezing. ITFS(I ,2) = 1 END IF C-Interface reaches other boundary IF(ITFS( I , 1 ) . EQ . 1 . AND . S( I , 1 ) . GE. YDIM) THEN ITFS( I , 1 ) = 0 S(I,1) = 0. END IF

PAGE 177

n o o o 162 IF( ITFS( I , 2) .EQ . 1 . AND.S( I ,2) . GE. YDIM) THEN ITFS(I ,2) = 0 S(I ,2) = 0. END IF C-Interface merges J^c^ , !x' GT ' 1 ' E ~ 4 ' AND ' ABS( s ^ I ’ 2 ^ _s ( I ’ 1 ^LE l E 8 ) then 1115(1 , 1 ) =0 ITFS(I ,2) = 0 S(I,1) = 0. S( I , 2) = 0. S0(I,1) = 0. S0(I ,2) = 0. END IF CIF 'cJt^ I n ‘ GT • ™ K * AND • S ( I , i ) .GT.S(I,2) .AND.S(I,2) .GT.O. ) THEN 5( I > 1 ) = 0 . S0(I,1) = o. S(I , 2 ) = o. SO(I ,2) = 0. ITFS( I , 1 ) = 1 ITFS( I ,2) = 0 END IF IF ^ t ( ^' LT, ™ k,and,s(i,2) gt s(i ’ 1 )and s ( i ’ 1 )gt -°0 then S(I ,2) =0. SO(I ,2) = 0. S(I,1) = 0. S0(I,1) = 0. ITFS( I , 1 ) = 0 ITFS(I ,2) = 1 END IF 10 CONTINUE RETURN END ==== SUBROUTINE FINDS ================== Find the interface position S at time = t C SUBROUTINE FIND1S(NY,S0,S,HFGCP,TMK,TI ,YI , AA7BBrCC, DD) DIMENSION YI(1),TI(1),AA(1),BB(1),CC(1),DD(1) ’ COMMON /ALL/XDIM , YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , HFG EPSS = l.E-06 EPST = l.E-05 DSDT = ABS(S-SO) SSO = S IF(DSDT.LE. 1 .E-08) THEN DDD = 0 . 025*YDIM ELSE DDD = DMAX1(DSDT, 0.0001 DO) END IF DO 100 J=1 ,NY AA( J) = 0.5

PAGE 178

BB( J) = 2.0 CC(J) = 0.5 163 100 200 10 CCONTINUE DD( 1 ) = 3.*(TI(2)-TI(1)) /DY DD(NY) = -3.*(TI(NY-1)-TI(NY))/DY DO 200 J=2,NY-1 DD(J) =3.*(AA( J)*(TI(J+1)-TI(J))+CC(J)*(TI(J)-TI(J-1)))/DY CONTINUE v v CALL TRID(NY, AA,BB,CC,DD) YL = SSO DDD RT = SSO + DDD IF(YL.LT.O. ) YL = 0 IF(RT.GT.YDIM) RT=YDIM CONTINUE YPRE = YL CALL SPLINE (NY, YI ,TI ,DD,YL,T1) CALL GREFT ( YL , SSO , YPRE , G I NT ) T2 = GINT*HFGCP FL = T1 + T2 TMK YLOC = RT YPRE = YLOC CALL SPLINE(NY,YI ,TI ,DD,RT,T1) CALL GREFT(YLOC, SSO, YPRE, GINT) T2 = GINT*HFGCP FH = T1 + T2 TMK IF(ABS(FH) . LT.EPST) GOTO 40 IF(FH*FL.GT.O. ) THEN YL = YL DDD RT = RT + DDD IF(YL.LT.O. ) YL = 0. IF(RT.GT.YDIM) then RT = YDIM GOTO 40 END IF END IF IF(FL.LT.O. ) THEN YH = RT FO = FL ELSE YH = YL YL = RT FO = FH END IF RT = 0.5*(YH+YL) YLOC = RT YPRE = YLOC DO = DDD D = DO CALL SPLINE(NY, YI ,TI ,DD,RT,T1 ) CALL GREFT(YLOC, SSO, YPRE, GINT) T2 = GINT+HFGCP

PAGE 179

o n o o F 164 = T1 + T2 TMK IF(ABS(F) .LT.EPST) GOTO 40 DF = (F-FL)/(RT-YL) RTO = RT FL = F IF(F.LT.O. ) THEN YL = RT ELSE YH = RT END IF CDO 20 L=1 ,200 IF(((RT-YH)*DF-F)*((RT-YL)*DF-F).GE.O. 1 . OR. ABS(2 . *F) . GT. ABS(DO*DF) ) THEN DO = D D = 0 . 5*(YH-YL) RT = 0.5*(YH+YL) ELSE DO = D D = F/DF RT = RT D END IF YLOC = RT YPRE = YLOC CALL SPLINE(NY,YI ,TI ,DD,RT,T1 ) CALL GREFT(YLOC,SSO,YPRE,GINT) T2 = GINT^HFGCP F = T1 + T2 TMK IF(ABS(F) .LT.EPST. OR. ABS(D) .LT.EPSS) GOTO 40 DF = (F-FL)/(RT-RTO) RTO = RT FL = F IF(F .LT.O. ) THEN YL = RT ELSE YH = RT END IF 20 CONTINUE 40 CONTINUE SO = S S = RT RETURN END == SUBROUTINE TRID =============== Tridiagonal solver SUBROUTINE TRID(M2,B,D,A,C1) DIMENSION A(1),B(1),D(1),C1(1) Ml = 1 MP1 = Ml+1 DO 100 I=MP1 ,M2

PAGE 180

o o o o 165 RR = B( I )/D( I-l ) D( I) = D(I)-RR*A(I-1) C1(I)= C1(I)-RR*C1(I-1) 100 CONTINUE RR = l./D(M2) Cl(M2)=Cl(M2)*RR DO 200 I=MP1 ,M2 J = M2-I+M1 RR = l./D(J) C1(J)=(C1(J)-A(J)*C1(J+1))*RR 200 CONTINUE RETURN END == SUBROUTINE SPLINE =========== Spline interpolation SUBROUTINE SPLINE(NY,Y,T,D,YY,TS1) DIMENSION Y(1),T(1),D(1) IF(YY.LE.Y(1)) THEN TS1 = T(l) RETURN END IF DO 10 J=2,NY IF(YY.GE.Y(J-1).AND.YY.LE.Y(J)) THEN DL = YY-Y(J-l) DR = Y( J)-YY DL2 = DL*DL DR2 = DR* DR H = Y(J)-Y(J-1) TERM1 = (3 . *H 2.*DR)*DR2*T(J-1)/H TERM2 = (3 . *H 2 . *DL)*DL2*T( J ) /H TERM3 = (H-DR)*DR2*D( J-l ) TERM4 = (H-DL)*DL2*D( J ) TS1 = (TERM1+TERM2+TERM3-TERM4) /H**2 GOTO 20 END IF 10 CONTINUE TS1 = T(NY) 20 CONTINUE RETURN END C C == =— ===== Subroutine GREFT =========== C Integration of GREEN’s function over t from t 0 to t SUBROUTINE GREFT(YLOC , YPO , YP1 ,GINT) COMMON /ALL/XDIM , YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , HFG EPS = l.E-04 IF(ABS(YLOC) .LT. 1 .E-08) THEN GINT = 0. RETURN

PAGE 181

o o o o 166 END IF PI = 3.1415927 DR = YP1 YPO DR2 = DR* DR IF(ABS(DR) . LT. 1 . E-08) THEN GINT = 0. RETURN END IF PIY = 0 . 5*PI/YDIM ALPHAT = ALPHA*DT SUM = 0. SUM1 = 0. DO 100 IT=1, 50000 BETA = (2*IT-1 )*PIY BRO = BETA* YPO BR1 = BETA*YP1 BXL = BETA*YLOC BETA 2 = BETA* BETA ABT = ALPHAT*BETA2 EPT = EXP(-ABT) DEN = ABT*ABT + BETA2*DR2 T1 = ABT/DEN T2 = BETA*DR/DEN TERM1 = T1*SIN(BR1 ) TERM2 = T2*C0S(BR1 ) TERM3 = T1*SIN(BR0) TERM4 = T2*COS(BRO) SUMI = TERM 1 -TERM2 -EPT * ( TERM3-TERM4 ) SUM = SUM + SUMI*SIN(BXL) SUMI = SUM1+SUMI TT = T1+T2 ERR = ABS(TT/SUM1 ) IF(ERR.LT.EPS) GOTO 200 100 CONTINUE IF(ERR.GE. 1 .E-3) THEN WRITE(*,*) ’No convergence !’ STOP No convergence in GreenS function’ END IF 200 CONTINUE GINT = 2 . *DR*SUM/YDIM RETURN END ===== SUBROUTINE FINDT2 ============ Find temp, correction due to interface SUBROUTINE FINDT2(NY,TMK,HFGCP,TI ,YI ,S0 S') DIMENSION YI (NY) ,TI (NY) REWIND ( 40 ) /XDIM ’ YDIM ’ M ’ ° Y ’ ° T ’ CK ’ CP ’ RH0 ’ ALPHA ’ HFG IF(S . GE. YDIM . OR. S. LE. 0 . ) RETURN YPRE = S

PAGE 182

o n o o sso = so DO 400 J=1 ,NY 167 YLOC = YI( J) CALL GREFT ( YLOC , SSO , YPRE , G I NT1 ) T2 = GINT1*HFGCP WRITE(40 , 10) YLOC , TI ( J) -TMK , T2 , T2+TI ( J ) -TMK TI( J) = TI(J) + T2 400 CONTINUE 10 F0RMAT(4F14.4) RETURN END ===== SUBROUTINE FINDT1 = FDM method for solving Tl SUBROUTINE FINDTl(NX,NY,ncy,HCOV,TINF,TIK,TO T TW 1 A, C, B, D ) DIMENSION TO(NX,NY) ,T(NX,NY) ,TW(NX) DIMENSION A(l) ,B(1) ,C(1) ,D(1) COMMON /ALL/XDIM , YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , HFG DX2 = DX*DX DY2 = DY*DY BIO = DX*HCOV/CK IF(HCOV . LE. 0. ) BIO = 0. COEX = 0 . 5*ALPHA*DT/DX2 COEY = 0 . 5*ALPHA*DT/DY2 N1=NY-1 C-SVEEP IN J-DIR. DO 300 1=1, NX DO 100 J=1 ,N1 JJ=J+1 A( J) = -COEY B(J) = 1.+ 2 . *COEY C(J) = -COEY IF(NX.EQ.l) THEN D2DTX2 = 0. ELSE IF(I.EQ.l) THEN D2TDX2 = 2.*C0EX*(T(I+1, JJ)-T( I , JJ)-BIO*(T( I , JJ)-TINF) ) ELSE IF(I.EQ.NX) THEN J D2TDX2 =-2.*C0EX*(T(I , JJ)-T(I-1 , JJ)+BIO*(T(I , JJ)-TINF)) ELSE D2TDX2 = C0EX*(T(I+1 , JJ)-2.*T(I,JJ)+T(I-1,JJ)) END IF V Â’Â’ END IF D( J) = T( I , JJ) + D2TDX2 100 CONTINUE T0( I , 1 ) = TW( I ) D( 1 ) = D( 1 ) + COEY*TW( I ) A(N1) = -2*C0EY CALL TRID(N1,A,B,C,D) DO 200 J=1 ,N1

PAGE 183

JJ = J+l T0( I , JJ) = D(J) 168 200 CONTINUE 300 CONTINUE C C-SVEEP IN I-D J1 = 2 DO 600 J=J1 ,NY DO 400 1=1, NX A( I) = -COEX B(I ) = 1.+ 2. *COEX C(I) = -COEX IF(J.EQ.NY) THEN D2TDY2 = 2.*C0EY*(T0(I , J-1)-T0(I , J)) ELSE END°IF = C0EY *C TO C 1 , J -*-i ) 2 . * T °( 1 , J)+TO( I , J-l ) ) D( I) = TO(I , J) + D2TDY2 400 CONTINUE IF(NX.GT.l) THEN C ( 1 ) = -2 . *COEX B(l) = 1.+ 2.*C0EX*(1.+BI0) D(l) = D( 1 ) + 2.*C0EX*BI0*TINF A(NX) = -2.*C0EX B(NX) = 1 . +2.*C0EX*(1 ,+BIO) D(NX) = D(NX) + 2. *COEX*BIO*TINF CALL TRID(NX,A,B,C,D) END IF DO 500 1=1, NX T(I,J) = D( I ) 500 CONTINUE 600 CONTINUE DO 700 1=1, NX T( I , 1 ) = TW( I ) 700 CONTINUE RETURN END C C==== SUBROUTINE OUTPUT ============ C C= SUBROUTINE OUTPUT ( NX , NY , TSTOP , Y , T ) DIMENSION Y(NX,NY) ,T(NX,NY) COMMON / ALL/XDIM, YDIM , DX , DY , DT , CK , CP COMMON /IO/IREAD, IOUT, IREST, IDUMP ,RHO, ALPHA, HFG 100 1000 WRITE( IOUT , *) NX,NY,0 DO 100 J=1 ,NY DO 100 1=1, NX X = ( I-1)*DX WRITE( IOUT, 1000) X , Y( I , J) ,T( I , J) CONTINUE F0RMAT(5E16 .6)

PAGE 184

o o o a 169 RETURN END ;=== SUBROUTINE FIND2S ============= Search for double interfaces SUBROUTINE FIND2S(NY,S01 ,S02 ,S1 ,S2,HFGCP,TMK, 1 TI,YI, AA, BB, CC, DD ) DIMENSION YI(1),TI(1),AA(1),BB(1),CC(1),DD(1) COMMON /ALL /XDIM , YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , HFG C= EPSS = l.E-06 EPST = l.E-05 551 = SI 552 = S2 DSDT1 = Sl-SOl DSDT2 = S2-S02 DIRO = S01-S02 IF(ABS(DSDTl).LE.l.E-08) THEN DDD1 = 0 . 005*YDIM IF(S01.LE.l.E-8) DDD1=0 . 1*YDIM ELSE DDD1 = DSDT1 END IF IF(ABS(DSDT2).LE.l.E-08) THEN DDD2 = 0 . 005*YDIM IF(S02 . LE . 1 . E-8) DDD2=0 . 1*YDIM ELSE DDD2 = DSDT2 END IF DO 100 J=1 ,NY AA( J) = 0.5 BB( J) = 2.0 CC(J) = 0.5 100 CONTINUE DD( 1 ) = 3.*(TI(2)-TI(1))/DY DD(NY) = 3.*(TI(NY)-TI(NY-1)) /DY DO 200 J=2 ,NY-1 DD( J) =3.*(AA(J)*(TI(J+1)-TI(J))+CC(J)*(TI(J)-TI(J-1)))/DY 200 CONTINUE CALL TRID(NY,AA,BB,CC,DD) YL1 = SSI RT1 = SSI DDD1 YL2 = SS2 RT2 = SS2 DDD2 RT1 = DMAX1(YI(1) ,RT1) RT2 = DMAX1(YI(1),RT2) DIR1 = RT1-RT2 IF(DIR0*DIR1 . LT . 0 . ) THEN IF(DIRO.GT.O. ) RT1 = RT2 RT2 = RT1 GOTO 40 END IF

PAGE 185

170 YPRE1 = RT1 YPRE2 = RT2 YLOC = RT1 CALL SPLINE(NY,YI,TI,DD,YLoc,Tll) CALL GREFT2 ( YLoc , SSI , ss2 , YPRE1 , ypr e2 , G INTI , g i nt2 ) T21 = HFGCP* (gintl gint2) FL1 = Til + T21 TMK YLOC = RT2 CALL SPLINE(NY, YI ,TI ,DD,Yloc,T12) CALL GREFT2(YLoc , ssl ,SS2,YPRE1 ,ypre2, gintl ,GINT2) T22 = HFGCP* (gintl gint2) FL2 = T12 + T22 TMK RT1 = SSI + DDD1 RT2 = SS2 + DDD2 IC0NT=O. CONTINUE ICONT = ICONT + 1 DIR1 = RT1-RT2 IF(DIR0*DIR1 . LT.O. ) THEN IF(DIRO.GT.O-) RT1 = RT2 RT2 = RT1 GOTO 40 END IF YLOC = RT1 YPRE1 = RT1 YPRE2 = RT2 CALL SPLINE( NY , YI , TI , DD , y loc , T1 1 ) CALL GREFT2 ( YLoc , SSI , ss2 , YPRE1 , ypre2 , G INTI , g i nt2 ) T21 = HFGCP* (gintl g int2) FH1 = Til + T21 TMK YLOC = RT2 CALL SPLINE(NY,YI ,TI ,DD,yloc,T12) CALL GREFT2 ( YLoc , SSI , ss2 , YPRE1 , ypre2 , G INTI , g i nt2 ) T22 = HFGCP* (gintl gint2) FH2 = T12 + T22 TMK IF(FH2*FL2.GT.O. ) THEN IF(ICONT.GT. 10) THEN DDD2 = -10*DDD2 ICONT = 0 END IF IF(RT2.GE.YDIM) DDD2 = -10.*DDD2 IF(RT2.LE.YI(1)) DDD2 = -10.*DDD2 IF(ABS(DDD2) . GT . 0 . 2*YDIM) DDD2=0 . 5*DDD2 RT2 = RT2 + DDD2 IF(RT2.GT.YDIM) RT2 = YDIM IF(RT2 . LT. YI (1 ) ) RT2 = YI(1) END IF IF(ABS(YLl-RTl).LT.l.E-8) RT1 = YL1 + l.E-05 if (ABS(YL2-RT2) . It. 1 .E-8) RT2 = YL2 + l.E-05 DO 20 L=1 ,50

PAGE 186

171 all = (fhl-f 1 1 )/( r tl-y 1 1 ) al2 = (fhl-f 11 )/( r t2-yl2) a21 = (fh2-f 12)/ (rtl-yll) a22 = (fh2-f 12)/ (rt2-yl2) if (abs(all) . It. 1 .e-8) then x2 = -f hl/al2 xl = ( f h2+a22*x2) /a21 else hj = a22-a21*al 2/all if (abs(h j) . It . 1 . e-8) then xl = -fhl/all x2 = -fh2/a22 else x2 = (a21*fhl/all fh2)/(a22 al2*a21/all) xl = -(fhl+al2*x2)/all end i f end if if (abs(xl)+abs(x2) . le.epss) goto 40 IF(ABS(X1 ) . GE. EPSS) THEN yll = rtl FL1 = FH1 END IF IF(ABS(X2) .GE.EPSS) THEN yl2 = rt2 f 12 = fh2 END IF if (abs(xl) .gt.0.5*ydim) xl=0 . l*xl*ydim/abs(xl ) if (abs(x2) .gt.0.5*ydim) x2=0 . I*x2*ydim/abs(x2) rtl = rtl + xl ' rt2 = rt2 + x2 YLOC = HT1 YPRE1 = RT1 ypre2 = rt2 CALL SPLINE ( NY, YI,TI,DD,yloc, Til) CALL GREFT2(YLoc,SSl ,ss2,YPREl ,ypre2,GINTl ,gint2) T21 = HFGCP* (gintl gint2) FH1 = Til + T21 TMK YLOC = RT2 CALL SPLINE(NY,YI,TI,DD,yloc,T12) CALL GREFT2(YLoc ,SS1 , ss2 , YPRE1 ,ypre2,GINTl ,gint2) T22 = HFGCP* (gintl gint2) FH2 = T12 + T22 TMK 20 CONTINUE 40 CONTINUE TS11 = Til TS12 = T12 TS21 = T21 TS22 = T22 SOI = SI s02 = s2 SI = RT1 s2 = rt2

PAGE 187

n o 172 if (ABS(sl-s2) .LT. 1 .E-6) then si = 0 . 5*( sl+s2) s2 = si end if RETURN END === GREENÂ’S FUNCTION FOR FINITE DOMAIN ====== SUBROUTINE GREFT2(YLOC,Y01 ,y02,Ypl ,yp2,Gl ,g2) COMMON /ALL/XDIM , YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , HFG EPS = l.E-04 IF(ABS(YLOC) . LT. 1 .E-08) THEN G1 = 0. g2 = 0. RETURN END IF PI = 3.1415927 DR1 = YP1 Y01 DR2 = YP2 Y02 DR12 = DR1*DR1 dr22 = dr2*dr2 PIY = 0 . 5*PI/YDIM ALPHAT = ALPHA*DT SUMgl = 0. sumg2 = 0. SUM1 = 0. sum2 = 0. DO 100 IT=1, 50000 BETA = (2*IT-1 )*PIY BR01 = BETA*Y01 BR11 = BETA*YP1 BR02 = BETA*Y02 BR12 = BETA*YP2 BXL = BETA*YLOC BETA2 = BETA*BETA ABT = ALPHAT*BETA2 EPT = EXP(-ABT) DENI = ABT*ABT + BETA2*DR12 DEN 2 = ABT*ABT + BETA2*DR22 Til = ABT/DEN1 tl2 = abt/den2 T21 = BETA* DR 1 /DENI T22 = BETA*DR2 /DEN2 TERM11 = T11*SIN(BR11) TERM21 = T21*C0S(BR11 ) TERM31 = T11*SIN(BR01) TERM41 = T21*C0S(BR01 ) TERM 12 = T12*SIN(BR12) TERM22 = T22*C0S(BR12) TERM32 = T12*SIN(BR02)

PAGE 188

o o o o 173 TERM42 = T22*C0S(BR02) SUM II = TERM11-TERM21-EPT*(TERM31-TERM41 ) SUMgl = SUMgl + SUMI1*SIN(BXL) SUM I 2 = TERM12-TERM22-EPT*(TERM32-TERM42) SUMg2 = SUMg2 + SUMI2*SIN(BXL) SUM1 = SUM1+SUMI1 sum2 = sum2+sumi2 TT1 = T11+T21 TT2 = T12+T22 IF( IT . LT . 20) GOTO 100 ERR = (ABS(TT1 )+ABS(TT2) )/ (SUM1+SUM2) IF(ERR.LT.EPS) GOTO 200 100 CONTINUE IF(ERR.GT. 1 .E-3) THEN WRITE(*,*) Â’No convergence !Â’ END IF 200 CONTINUE G1 = 2 . * DR 1* SUMgl /YDIM G2 = 2 . *DR2*SUMg2/YDIM RETURN END == =====SUBROUTINE FINDTT======= === Solve for V-problem for double interfaces 400 10 SUBROUTINE FINDTT(NY , TMK , HFGCP , SOI , S02 , SI , S2 , TI , YI ) DIMENSION TI (1 ) , YI ( 1 ) Â’ Â’ Â’ ; COMMON /ALL/XDIM, YDIM , DX , DY , DT , CK , CP , RHO , ALPHA , HFG COMMON /SS/TS11 ,TS12,TS21 ,TS22 REWIND (40) YPRE1 = SI SSI = SOI YPRE2 = S2 SS2 = S02 DO 400 J=1 ,NY YLOC = YI( J) CALL GREFT2(YL0C,SS1 ,ss2,YPREl,ypre2,GINTl,gint2) T2 = HFGCP* (gintl gint2) write(40 , 10) yloc,tI( j)-tmk,t2,t2+TI( j)-tmk TI(J) = TI( J) + T2 CONTINUE F0RMAT(4F14 .4) RETURN END C========== ==== END OF PROGRAM

PAGE 189

REFERENCES Alexiades , V, and Solomon, A. D. , 1993, Mathematical Melting and Freezing Processes, Hemisphere Publishing Washington, DC. Modeling of Corporation, Allen, D. N. de G. , relaxation methods differential equations two dimensions,” Quart. and Severn, R. T. , 1962, "The application of to the solution of non-elliptic partial ; III. Heat conduction, with change of state, in J. Mech. Appl. Math., Vol. 15, pp. 53-62. Anderson, D. A., Tannehill, J. C. , and Pletcher, R. Computational Fluid Mechanics and Heat Transfer, Hemisphere Corporation, Washington, DC. H., 1984, Publishing Barozzi, G. S. , and heat transfer probl pipe,” Trans. ASME, Pagl iarini , G., 1985, ”A method to solve conjugate ems: The case of fully developed laminar flow in a J. Heat Transfer, Vol. 107, pp. 77-83. Basu, B. and Date W. A., 1990, "Numerical study of steady state and transient laser melting problems -I. Characteristics of flow field and ea ransfer, Int. J. Heat Mass Transfer, Vol. 33, pp. 1149-1163. annl ! ” A refinement of the heat balance integral method 1357-1362° a mel S P roblem ’” InL JHeat Mass Transfer, Vol. 21, pp. Bell, G. E. , 1979, "Solidif icati on pipe, Int. J. Heat Mass Transfer, Vol. of a liquid about a cylindrical 22, pp. 1681-1686. B ' • A ' ’ 19?4 ’ ” The embedding technique solidification problems,” Moving Boundary Problems Diffusion, Proc. Conf . Univ. Oxford, March 25-27, pp. in melting and in Heat Flow and 150-172. Bonnerot, R. , and Jamet, P., 1974, ”A second order finite method for one-dimensional Stefan problem,” Int. J Numer. Meth. element Engng. , Boucheron E A., and Smith, R. N. , 1985, ”An enthalpy formulation of r® ,® lm P 1 ® ^orithrn for Phase change heat transfer problems,” ASME / AIChE Natl. Heat Transfer Conf., 23rd. Denver, August 4-7, 85-HT-7. 174

PAGE 190

175 Budhia, freezing H, and Kreith, F. , 1973, ’’Heat transfer with melting and in a wedge,” Int. J. Heat Mass Transfer, Vol. 16, pp. 195-211. Cars law, H. S. , and Jaeger, J. C. , Oxford University Press, Oxford, U. 1959, Conduction of Heat in Solids, K. Choi, C. melting variant Universi Y v ’ ”F Xac ^ and numeric al solution of oneand two-phase and solidification problems imposed with constant or timetemperature and flux conditions,” Ph.D. Dissertation, ty of Florida. Choi, C. Y and Hsieh, C. K. , 1992, ’’Solution of Stefan imposed with cyclic temperature and flux boundary conditions Heat Mass Transfer, Vol. 35, pp. 1181-1195. problems ” Int. J. Chou, F. C., Lien, V. Y. , and Lin, S experiment of non-Darcian convection in sphere channels— 1. Forced convection,” Vol. 35, pp. 195-205. H. , 1992, ’’Analysis and horizontal square packedInt. J. Heat Mass Transfer, Chuang, Y. K. , and Szekely, J., 1971, ”0n the use for solving melting or solidification problems Transfer, Vol. 14, pp. 1285-1294. of Green’s functions Int. J. Heat Mass Crank, J., 1981, ”How to deal with moving problems,” Numerical Method in Heat Transfer, K., and Zienkiewicz (eds.), Wiley, New York, pp boundaries Lewis, R. 177-200. in V., thermal Morgan , Crank, J. , Oxford. 1984, ’’Free and Moving Boundary Problems,” Clarendon Press, Duda, J. L., and Vrentas, J. S. , 1969, diffusion-controlled moving boundary probl 24, pp. 461-470. ’’Perturbation solutions of ems,” Chem. Eng. Sci., Vol. Duraunkaya, Z. , and Nair, S. , finite domain,” J. Appl. Mech. , 1990, ”A moving boundary problem in a Vol 57, pp. 50-56. Elliot, C. M., and Ockendon, for Moving Boundary Problems Advanced Publishing Program, J* R. , 1982, Weak and Variational , Research Notes in Mathematics 59, Boston . Methods Pitman Furzeland, R. M, 1980, ”A comparative study moving boundary problems”, J. Inst. Math. Appl., of numerical methods Vol. 26, pp. 411-429. for Gniel inski , V. ,1976, ’’New equations for heat and mass turbulent pipe and channel flow,” Int. Chem. Eng., Vol. transfer in 16, pp. 359Goodman, T.R., 1958, problems involving a pp. 335-341. The heat-balance integral and its application to change of phase,” ASME J. Heat Transfer, Vol. 80,

PAGE 191

176 Goodman, T. R. , 1964, ’’Application of integral nonlinear heat transfer,” Adv. Heat Transfer, Vol. methods to transient 1, pp. 71-79. T f'M a r d s i e i’ J J ” 1960 ’ ” The meitin « ° f finite «^bs,” AS ME J. Appl. Mech. , Vol. 27, pp. 16-24. Goodrich, L. E., 1978, ’’Efficient numerical dimensional thermal problems with phase change Transfer, Vol. 21, pp. 615-621. technique ” Int. J. for oneHeat Mass Gupta, R. S., and Kumar, method for one-dimensional Eng., Vol. 23, pp. 101-109. D. , 1980, ”A modified variable time step Stefan problem”, Comp. Meth. Appl. Mech. Gupta, R. S., and Kumar, D. , one-dimensional Stefan problem Heat Mass Transfer Vol. 24, pp. 1981, ’’Variable time step methods for with mixed boundary condition,” Int J 251-259. Hastaoglu, M. A., 1986, ”A numerical solution to problem Application to melting and solidification,” Transfer, Vol. 29, pp. 495-499. moving boundary Int. J. Heat Mass Hastaoglu, M A., and Baah , C. A., 1991, ’’Desublimation: A boundary problem and numerical solution,” Num. Heat Transfer, moving Part A , Hill, J. M., 1987, One-dimensional Stefan Problems: Longman Scientific and Technical, Essex, England. An Introduction, Hong, C. P. , Umeda , T., and Kimura, Y. , 1984, ” casting solidification: Part I, the coupling of and finite difference methods for solidification Transaction B, Vol. 15B, pp. 91-99. Numerical methods for the boundary element problems,” Metallurgical Hsieh , C. K. , 1993, ’’Unification boundary element for the solution Transfer, in press. of source-and-sink method and of potential problems,” J . Heat Hsieh, C. K., Choi, C. Y. , 1992a, ’’Solution of melting and solidification problems imposed with variant temperature and flux boundary conditions, Vol. 114, pp. 524-528. oneand constant ” J. Heat two-phase or timeTransfer , Hsieh, C. K., Choi, C. Y. , 19926, energy storage for solar energy Vol. 114, pp. 203-276. A general analysis of phase change applications,” J. Solar Energy Eng., Hsieh, C. K., Choi, C. Y. , and Stefan problems by a boundary Technology VII, C. A. Brebbia and Mechanics Publications, Southampton Kassab, A. J., 1992, ’’Solution of element method, Boundary Element M. S. Ingber , eds., Computational , UK, pp. 473-490.

PAGE 192

177 Hsu, C. F., Sparrow, E. M. , and Patankar, S. solution of moving boundary problems by boundary control-volume-based finite-difference scheme,” Transfer, Vol. 24, pp. 1335-1345. V., 1981, ’’Numerical immobilization and a Int. J. Heat Mass Kays, V. M. , and Crawford, M. E. , 1980, Convective Transfer, McGraw-Hill Publishing Company, New York. Heat and Mass Kolodner , I . I . , application to Mathematics, Vol. 1956, ’’Free boundary problems for heat equation with problems of change phase,” Comm. Pure Applied 9, pp. 1-31. Lacroix, M. , and Voller, V. R. , 1990, ’’Finite di solidification phase change problems: Transformed Numerical Heat Transfer, B, Vol. 17, pp. 25-41. f f erent versus solutions of fixed grids,” Lane, A. George, 1986, ’’Solar Heat Storage: Latent Vol. II: Technology, CRC press, Boca Raton, Florida. Heat Material,” Langford, D. , 1973, ’’The heat balance Mass Transfer, Vol. 16, pp. 2424-2428. integral method,” Int. J. Heat Liang, S. J., involving thin 1993, ”A study of free and moving boundary problems crystal growth,” Ph.D. Dissertation, Univ. of Florida. Lightfoot, N. M. H., 1929, ’’The solidification of molten London Math. Soc., Vol. 31, pp. 97-116. steel,” Proc. Murray, V. D. , and Landis, F. , 1959, of transient heat-conduction problems Part I -Method of analysis and sample 81, pp. 106-112. Numerical and machine solution involving melting and freezing; solutions,” J. Heat Transfer, Vol. Noble, B. , 1975, ’’Heat balance method in melting problems, boundary problems in heat flow and diffusion, Ockendon, J. Hodgkins, W. R. eds., Clarendon Press, Oxford, pp. 208-209. ” Moving R. , and Nyros, J. G., and Hsieh, C. K. , method and complextemperature problems imposed with cyclic Submitted to J . Heat Transfer. 1993, Combination of source-and-sink method for the solution of Stefan temperature and flux conditions,” Oberkampf, V. L. , 1976, ’’Domain mapping partial differential equations,” Int. J pp. 211-223. for the numerical solution of Numer. Meth. Engng . , Vol. 10, Ockendon, J. R. , 1974, ’’Moving boundary problems diffusion,” Proc. Conf. Univ. Oxford, March 25-27. in heat flow and O’Neill, K., 1983, ’’Boundary boundary phase change problems pp. 1825-1850. integral equation solution of moving Int. J. Numer. Meth. Engng., Vol. 19,

PAGE 193

178 Pa^el » P. D. , 1968 , ’’Interface condition in heat conduction with change of phase,” .47.4,4 J., Vol. 6, pp. 2454-2456. problems Pham, q. T., 1985, ”A fast, unconditionally stable finite scheme for heat conduction with phase change,” Int J Transfer , Vol. 28, pp. 2079-2084. difference Heat Mass Pham, Q. T., 1986, ”The use of lumped capacitance element solution of heat condition with phase change ” Transfer , Vol. 29, pp. 285-291. in the finiteInt. J. Heat Mass Poots, G., 1962, ”An approximate treatment of heat involving a two-dimensional solidification front ” Transfer, Vol. 5, pp. 339-348. conduction problem Int. J. Heat Mass Rathjen, K. or freezing A., and Jiji, L. M., 1971, ”Heat conduction with melting in a corner,” J. Heat Transfer, Vol. 93, pp. 101-109. Riaz, M. 1978, ’’Transient analysis of packed-bed systems,” Solar Energy, Vol. 21, pp. 123-128. thermal storage Rubinstein, L. I., 1971, Math. Monogr., No. 27. ’’The Stefan problems,” Am. Math. Soc. Trans. Saitoh, T., 1978, ’’Numerical method for multidimensional freezing problems in arbitrary domains,” J. Heat Transfer, Vol. 100, pp. 294-299 Sadegh, A. M., Jiji, L. M. , and Veinbaum, S., 1987 , ’’Boundary integral In^J 1 ^ / G M niq ^ e W1 / h application to freezing around a buried pipe,” Int. J. Heat Mass Transfer, Vol. 30, pp. 223-232. Sarler , B. , Alujevic, A., and Kuhn, multidimensional Stefan problems by BEM,” Mathematik und Mechanik, Vol. 71, pp. 611-614 , . G., 1991 , ’’Solution of Zeitschrift fur Angewandte Schumann, T. E. V. , 1929, ’’Heat transfer: a porous prism,” J . Franklin Inst., Vol. 208, pp. liquid flowing 405-416 through a !*•» Pang ’ Y> ’ Hunter, G. B. , Wei, D. Y. , and Che Modeling of turbulent transport during continuous ingot J. Heat Mass Transfer, Vol. 35, pp. 1229-1245. , M . -H. , 1992 , casting,” Int. Sikarskie, D. L., and Boley, B. A., 1965, ’’The solution two-dimensional melting and solidification problems,’ Struct., Vol. 1, pp . 207-234. of a class of Int. J. Solid Solomon, A. D. , 1978, Meger 1 in’s method for Moving Boundary Problems P.T. eds., Academic Press The application and extendabi 1 ity of the solving parabolic free boundary problems,” , Wilson, D. G., Solomon, A. D., and Boggs, , New York, pp. 187-202.

PAGE 194

Tamma , K . K , perspectives applications Vol. 30, pp. 179 and Namburu, R. R. , 1990, "Recent advance, trends and new via enthalpy-based finite element formulations for S03S 8°20 i
PAGE 195

BIOGRAPHICAL SKETCH Hongjun Li was born in Beijing, China, on September 15, 1961, one hundred thirty years after the phase-change problems were first studied. He received his B.S. degree in 1982 and M.S. degree in 1985 in Aerospace Engineering from Northwestern Polytechincal University, Xian, China. He finished his Ph.D. requirements in 1988 and was assigned to fulfill an international cooperation plan in W. Germany before he could make the final step — oral defense — to get the cap. After completing the cooperation contract, he forfeited his Ph.D. degree in China and came to thq United States in 1989, where he began his new education career starting from pre-elementary education: English, in the University of Florida. When he was 30-year old, he received his second M.S. degree in Mechanical Engineering from the University of Florida. In the same year, he had a lovely daughter Angela Y. Li. Then he was admitted to Ph.D. candidacy in 1992. H. Li married Yu Luo in 1986 in Xian, China. 180

PAGE 196

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Chung K. Hsi4h, Chairman Professor of Mechanical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. I I I * I ' Ot Tr, m -1 Yogi D.Goswami, Cochairman Professor of Mechanical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. We Professor of Aerospace Engineering, Mechanics and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. James F. Klausner Assistant Professor of Mechanical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Renwei Mei Assistant Professor of Aerospace Engineering, Mechanics and Engineering Science

PAGE 197

This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy . December 1993 Dean, College of Engineering Karen A. Holbrook Dean, Graduate School