SOURCEANDSINK 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 MultiDimensional 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 PHASECHANGE 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 Longtime 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 zdirection .......................... 58
3.4 System analyzed in cylindrical coordinates: Interface moves in rdirection .......................... 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 timevariant temperature condition .................................. 102
viii
5.5 Interface positions for timevariant 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 meltingfreezing 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
coalescingfront 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 Postprocessor for a two dimensional TES unit ............ 142
5.32 Postprocessor 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 ydirection 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 xdirection 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
SOURCEANDSINK 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 sourceandsink 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
Phasechange problems involving melting and freezing have found numerous applications in science and engineering. Phasechange 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 volumetoenergy 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.
Phasechange 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 phasechange 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 phasechange problems thus pose a real challenge for analysis.
Phasechange problems that have been solved in the literature have often been treated with various degrees of simplification. Earliest phasechange 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 timevariant 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.
Phasechange 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 phasechange temperature is not distinct. There is a crystallinelike 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 phasechange 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.
Phasechange 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 finitedifference method in one direction and a sourceandsink 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 phasechange interface
4
moves rapidly. The sourceandsink 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 sourceandsink 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 NeumannStefan formulation. The work involved in the solution can thus be reduced considerably for the sourceandsink method. In addition, in the traditional method of solution, the interface Stefan condition must be given as a separate condition, whereas in the sourceandsink method, this condition has been incorporated into the governing heat diffusion equation. The sourceandsink method is thus highly efficient in the solution of the phasechange problems.
In the present work, the sourceandsink 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 longtime solutions. An apparatus has also been built for experimentally testing its accuracy for shorttime 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 phasechange 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 phasechange 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 phasechange 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, diffusiondriven, phasechange solution scheme that combines an analytical solution in one direction and a numerical solution in another direction, only mathematical methods for diffusiondriven 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 phasechange 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 PhaseChange Problems
Analytical Solutions
Numerical Solutions
Exact Solutions (Neumann Solution)
Approximate
Analytical Solutions
 Power Series Solution Methods
 Solution of Integrodifferential Equations
(SourceandSink 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, phasechange 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 phasechange 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 shorttime 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.
Phasechange 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 sourceandsink method. Chuang and Szekely (1971) have used this method in the development of a numerical scheme for the solution of Stefan problems. Recently, the sourceandsink 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 sourceandsink 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 sourceandsink 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 sourceandsink 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 sourceandsink 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 timevariant 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 phasechange temperature (i.e., onephase problem) and imposed with a timeinvariant 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 higherorder 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 RungeKutta numerical integration algorithm.
2.1.3 Other approximate solutions
Stefan problems can also be solved by using quasisteady and quasistationary approximations. In the quasisteady approximation, the unsteady terms in the governing equations are dropped. This
13
results in a temperature that is good for longtime solution. As for the quasistationary 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 shorttime solutions.
In the asymptotic method of solution, these quasisteady and quasistationary limits are used for inner and outer expansions in the construction of temperature and interface position. Duda and Vrentas (1969) have demonstrated that the quasistationary approximation is actually the zeroth order solution of a regular perturbation, which is, indeed, good for a shorttime solution. On the other hand, Weinbaum and Jiji (1977) have shown that the lowest order term of a longtime perturbation actually leads to a quasisteady limit and is thus good for a longtime 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 firstorder 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.
Phasechange problems have also been solved by a Fourier transform (Ockendon, 1974). In this method, the Stefan problems are reduced to an initialvalue problem. In a similar fashion, a Laplace transform can be used to change the problems into boundaryvalue
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 phasechange 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 phasechange 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 timevariant. 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 thirdorderaccurate finite element method with an adaptive grid. In their work, the computational grid is reconstructed 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 fronttracking methods and have been used extensively in the numerical solution of phasechange 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 frontfixing or boundaryfixing. Furzeland (1980) and Lacroix and Voller (1990) have made extensive studies for comparison of the fixedgrid method and the transformedgrid 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 bodyfitted 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 energyconserving finitedifference scheme in the analysis of freezing outside a cooled tube. Shyy et al. (1992) used bodyfitted 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 RankineHugoniot equation in the shockwave analysis, the weak formulation can be used to solve phasechange 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 phasechange 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 phasechange temperature, the nonlinearity behavior of the enthalpy near the phasechange 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
Phasechange 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 timedependent 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 phasechange 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 timestep limitations. Sadegh et al. (1987) applied BEM to the solution of solidification around a buried pipe in a semiinfinite 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 phasechange 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 sourceandsink method and their unification is discussed in a great detail.
CHAPTER 3
NUMERICAL ANALYSIS
In this chapter, the solution of Stefan problems by a sourceandsink 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 phasechange 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 liquidphase 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 outwarddrawn 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 sourceandsink 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(xS(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)dr (3.10) rt
0t
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 semiinfinite domain. In the problems they solved, the medium was initially at a uniform temperature. The boundaries were imposed with constant or timevariant temperature and flux conditions. To account for possible initial subcooling or superheating of the medium, the solutions were divided into two stages: a premelt stage and a phasechange stage. The solutions for these stages were then combined into a single equation as
t
T(z,t) To(z,t) ï¿½H(tto)L dS(r) (, t IS(r),r)dr (3.11)
where to is the time when phase change took place, and H(tto) is the Heaviside function given by
0 for t< to
H(tto) = {
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 firstkind boundary condition
E~){2a tr
;qb(T) for the secondkind 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+(n1)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(Sc, t)/at 1o
idt OT(S+c, t)/ax o= OT(S, t)/ax (3.19)
T(Sc,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)/OxOT(S ,t)/Ox 1
dt [T(S+ct)i OT(S+,t)/OxOT(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)/taT(S,t)/Ot (3.23) dt OT(S+,t)/OxT(Sc,t)/x (3LO Equation (3.18) can then be expressed as
[ O T11 _ _pLl[OT(S+ct)/tOT(Sc,t)/Otl
5x Ox]=s  k [aT(S+,t)/OxOT(Sct)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(SCt) _U(S+Ct) OU(S,t)0 as c40 (3.25)
Ox Ox  Ot Ot
Equations (3.23) and (3.24) can then be reduced to
dS _ OV(S+E,t/OtOV(Sf,t)/Ot (3.26) dt OV(S+,,t)/OxOV(S,,t)/Ox (O
[OI' aV pILI dS(t) (3.27) [Ox 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 Uproblem 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 multidimensional 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 OX)
t > to, z E
where xd is the SturmLiouville 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, 5x '
t> to, XE Q
V(X, to) = 0
[a Ov J [ a Ox 5I1=
_pill dS(x)
k dt
(3.40)
The Uproblem, 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 Uproblem.
To solve the Vproblem, the sourceandsink 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 Tproblem is decomposed into a Uproblem and a Vproblem. The Uproblem satisfies the nonhomogeneous initial and boundary conditions (see Eqs.(3.33) through (3.35)), while the Vproblem 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 Vproblem, which incorporates (3.36) and (3.40), a reduction of one equation. Then the Uproblem 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 Uproblem 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 Vproblem 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 multipoint GaussQuadrature 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 NewtonRaphson 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 slowstart 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 Stefanproblem examples in Cartesian, cylindrical, and spherical systems. For the examples in Cartesian system, both semiinfinite 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
Uproblem:
The governing equation for the Uproblem 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 xj1 + ( +=x)t/Ax2.  He2rXej+, th _seXUi 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 secondorder accurate in space and firstorder accurate in time and is often referred to as the simple implicit method. Equation (3.45b), however, is secondorder accurate in both space and time and is commonly known as the CrankNicolson 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 backwardsubstitution, an efficient numerical procedure.
Vproblem:
The solution of the Vproblem 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 semiinfinite 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 (secondkind) and temperature (firstkind) 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 = (2n1)
2Hi' 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) =
(2n1) 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,)
ekn' [(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
Uproblem:
The governing equation for the Uproblem 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 xO (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.
Vproblem:
The solution of the Vproblem 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) =  En  )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
Uproblem:
The governing equation for the Uproblem 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.
Vproblem:
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(Eeak3n(tlT)sin(,3nx')sin( nx) (3.64) n=1
For a temperature condition imposed on the spherical surface, 1_2 #nwn NRs' =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 Kh 1
k Rs
42
The solution for the Vproblem 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 ckto 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 Uproblem 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)dr
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(1cosA) (3.76a) if the simple implicit method (3.45a) is used, and ,n+1  1O(1cosA) n (3.76b)
1 1+0(1cosA)
if the CrankNicolson 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(xS) 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
d1 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 1O(1cosA) 1).
OEcosA) (1 (3.87b) 1+0(1OA
for the CrankNicolson 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(lc osA)]
and
A  (1cosA) (1H (3.88b) 1+O(1cosA)
for Eqs.(3.45a) and (3.45b), respectively.
Since
1+20(1cosA) , 1+0(1cosA) 1
for all 0 and A, the stability condition hinges on the value of 1i11. 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 11H[<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 Uproblem is solved by a finite difference method, while as the Vproblem is solved by an analytical method. The numerical scheme will be stable and consistent If both the U and the Vproblems are stable and consistent. It has been well established that the finitedifference equations, (3.45), for the Uproblem are consistent with their differentialequation counterpart in (3.33). The solution of the Vproblem is exact up to Eq.(3.43). In (3.44), the velocity of the interface is approximated as a finitedifference. 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 MultiDimensional 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 twodimensional
Cartesian coordinate system
l0T =2T 02T L O!y
0t X2 + 6(yS(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 Uproblem
(3.97)
(3.98)
i=1,2, ...
(3.99)
For the Vproblem, the governing equation is
10 V 02V + 2V L Oy a at 2'a2 + 6(yS(z,t))a
t>to, (X,y)EQ
The initial, boundary, and interface conditions for the Vproblem 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 Uproblem 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 ydirections are calculated alternately. The finite difference equations are
Sx il 1l_, J + (1l+20x)6,junl2 il,jnl2= _v m 1 +(120Y)U!i +O T!'~ l
0 U!+ n1l n+l 0Un31 + (12 )+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 Vproblem can be expressed again in terms of a Green's function as
t1 H
r=t0 x0
The interface position curve S(x,tl) is then determined implicitly by the following equation
tl H
TmU(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(xAx, 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 rewritten as follows
1 aV(Y, t) _ a2Vi(y, t) Vi+l(y, t)  2Vi(y, t) + Vjl(y, t)
a t  y2 + Ax2
6(ySi(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 Vproblem, the solution of V can be written as
t1
v (y~t C:,)(
r=at0
ti D
+af drif gi(y',r)Gi(y,tl y',r)dy' (3.109) rtO y'=O
where
gi(y',r) = V , +i(y', 7)  2Vi(y',7) + VjI (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) rto)+..., tO_ r
If higher order terms are neglected, then g can be expressed as
Og(y t o)
M3Y,T7) = M,yto)+ at (rto) (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 rewrite (3.113) as
A 2 6y s.
a t C abt 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+12wj+wj.._) (3.115)
where
t1
W f ( O Sa,(o)l~ iIS
w = jarto)t " i(y,tlS(t),r)dr, K=i,iï¿½1 (3.116)
r~to
In practice, a firstorder approximation is sufficiently accurate to solve the problem at hand. For the firstorder 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 addressedthe phasechange interface may be moving in zdirection (Figure 3.3), and the phasechange interface may be moving in rdirection (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 zdirection
I r
RLiquid
Solid
0
Figure 3.4 System analyzed in cylindrical coordinates:
Interface moves in rdirection
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 Uproblem 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+j2Ui, +Uilj)
+ Or'(Ui, j+1  Ui, j1) + 20r(Ui, j+1  2Ui,j + Ui, j1)
(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 Vproblem 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 zdirection, the governing equation becomes
lg 02V lOI )+ L a
aPt r az rZr r(zS(r,t)) (3.129)
The governing equation is split into two directions, where a finite difference method is used in the rdirection. This reduces the V equation to the following IaVj =02VJ L OZ (a =  a (zSj(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 firstorder approximation of gj leads to a very simple approximation as gj(z,r)  gj(z, to) = 0 (3.132)
The Vproblem can be solved as
V =ZLTl = i Gj(z,tlISj,,r)dr (3.133) r=t0
Similarly, the secondorder 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 rdirection, the governing equation of the equivalent problem is given by
1a8  =2 raVY)+k. L 6(rS(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 zdirection. The V equation becomes
1OVi(r, t) 1 _o_[_o t)(+,Lt)a ] (rS(t)) + g,(r,t) (3.138)
where
gi(r,t) = V,+l(r ,t)  2V,(r, t) + Vi,(r, t) (3.139) Az2
The firstorder accurate solution of the Vproblem is readily obtained as t
Wi(r,ti) =LC 'OSi(,)Gj(r, tltISi, 7)d7 (3.140) .r=t0
and the secondorder 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 (rto)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 Uproblem 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 Uproblem can be derived as
U.+tI  = 20x(Ui+1,, ,j  2Uj, jk+ U ,j, k) " 20y(Ui, j+1,k  2Ui, j,k + U, j1,k) + 20z(Ui, j, k+1  2Ui, j, k + Ui, j, k1) (3.146)
where O, O., and 0. have been defined previously. This equation can be solved by using an iteration procedure such as GaussSeidel and SOR method, or by a modified ADI method.
Again, the Vproblem 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 ydirections, Eq.(3.147) can be rewritten as
Ly,,, 2Vi, +L L t(zS) +g, j(z,t) (3.148)
aat acz2 ct
where
gi, j(z, t) = Vj+l,j2V,i, +Vi,I V, j+1 2VI, j+Vi, j1 (3.149) Ax2 Ay2
With the help of the homogeneous initial condition for the Vproblem, the firstorder 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 Vproblem 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 secondorder accurate approximation of gi,j leads to a solution as
_i ass.,
Vi, jZ' t' = L (Z ztj I Si, j, r)dr
r=t0
+ aL 2(Wi+l,,j_2wi, j + Wi1, 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 multidimensional 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 Vproblem 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 Uproblem 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 phasechange problems solved in the preceding sections by the sourceandsink 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 timedependent phasechange 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 NavierStokes equations and the energy equation. However, a direct numerical simulation of the NavierStokes 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 phasechange 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 steadystate conjugate problem without phasechange, 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 rewritten 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 phasechange temperature, and the Ay is changed to S.
The energy equation (3.154) can then be written as
OT, OUT T1T1 5t Ox + PcjRd 
(3.163)
This equation can be solved by using a firstorder 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 phasechange 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 opensystem 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 phasechange 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 surfacetovolume 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(T8T=k aT = (3.174) h(sT)= Par =p=
In this effort, the convective heat transfer coefficient h for the packed spheres is obtained by curvefitting 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(TTs) = 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 PHASECHANGE 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 phasechange 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 phasechange 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

PAGE 1
SOURCEANDSINK 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 MultiDimensional 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 PHASECHANGE 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 Longtime 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 rdirection 58 3.5 A schematic drawing for the thermal energy storage unit . . 68 3.6 Simplified system for analysis 68 37 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 timevariant 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 timevariant flux condition 102 Interface positions in finite domain 104 Temperature history at the insulated wall ( xH ) 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 meltingfreezing 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 xH /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 coalescingf 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 Postprocessor for a two dimensional TES unit 142 5.32 Postprocessor 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 ydirection 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 SOURCEANDSINK 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 sourceandsink 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 diff 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 Phasechange problems involving melting and freezing have found numerous applications in science and engineering. Phasechange 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 volumetoenergy 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. Phasechange 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 phasechange 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 phasechange problems thus pose a real challenge for analysis. Phasechange problems that have been solved in the literature have often been treated with various degrees of simplification. Earliest phasechange 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 timevariant 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. Phasechange 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 phasechange temperature is not distinct. There is a crystallinelike 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 phasechange 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. Phasechange 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 finitedifference method in one direction and a sourceandsink 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 phasechange interface
PAGE 19
4 moves rapidly. The sourceandsink 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 sourceandsink 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 NeumannStefan formulation. The work involved in the solution can thus be reduced considerably for the sourceandsink method. In addition, in the traditional method of solution, the interface Stefan condition must be given as a separate condition, whereas in the sourceandsink method, this condition has been incorporated into the governing heat diffusion equation. The sourceandsink method is thus highly efficient in the solution of the phasechange problems. In the present work, the sourceandsink 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 longtime solutions. An apparatus has also been built for experimentally testing its accuracy for shorttime 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 phasechange 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 phasechange 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 phasechange 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, diffusiondriven, phasechange solution scheme that combines an analytical solution in one direction and a numerical solution in another direction, only mathematical methods for diffusiondriven 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 phasechange 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 Â•49 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, phasechange 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 phasechange 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 shorttime 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. Phasechange 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 sourceandsink method. Chuang and Szekely (1971) have used this method in the development of a numerical scheme for the solution of Stefan problems. Recently, the sourceandsink 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 sourceandsink 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 sourceandsink 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 sourceandsink 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 sourceandsink 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 timevariant 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 phasechange temperature (i.e., onephase 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 higherorder 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 RungeKutta numerical integration algorithm. 2.1.3 Other approximate solutions Stefan problems can also be solved by using quasi steady and quasistationary approximations. In the quasisteady approximation, the unsteady terms in the governing equations are dropped. This
PAGE 28
13 results in a temperature that is good for longtime 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 shorttime solutions . In the asymptotic method of solution, these quasisteady 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 shorttime solution. On the other hand, Weinbaum and Jiji (1977) have shown that the lowest order term of a longtime perturbation actually leads to a quasisteady limit and is thus good for a longtime 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 firstorder 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. Phasechange problems have also been solved by a Fourier transform (Ockendon, 1974). In this method, the Stefan problems are reduced to an initialvalue problem. In a similar fashion, a Laplace transform can be used to change the problems into boundaryvalue
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 phasechange 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 phasechange 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 timevariant. 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 thirdorderaccurate finite element method with an adaptive grid. In their work, the computational grid is reconstructed 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 fronttracking methods and have been used extensively in the numerical solution of phasechange 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 frontfixing or boundaryfixing. Furzeland (1980) and Lacroix and Voller (1990) have made extensive studies for comparison of the fixedgrid method and the transf ormedgr 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 bodyfitted 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 energyconserving finitedifference scheme in the analysis of freezing outside a cooled tube. Shyy et al . (1992) used bodyfitted 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 RankineHugoniot equation in the shockwave analysis, the weak formulation can be used to solve phasechange 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 phasechange 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 phasechange temperature, the nonlinearity behavior of the enthalpy near the phasechange 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 Phasechange 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 timedependent 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 phasechange 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 timestep limitations. Sadegh et al. (1987) applied BEM to the solution of solidification around a buried pipe in a semiinfinite 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 phasechange 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 sourceands 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 sourceandsink 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 phasechange 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 liquidphase 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 outwarddrawn 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 phasechange 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 sourceandsink 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 TG( 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 semiinfinite domain. In the problems they solved, the medium was initially at a uniform temperature. The boundaries were imposed with constant or timevariant temperature and flux conditions. To account for possible initial subcooling or superheating of the medium, the solutions were divided into two stages: a premelt stage and a phasechange stage. The solutions for these stages were then combined into a single equation as
PAGE 39
24 T(z, t) = T 0 (Â«, t) + H{tt 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(tt 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 tT ^( r ) for the firstkind boundary condition for the secondkind 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 Tt F(r)d T = nÂ— 1 .n1 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 Tt" 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(Se, 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 xS dT(S Â± e, t) dx dx dT(SÂ±e,t) ^ xS 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(Se,t)/dt dT(SÂ€,t)/dx e^o eÂ— >o (3.19)
PAGE 43
28 which can be recast as dT{Se,t) fdT(Se,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(Se,t ) ] dx dx lo ^ (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)/dxdT(Se,t)/dx dt [, dT(S+e,t)/8x ' dTls+e,t)/dxdT(Se,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)/dtdT(Se,t)/dt dT(S+c,t)/dxdT(Se,t)/dx eÂ»o (3.23) Equation (3.18) can then be expressed as Â’dT s dT, P\L  dT(S+Â€,t)/dtdT(SÂ€,t)/dt ' dx dx 1 CO II H dT(S+e,t)/dxdTise,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(Sc,t) _5U(S+e,t) dU(Sc,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)/dtdV(Se,t)/dt dV(S+e,t)/dxdV(Se,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 multidimensional 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 SturmLiouvi 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 Â‘XS (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 Uproblem, the sourceandsink 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)=^ ^ pS 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 Tproblem 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 Fproblem, 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 Fproblem 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 multipoint GaussQuadrature 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 NewtonRaphson 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 slowstart 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 Stefanproblem examples in Cartesian, cylindrical, and spherical systems. For the examples in Cartesian system, both semiinfinite 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 secondorder accurate in space and firstorder accurate in time and is often referred to as the simple implicit method. Equation (3.45b), however, is secondorder accurate in both space and time and is commonly known as the CrankNicolson 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 backwardsubstitution, 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 (xsy (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 (2nl)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\ ^rS(T)G(x,t 1 \S(r),T)dr (3.58) Tt 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. Kproblem: 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) TmU{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( 1cosA) (3.76a) if the simple implicit method (3.45a) is used, and n+1 lg(lcosA) Â„ 1 l+0(lcosA) (3.76b) if the CrankNicolson 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 _ (1n) 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 CrankNicolson 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 (in) e n [1+20(1 cosA)] (3.88a) A l0(lcosA) l+0( 1cosA) ( 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 Kproblems are stable and consistent. It has been well established that the finitedifference equations, (3.45), for the {/problem are consistent with their differentialequation counterpart in (3.33). The solution of the Fproblem is exact up to Eq.(3.43). In (3.44), the velocity of the interface is approximated as a finitedifference. 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 MultiDimensional 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 twodimensional 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(yS(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 ydirections are calculated alternately. The finite difference equations are vr+ ,Y?+ (lwpti* 3 + 1 VWi + 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(xAx, 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 rewritten 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 % KySi(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 rewrite (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 firstorder approximation is sufficiently accurate to solve the problem at hand. For the firstorder 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 phasechange interface may be moving in zdirection (Figure 3.3), and the phasechange interface may be moving in rdirection (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 rdirection
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 zS _ 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 aaAt 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(zS(r,t)) (3.129) The governing equation is split into two directions, where a finite difference method is used in the rdirection. 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 firstorder 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 Fproblem 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 secondorder 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 rdirection, 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 zdirection. 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(Â».02y,(r > 0+V 1 1 (r > 0 ' Â’ A z 2 (3.139) The firstorder accurate solution of the Vproblem is readily obtained as l l ^(^i) = Â£ ^5.(r)G t (r,t 1 5,.,r)dr TÂ—t r ds and the secondorder accurate solution can be derived (3.140) ^>(Mi)=Â£j ^5,(r)G  .(r,l 1 5  .,r)d7 + ~Â£? (Â“i+l 2w i + Â”il) 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 + Uil,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 GaussSeidel and SOR method, or by a modified ADI method. Again, the Fproblem 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 ydirections , Eq. (3.147) can be rewritten as l^.i a dt d 2 V. dz T^ + 53 f S(zS) +g ifj (z,t) (3.148) where 9i, ]{*>*) = Az 2 + Ay 2 (3.149) With the help of the homogeneous initial condition for the Fproblem, the firstorder 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 secondorder 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 multidimensional 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 phasechange problems solved in the preceding sections by the sourceandsink 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 timedependent phasechange 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 NavierStokes equations and the energy equation. However, a direct numerical simulation of the NavierStokes 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 phasechange 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 steadystate conjugate problem without phasechange, 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 rewritten 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 phasechange temperature, and the Ay is changed to 5. The energy equation (3.154) can then be written as 8T f dT f TfT , ~df + U 1h + p fCf R"d = 0 (3.163) This equation can be solved by using a firstorder 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 phasechange 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 opensystem 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 phasechange 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 surfacetovolume 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 VPfCfRs (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\rR p the convective heat transfer coefficient h obtained by curvefitting the experimental as (3.174) for the data of Nu =11 + 0.0 6Pe (3.175) where, Pe is the Peclet number ( RePr ), 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 Nuk 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 .) = Â» (3180) 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 + Qy + $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,nEout = 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 PHASECHANGE 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 phasechange 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 phasechange 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 doublewalled 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. Copperconstantan 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 semitransparent 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 phasechange 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 onedimensional 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, onephase melting in a semiinfinite 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, twophase melting in a semiinfinite 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, twophase melting in a finite domain 4 Stefannumber 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, twophase melting in a finite domain 5 Test of oscillating interface for longtime 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 anNeumann 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, timevariant temperature and flux conditions are imposed on a semiinfinite medium (aluminum), which is initially subcooled (see Table 5.1). A twophase 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 interfacepositionversustime results are shown in Figures 5.4 and 5.5. In Figure 5.4, the medium is imposed with a timevariant temperature condition, while in Figure 5.5, the medium is imposed with a timevariant 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 timevariant 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 = TmTj T u,~Tm x2H' Â— TT 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 semiinfinite 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 semiinfinite 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 semiinfinite 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 preheated prior to melting. Thus, for the final stage of melting, the medium is fully preheated to its melting point; the problem reduces to a onephase 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 (xH) are identical. Clearly, because of the boundary conditions imposed, the interface position rapidly reaches a quasisteady state. The solid curve fluctuates about the dashed curve, which eventually reaches the midplane 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 longtime 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 followup 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 onephase StefanNeumann 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 StefanNeumann 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 einterface 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 meltingfreezing 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 timevariant 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 refreezes 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 twodimensional 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 ^=397100 Â§ at t/=0, other boundaries insulated, T, =T m = 327 wax, onephase 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, twophase 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 spacevariant 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 StefanNeumann 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 xdirection of the V problem, Eq.(3.107), does not introduce any instability in the solution. In twodimensional problems, the time steps are determined
PAGE 134
119 T w =T w( X ' t ) Figure 5.16 System used in twodimensional 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 Longtime solutions Case 2 is intended for testing the convergence of the results to a longtime 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 steadystate 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 longtime solution
PAGE 138
Interface position, S (cm) 123 Figure 5.20 Comparison of numerical and analytical solutions for a twodimensional problem at long time
PAGE 139
124 This steadystate 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 sourceandsink 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 shorttime 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 twodimensional phase change
PAGE 141
126 from the experiment described in Chapter 4, while the bottom figure shows the positions of the phasechange interface. The y axis for this lower figure is intentionally drawn upsidedown 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 phasechange 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 timevariant 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 twodimensional 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 refreezes unevenly in the x direction, giving rise to a two dimensional multiplecoalescingfront problem. The wax melts and resol 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 prestored 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/mK ) 0.6 Heat capacity, c ( kJ/kgK ) 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/mK) 0.58 Heat capacity, c f (kJ/kgK) 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 2gpm flow rate is shown in Figure 5.25. It should be noted that the plots are rescaled 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 phasechange 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 surfacetovolume 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 postprocessors are developed using a QuickBASIC language. Sample outputs of the postprocessors 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 Postprocessor for a two dimensional TES unit
PAGE 158
143 Figure 5.32 Postprocessor 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 sourceandsink method has been developed for the solution of multi dimensional phasechange problems. The method has been shown to be accurate by comparison with exact analysis and experimental measurements. In the sourceandsink 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 multidimensional 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 shorttime and longtime 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 premelting, melting, prefreezing, and freezing stages. Furthermore, the initial condition can be uniform and nonuniform, and the domain of solution can be finite and semiinfinite. In this work, the effects of subcooling on the interface positions in a finite medium have been studied. Problems with oscillating interfaces resulting from timevariant 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 phasechange experiment was conducted for testing the accuracy of the shorttime solutions of the numerical method. An analytical solution was used for the test of accuracy for longtime 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 opensystem 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 postprocessors for graphical visualization of the TES system performance were also developed. The following are recommended for further work: 1. The sourceandsink method with superposition technique will be extended to the solution of three dimensional phasechange 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 sourceandsink 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 sourceandsink 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 (2nl)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(rt 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 ( 2nl)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) = ij ^2 fie 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/Twodimensional phase change problems, Combined FDM and SSM method, Time marching scheme, Combined NewtonRaphson/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(NX1 ) END IF DY = YDIM/FL0AT(NY1 ) 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) = (J1)*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: Xvariant T. B.C. IBC>0: Constant T. B.C. IBC<0: timevariant 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 = TIB2TIB1 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 ) )*(TIMETIB1 )/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( 11 ) )/ (XBC( 1+1 )XBC( 11 ) ) SIG*Y2(Il)+2. = (SIG1 . ) /P = (6.*((TB(I+1)TB(I))/(XBC(I+1)XBC(I))(TB( I )TB( 11 ) )/(XBC( I )XBC( 11 ) ) )/ (XBC(I+1)XBC(I1))SIG*U(I1))/P CONTINUE QN = 0.5 UN = (3. / ( XBC ( NBC ) XBC( NBC1 )) )*(YPN(TB( NBC) TB( NBC1 ) ) / ( XBC (NBC) XBC ( NBC1 ) ) ) Y2(NBC) = (UNqN*U(NBCl))/(QN*Y 2 (NBCl)+l. ) DO 20 I=NBC1 ,1,1 Y2( I ) = Y2( I )*Y2( 1+1 ) + U(I) CONTINUE DO 50 1=1, NX XX = (I1)*DX IF(XX.LE.XBC(1)) THEN TWC = TB( 1 ) + 273.0 ELSE DO 30 J=2,NBC IF(XX.GE.XBC(J1). AND . XX . LE . XBC ( J ) ) THEN H = XBC( J)XBC( Jl ) A = (XBC( J)XX)/H B = (XXXBC( Jl ) ) /H TERM1 = A*TB( Jl )+B*TB( J) TERM2 = ( ( A**3A)*Y2( Jl ) + (B*=t<3B)*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 CSlow 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* (TIME0 . 5*DT) /PRID) TWO(I) = T( I , 1 ) Â’Â’ Â’ DO 100 J=1 ,NY T0( I , J) = T( I , J) 100 CONTINUE C CSolve for Uproblem 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 .E8) GOTO 600 C CCheck 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 Vproblem, 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 Vproblem, 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 . (SiMelting, S2Freezing) 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 Vproblem, 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),(ii)*dx CONTINUE ' 600 700 900 1000 CONTINUE IP = MOD(NCY, ITOUT) DO 700 1=1, NX IF(IP.EQ.O) WRITE(30, 1000) CONTINUE (I1)*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 CTw > Tm : Melting. ITFS(I , 1) = 1 r T . ELSE ! f ( tw ( i )Â™ K lt 0.5.and.TW0(I)TMK.ge.0.5) then tlw < Im : Freezing. ITFS(I ,2) = 1 END IF CInterface 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 CInterface 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.E06 EPST = l.E05 DSDT = ABS(SSO) SSO = S IF(DSDT.LE. 1 .E08) 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(NY1)TI(NY))/DY DO 200 J=2,NY1 DD(J) =3.*(AA( J)*(TI(J+1)TI(J))+CC(J)*(TI(J)TI(J1)))/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 = (FFL)/(RTYL) RTO = RT FL = F IF(F.LT.O. ) THEN YL = RT ELSE YH = RT END IF CDO 20 L=1 ,200 IF(((RTYH)*DFF)*((RTYL)*DFF).GE.O. 1 . OR. ABS(2 . *F) . GT. ABS(DO*DF) ) THEN DO = D D = 0 . 5*(YHYL) 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 = (FFL)/(RTRTO) 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( Il ) D( I) = D(I)RR*A(I1) C1(I)= C1(I)RR*C1(I1) 100 CONTINUE RR = l./D(M2) Cl(M2)=Cl(M2)*RR DO 200 I=MP1 ,M2 J = M2I+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(J1).AND.YY.LE.Y(J)) THEN DL = YYY(Jl) DR = Y( J)YY DL2 = DL*DL DR2 = DR* DR H = Y(J)Y(J1) TERM1 = (3 . *H 2.*DR)*DR2*T(J1)/H TERM2 = (3 . *H 2 . *DL)*DL2*T( J ) /H TERM3 = (HDR)*DR2*D( Jl ) TERM4 = (HDL)*DL2*D( J ) TS1 = (TERM1+TERM2+TERM3TERM4) /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.E04 IF(ABS(YLOC) .LT. 1 .E08) 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 . E08) 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*IT1 )*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 * ( TERM3TERM4 ) 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 .E3) 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=NY1 CSVEEP IN JDIR. 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(I1 , JJ)+BIO*(T(I , JJ)TINF)) ELSE D2TDX2 = C0EX*(T(I+1 , JJ)2.*T(I,JJ)+T(I1,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 CSVEEP IN ID 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 , J1)T0(I , J)) ELSE ENDÂ°IF = C0EY *C TO C 1 , J *i ) 2 . * T Â°( 1 , J)+TO( I , Jl ) ) 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 = ( I1)*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.E06 EPST = l.E05 551 = SI 552 = S2 DSDT1 = SlSOl DSDT2 = S2S02 DIRO = S01S02 IF(ABS(DSDTl).LE.l.E08) THEN DDD1 = 0 . 005*YDIM IF(S01.LE.l.E8) DDD1=0 . 1*YDIM ELSE DDD1 = DSDT1 END IF IF(ABS(DSDT2).LE.l.E08) THEN DDD2 = 0 . 005*YDIM IF(S02 . LE . 1 . E8) 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(NY1)) /DY DO 200 J=2 ,NY1 DD( J) =3.*(AA(J)*(TI(J+1)TI(J))+CC(J)*(TI(J)TI(J1)))/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 = RT1RT2 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 = RT1RT2 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(YLlRTl).LT.l.E8) RT1 = YL1 + l.E05 if (ABS(YL2RT2) . It. 1 .E8) RT2 = YL2 + l.E05 DO 20 L=1 ,50
PAGE 186
171 all = (fhlf 1 1 )/( r tly 1 1 ) al2 = (fhlf 11 )/( r t2yl2) a21 = (fh2f 12)/ (rtlyll) a22 = (fh2f 12)/ (rt2yl2) if (abs(all) . It. 1 .e8) then x2 = f hl/al2 xl = ( f h2+a22*x2) /a21 else hj = a22a21*al 2/all if (abs(h j) . It . 1 . e8) 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(sls2) .LT. 1 .E6) 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.E04 IF(ABS(YLOC) . LT. 1 .E08) 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*IT1 )*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 = TERM11TERM21EPT*(TERM31TERM41 ) SUMgl = SUMgl + SUMI1*SIN(BXL) SUM I 2 = TERM12TERM22EPT*(TERM32TERM42) 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 .E3) 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 Vproblem 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 nonelliptic partial ; III. Heat conduction, with change of state, in J. Mech. Appl. Math., Vol. 15, pp. 5362. 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. 7783. 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. 11491163. annl ! Â” A refinement of the heat balance integral method 13571362Â° 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. 16811686. B ' Â• A ' Â’ 19?4 Â’ Â” The embedding technique solidification problems,Â” Moving Boundary Problems Diffusion, Proc. Conf . Univ. Oxford, March 2527, pp. in melting and in Heat Flow and 150172. Bonnerot, R. , and Jamet, P., 1974, Â”A second order finite method for onedimensional 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 47, 85HT7. 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. 195211. 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 twophase 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. 11811195. problems Â” Int. J. Chou, F. C., Lien, V. Y. , and Lin, S experiment of nonDarcian convection in sphere channelsÂ— 1. Forced convection,Â” Vol. 35, pp. 195205. 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. 12851294. 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. 177200. in V., thermal Morgan , Crank, J. , Oxford. 1984, Â’Â’Free and Moving Boundary Problems,Â” Clarendon Press, Duda, J. L., and Vrentas, J. S. , 1969, diffusioncontrolled moving boundary probl 24, pp. 461470. Â’Â’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. 5056. 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. 411429. 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. 335341. The heatbalance 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. 7179. 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. 1624. Goodrich, L. E., 1978, Â’Â’Efficient numerical dimensional thermal problems with phase change Transfer, Vol. 21, pp. 615621. technique Â” Int. J. for oneHeat Mass Gupta, R. S., and Kumar, method for onedimensional Eng., Vol. 23, pp. 101109. D. , 1980, Â”A modified variable time step Stefan problemÂ”, Comp. Meth. Appl. Mech. Gupta, R. S., and Kumar, D. , onedimensional Stefan problem Heat Mass Transfer Vol. 24, pp. 1981, Â’Â’Variable time step methods for with mixed boundary condition,Â” Int J 251259. Hastaoglu, M. A., 1986, Â”A numerical solution to problem Application to melting and solidification,Â” Transfer, Vol. 29, pp. 495499. 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, Onedimensional 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. 9199. Numerical methods for the boundary element problems,Â” Metallurgical Hsieh , C. K. , 1993, Â’Â’Unification boundary element for the solution Transfer, in press. of sourceandsink 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. 524528. oneand constant Â” J. Heat twophase or timeTransfer , Hsieh, C. K., Choi, C. Y. , 19926, energy storage for solar energy Vol. 114, pp. 203276. 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. 473490.
PAGE 192
177 Hsu, C. F., Sparrow, E. M. , and Patankar, S. solution of moving boundary problems by boundary controlvolumebased finitedifference scheme,Â” Transfer, Vol. 24, pp. 13351345. V., 1981, Â’Â’Numerical immobilization and a Int. J. Heat Mass Kays, V. M. , and Crawford, M. E. , 1980, Convective Transfer, McGrawHill 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. 131. Lacroix, M. , and Voller, V. R. , 1990, Â’Â’Finite di solidification phase change problems: Transformed Numerical Heat Transfer, B, Vol. 17, pp. 2541. 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. 24242428. 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. 97116. steel,Â” Proc. Murray, V. D. , and Landis, F. , 1959, of transient heatconduction problems Part I Method of analysis and sample 81, pp. 106112. 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. 208209. Â” 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 sourceandsink method for the solution of Stefan temperature and flux conditions,Â” Oberkampf, V. L. , 1976, Â’Â’Domain mapping partial differential equations,Â” Int. J pp. 211223. for the numerical solution of Numer. Meth. Engng . , Vol. 10, Ockendon, J. R. , 1974, Â’Â’Moving boundary problems diffusion,Â” Proc. Conf. Univ. Oxford, March 2527. in heat flow and OÂ’Neill, K., 1983, Â’Â’Boundary boundary phase change problems pp. 18251850. 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. 24542456. problems Pham, q. T., 1985, Â”A fast, unconditionally stable finite scheme for heat conduction with phase change,Â” Int J Transfer , Vol. 28, pp. 20792084. difference Heat Mass Pham, Q. T., 1986, Â”The use of lumped capacitance element solution of heat condition with phase change Â” Transfer , Vol. 29, pp. 285291. in the finiteInt. J. Heat Mass Poots, G., 1962, Â”An approximate treatment of heat involving a twodimensional solidification front Â” Transfer, Vol. 5, pp. 339348. 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. 101109. Riaz, M. 1978, Â’Â’Transient analysis of packedbed systems,Â” Solar Energy, Vol. 21, pp. 123128. 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. 294299 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. 223232. Sarler , B. , Alujevic, A., and Kuhn, multidimensional Stefan problems by BEM,Â” Mathematik und Mechanik, Vol. 71, pp. 611614 , . 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 405416 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. 12291245. , M . H. , 1992 , casting,Â” Int. Sikarskie, D. L., and Boley, B. A., 1965, Â’Â’The solution twodimensional melting and solidification problems,Â’ Struct., Vol. 1, pp . 207234. 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. 187202.
PAGE 194
Tamma , K . K , perspectives applications Vol. 30, pp. 179 and Namburu, R. R. , 1990, "Recent advance, trends and new via enthalpybased 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 phasechange 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 preelementary education: English, in the University of Florida. When he was 30year 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
