THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING DIFFERENTIAL ANALYZER SOLUTION OF TEE DIFFUSION EQUATION WITH A FREE BOUNDARY Orlando Joseph Manci, Jr. A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The University of Michigan 1960 June, 1960 IP-439

Doctoral Committee: Professor Robert M. Howe, Chairman Professor Robert C. F. Bartels Assistant Professor Edward 0. Gilbert Associate Professor Elmer G. Gilbert Associate Professor Donald T. Greenwood ii

ACKNOWLEDGEMENTS Many individuals and organizations have contributed materially to the preparation of this dissertation, and I am indebted to all of them. I would like to express my special gratitude to the following: The United States Air Force for making this continuation of my education possible. The members of my doctoral committee Professor R. M. Howe, Chairmang Professor R. C. F. Bartelsj Associate Professor D. T. Greenwood; Assistant Professor E. O. Gilbert) and Associate Professor E. G. Gilbertfor the incentive which they provided as well as their help and encouragement throughout the preparation of this dissertation. The staff and faculty of the Air Force Institute of Technology and especially to Professor H. Co Larsen, Head of the Department of Aeronautical Engineering, for his encouragement, support, and suggestions Associate Professor R. T. Harling and Mr. Floyd Hild, for making analog computer time readily available and for the excellent maintenance and service on the computers; Major V. S. Haneman, formerly of the Department of Aeronautical Engineering and now of the Air Force Ballistics Division; Captain C. A. Anderson, Department of Aeronautical Engineering; Dr. A. J. Shine, Head of the Department of Mechanical Engineering; Dr. H. D. Colson, Captain J. Jones, and Captain C. Richard, Department of Mathematics; and Dro W. L. Lehmann, Head of the Department of Physics; and Dr. G. John, Department of Physics; for their many helpful discussions of the problem and suggestions; Captain B. S. Morgan, Department of Aeronautical Engineering, for his special interest and assistance; and to Captain J. S. Goodman and Associate Professor I. W. Surrat, Department of Humanities, for reading and correcting my rough draft -- any grammatical errors remaining are my own. Mr. M. Tikson, Chief, Dr. M. S. Edwards, Chief Analyst, and Captain B. Rodstein of the Digital Computation Branch, Wright Air Development Division, for the numerical computation of portions of the solution of the warming phase of the problem. Dr. M. G. Scherberg, Aeronautical Research Laboratory, for making digital computer time available. Mr. L. M. Warshawsky, Chief of the Analog Computer Branch, Wright Air Development Division, for making special equipment available for this work. The University of Michigan Industry Program of the College of Engineering and its personnel for rendering this dissertation in an excellent form on very short notice. Finally, to my wife for her encouragement and her patience throughout this period of study and preparation. iiei

TABLE OF CONTENTS Page ACKNOWLEDGEMENT................................................ ii LIST OF TABLES........o..o..................o................ vi LIST OF FIGURES...O..........O...................... vii NOMENCLATURE.......................... o.o...o...oo.........o x CHAPTER I INTRODUCTION................. o~............e.... 1 Backsgroundo.............................. 1... II THE SEMI-INFINITE SLAB..................................6 The Physical Problem................................ 6 The Mathematical Model............................... 8 The Equations in Non-Dimensional Form................ 10 Reduction of the Problem to a Finite Space Interval and Removal of the Moving Boundary.......... 16 Steady State Boundary Velocity for a Constant Heating Rate.............................. 18 The Finite Difference Approximations............... 20 Interior Station Finite Difference Approximations.... 22 Left Boundary Approximation.......................... 24 Right Boundary Finite Difference Approximation....... 24 Determination of the Temperature of the Heated Boundary................... 25 Determination of the Boundary Velocity and Position During Melting..................................... 27 The Finite Difference Equations..................... 28 Analysis of the Equations During the Warming Phase,.. 30 Analysis of the Melting Phase........................ 38 Differential Analyzer Scaling and Circuits........... 41 III COMPUTER SOLUTION FOR THE SEMI-INFINITE SLAB............ 47 Comparison of the Boundary Velocities................ 53 General Solution Curves for Constant Heating Rate.... 62 Solutions for a Time Varying Heating Rate............ 62 IV THE FINITE SLAB......................... 68 The Physical Problem................................ 68 The Mathematical Model............................... 69 iv

TABLE OF CONTENTS (CONT'D) Page Choice of Scaling Constants K1, K2, K3............... 71 Problem Variables.................................... 71 Removal of the Moving Boundary........................ 72 The Nonlinear Space Variable Transformation.......... 73 The Equations in Finite Difference Form...........75 Scaling for the Differential Analyzer................. 77 Nonlinear Space Variable Transformations............ 79 Computer Solution Results...,............... 84 V CONCLUSION................. 90 Description of the Solution......................... 90 Extensions of the Solution............................ 92 APPENDIX A - DIMENSIONAL ANALYSIS........................ 95 APPENDIX B - STEADY STATE SOLUTION FOR THE SEMI-INFINITE SLAB WITH A CONSTANT HEATING RATE.............. 100 APPENDIX C - FINITE DIFFERENCE APPROXIMATIONS AND ERRORS........ 103 APPENDIX D - EIGEN VALUES OF THE MATRIX A....................... 113 APPENDIX E - GENERAL SOLUTION CURVES FOR CONSTANT HEATING RATE AND UNIFORM INITIAL TEMPERATURE...,........... 115 APPENDIX F - SOLUTIONS FOR HEATING RATES OF THE FORM H = h(l-cos2it)....... 124 APPENDIX G - SOLUTION CURVES FOR THE FINITE SLAB................ 131 APPENDIX H - THE EXACT VALUE OF tm, FINITE SLAB................ 148 BIBLIOGRAPHY................... oo oo*.....o... e o 149 v

LIST OF FIGU~T'RES Figure Page Physical Schermatic Semi.infinite Solid......................... 6 2 The (x.9) Domairn......................... 9 3 The (x,t) Domain.................................... 18 4 Values of X'-n IDring Warming for Sevreral Values of r......... 21 5 Finite Difference Statior!s in the Interior of the Slab..... 22 6a-d Temperatures During the Warming Phase vrs. Time............... 34- 37 7 Computer Circuit f+c,r Determining the Temperature at the nth Station........................................................COOO OOO O cc 44 8 Computer Circuits for the Determination of Ur................ 45 9 Computer Circuit for the Determination of t;hle Boundary Velocity and Position....................... 46 10 Steady State Boundary Ve locity v1so Heat iLng Rate.............. 52 la Boundary VeL. oc. tiy a.nd. Psi.tio o ovo ooo o o o e, o o.o.o o 54 llb Tempexraturer- a, the Firn. te Difference T o l tat ions vs Tie.... 5 12a Boundary Vel ociy ande Positio,n v Tie. j.. 56. 12b Temrperature at heZ FinLte Difference Stations vs. Time... 57 13a Boundary Velocity and Po sition vs. Time... c.............ODO 58 13b Temperature at the tr'';z te Difterence avs. Time.D....... 59 14a Boundary Velocity and Posz.tion vs. i-m O..e........D 60 14b Temperature at the Finitee Difference Stations vs. Time........ 6-L 15 Physical Bchemat c.: Finite Slab...... o..O....O............O... 68 16 The (xt) Domain, Finite Slab........................... 69 vii

LIST OF FIGURES (CONT'D) Figure Page 17 The (x,t) Domain, Finite Slab............................. 73 18 The Nonlinear Transformation.............................. 74 l9a-d Computer Schematic for the nth Station.................... 80-82 El Temperature vs. Time.................................... 116 E2 Temperature vs. Time................................. 117 E3 Temperature vs. Time.............................. 118 E4 Temperature vs. Time..... 119 E5 Boundary Velocity and Position vs. Time................... 120 E6 Boundary Velocity and Position vs. Time................... 121 E7 Boundary Position and Velocity vs. Time................. 122 E8 Boundary Velocity and Position vs. Time.................. 123 Fl The Heating Rate Pulse............................. 125 F2 Temperature vs. Time.................................... 126 F3 Temperature vs. Time.................................. 127 F4 Boundary Velocity and Position vs. Time................... 128 F5 Temperature vs. Time...................................... 129 F6 Boundary Velocity and Position vs. Time................... 130 Gla Temperature vs. Time.............................. 132 Glb Boundary Velocity and Position vs. Time................... 133 G2a Temperature vs. Time...................................... 134 G2b Boundary Velocity and Position vs. Time................. 135 viii

LIST OF FIGURES (CONT'D) Figure Page G3a Temperature vso Time.e......... o..............o... 1.36 G3b Boundary Velocity and Position vs. Time1... e......... 1.37 G4a Temperature vs. Time....................................o.. 138 G4b Boundary Velocity and Position vs. Time.O.0..............O. 139 G5a Temperature vs. Time....o.................................... 140 G5b Boundary Velocity and Position vso Ti.me~.......,.......... 1.4_ G6a Temperature vs. Timeo.................o.o...o.,....... l.... 142 G6b Boundary Velocity and Position vs. Time..........,-......... 143 G7a Temperature vs. Time.................................... 1.44 G7b Boundary Velocity and Position ar.d Heat Pulse vs. Time...... 1.45 G8a Temperature vs. Time........................................ 1i46 G8b Boundary Velocity and Position and Heat Pulse vso Time.,..... 1.47 ix

NOMENCILATURE (-) indicates a quantity in the real dimensional space, (x,t). [Equations (1)] C() indicates a quantity in the non- dimensional space, (xt). [Equations (2)] (A) indicates a quantity in the (,t") space for the finite slab. [Equations (38)] indicates the inverse of the quantity in the parentheses. indicates approximately equal quantitieso b steady state boundary speed. [Equation (11)1 B a matrix of eigenvectors of the matrix A. [Equation (27)] c specific heat of the conducting medium. [Equations (1)] c maximum value of the specific heat of the conducting medium. 0 [Equations (2)] C constant in the nonlinear power transformation of the space variable. [Equations (45)] D length dimension (Appendix A)o dim( ) dimension of the quantity in the parentheses (Appendix A). e the base of the natural logarithms. [Equations (12)] [x]t kit e a diagonal matrix with elements e i [Equation (27)] f(x) functional notation for the nonlinear space variable transformation. [Equations (40)] f(t) a column vector describing the boundary conditions in Equation (26). Fn heat flux passing the nth finite difference station (Appendix C), the temperature at x=; in the (x,t) and (x,t) spaces. [Equations (1)] go the temperature at x=O in the (x,t) space. [Equations (9)] g(x) initial temperature distribution. [Equations (1)]

NOMENCLATURE (CONT D) g a column vector of initial temperatures. [Equation (27)] Gn initial temperature distribution in the computer variables. [Equations (33)] h maximum amplitude of the heat pulse H=h(l-coswt)(Page 62 ) H(t) heat flux entering the conducting medium (Figure 1). A>) the heating rate in the compute variables. [Equations (33)] J Jacobian of a transformation. k exponent in the nonlinear space variable transformations. [Equations (45) and (46)];k thermal conductivity of the conducting medium. [Equations (1)] ko maximum value of the thermal conductivity of the conducting medium. [Equations (2)]. K1 scaling constant for the space variable. [Equations (2)] K2 scaling constant for the time variable. [Equations (2)] K3 scaling constant for the heating rate. [Equations (2)] e(t) position of the moving boundary as a function of time (Figure 1). L latent heat of fusion of the conducting medium. [Equations (1)] M mass dimension (Appendix A). n an index designating the nth finite difference station (Appendix C). q a point heat source in the conducting medium. [Equations (1)] Q The heat added to the finite slab (Page 88 ). r the number of finite difference increments. [Equations (13)] s the initial thickness of the finite slab. [Equations (36)] t the independent time variable. t1 the time at which melting starts. [Equations (1)] xi

NOMENCIAITRE (CONTTD) t2 the time at which melting stops. [Equations (1)] t3 the time after which the boundary velocity may be considered to be linear in time (Page 19). tm the time at which the finite slab is completely melted. [Equations (36)] T time dimension (Appendix A). u the dependent temperature variable (Figure 1). uc the melting temperature of the conducting medium. [Equations (1)] u(t) a column vector of the temperatures at the finite difference stations. [Equation (26)] Un the computer temperature variable. [Equations (33)] x the independent space variable (Figure 1). finite difference increment in x. [Equations (13)] temperature dimension (Appendix A). O(x,t) an arbitrary function (Appendix C). Xi the eigenvalues of the matrix A. [Equation (27)] [X] the non-singular diagonal matrix of the eigenvalues of the matrix A. [Equation (27)] X(T) the boundary position in the computer variable. [Equations (33)] it the constant 3.14159... in Equations (29), (30), and (47). ET a dimensionless ratio in Appendix A. p the density of the conducting medium. [Equations (1)] Po the maximum value of the density of the conducting medium. [Equations (2)] T the independent computer time variable. [Equations (33)] xii

NOMENCLATURE (CONT' D) Oc the non-dimensional specific heat function. [Equations (2)] the non-dimensional conductivity function. [Equations (2)] the non-dimensional density function. [Equations (2)] xiii

CHAPTER I INTRODUCTION Background The problem of the conduction of heat in a solid which is freezing at its boundary was first posed by Fourier in connection with his study of the solidification of the eartho(4) Crank(4) and Evans, Isaakson, and Macdonald(l2) list several other problems which are governed by essentially the same equations. Recently9 these problems have received renewed interest in such diverse areas as the tarnishing of metals(4) and the ablation of material from bodies entering the earth's atmosphere. The first sol.ution of the freezing problem was obtained by Neumann about 1840.(2,6) This solution was for a semi-infinite body of liquid which initially existed at the'temperature of fusion. In 1891, Stefan(38) solved the problem again. for the semi-infinite body, for two particular initial conditions. This freezing or melting problem has sirce (6) been called the Stefan probl.em;() other similar problems have been called Stefanlike problems; (12) and the published literature is profuseo The article by Brillo-uin(l) contains an excellent bibliography of the li'terature up'to 1929o The existence and the uniqueness of the solution to this problem (11) Douglas (?) have been shown for very general circumstances by Evans, Douglas(7) Miranker,(30) Miranker and Frisch.,(31) and Kyner. (26) The general. solution has not yet been obtained, and it is not in the offing because the problem is nonlinear, as LandauL showso(27) -1

-2Several schemes for solving the problems at least approximately, have been devised. These analytical solutions are for the semiinfinite body and apply only for particular initial and boundary conditions. (All of these solutions provide for the two phases to remain in contact after the change of phase.) Datseff(6 has developed a step-by-step solution which is very laborious and not completely general. Kolodner(24) has reduced the problem by eliminating the diffusion equation and reducing the problem to that of solving a nonlinear, ordinary, integro-differential equation. Friedman(l5) has refined a method developed by Rubinstein(36) which reduces the problem to another nonlinear, integral equation that is solved by a method of successive approximations. Douglas and Gallie(8) have developed an iterative process based on a difference analog for a simple case which does not allow for heat flow in the material before the change of phase. Ehriich(lO) has developed an iterative scheme again based on a difference analog which does provide for heat conduction before the change of phase, and Crank(5) has used an iterative process which allows this heat conduction. Crank, in this same article, uses Lagrange's interpolation formula in a second method to allow the variation of finite difference increment size when the boundary passes through a given finite difference cell. Longwell(28) has devised a graphical solution using the Schmidt construction. The ablative problem, in which the melted material is removed instantaneously, was first studied, according to Landau, (27) by Soodak.(42) Landau(27) has developed a solution for the melting of a semi-infinite solid with perfect ablation, in which none of the input energy is used to heat, the

liquid. In this solution he considers the case of a prescribed temperature at infinity, a constant heating rate at the finite boundary, and a uniform initial. temperature distribution. He solves the warming phase, before melting starts, by means of the classical analytical solution as given in (2) Carslaw and Jaeger. The transient portion of the melting phase is solved numerically. This transient solution is then matched to the steady state solution, which he obtains analytically, by means of an empirical. equation. (16) (40) Goodman and Sutton have studied the semi-infinite body with imperfect ablation in which part of the input energy is used to heat the liquid. Goodman and Shea(l7) have applied the heat integral method to the problem for finite slabs, and Citron(3) has developed a method of successive approximations for the solution of the finite slab melting with ablation. In all of the above methods some form of numerical computation is required. This ranges from the use of large-scale digital computers as required by the iterative schemes through the numerical solution of a Volterra equation of the second type required by Kolodner's solution to the use of a desk calculator in Citron's method. The analog computer has also been used in recent investigations of this problem. Martini(29) proposed a scheme in which the stations in the finite difference approximation would be switched out of the circuit as the moving boundary passed them. Mirray(32) also uses a similar technique in his second method while in his first method he uses a moving network of stations. Aldrich and Paynter(41) designed a special analog computer for

the study of frost penetration in the soil. Otis(34) developed a passive network for the solution of the finite thickness slab with perfect ablation. In order to eliminate some of the nonlinear equipment in this computer; Otis employed a nonlinear time transformation which requires an infinite time to melt the slab completely. More recently, Sunderland(39) has used passive networks to solve the melting problem in both finite and semiinfinite slabs for the ablative and non-ablative cases. He attached the coordinate system to the moving boundary in the same manner as Landau (27) In the past year, Murray(32_ 33) has used the electronic differential analyzer to solve the melting problem for the finite slab without ablation. His method employs a network of stations which move within the slab as melting progresses. The spaces between the stations change in size as melting progresses but remain equal to each other in each phase. Murray compares his results with solutions obtained by means of approximations suggested by Eyres, et al.(13) This article by Eyres et al, is the earliest reference that this author has found to suggest the use of a space network in which the increment size varies with time after melting starts. The solution of the melting problem presented in this dissertation requires the use of the differential analyzer to solve a set of ordinary) time-varying, nonlinear differential equations. The digital differential analyzer could also be used or the equations could be solved numerically. This method provides for the solution of the melting of the semi-infinite and finite slab problems with perfect ablation. The original contributions

-5contained are believed to be: 1. a new set of non-dimensional parameters which separates the important variables of the problem; 2. the reduction of the semi-infinite interval to a finite interval in order to use the finite difference techniques and the differential analyzer; 3. the development of a method to make certain that the finite difference approximation solution approaches the steady state solution for the continuous case for a constant heating rate in the semi-infinite slab; 4. the use of non-uniform finite difference increments to improve the accuracy of the solution in the finite slab solution. The solution presented here also provides for arbitrary heating rates (which may be functions of the solution as well as of time) and arbitrary initial temperature distributions. Although only the ablative case is discussed here, the solution may be extended to the non-ablative case. This would require only slight- modifications but would approximately double the amount of computing equipment needed. The semi-infinite medium is considered in detail first; this discussion is followed by a consideration of the finite thickness medium. Diffusion in one dimension is assumed in both cases.

CHAPTER II THE SEMI-INFINITE SLAB The Physical Problem The physical problem is here described in terms of heat conduction in one dimension with melting at one boundary. The conducting medium consists of a slab which extends from negative infinity to positive infinity in two dimensions. In its third dimension, x, it extends from a boundary which is a positive finite distance, ~(t), from the origin to positive infinity. Heat enters the slab at x = (t ). This configuration is shown schematically in Figure 1. Conducting Medium x Figure 1. Physical Schematic, Semi-infinite Solid. Let x = the real space variable. t = the real time variable. u(f,t) = the temperature within the conducting medium. H t) = the heat flux entering the left boundary. e(t) = the position of the left boundary. -6

-7Also let 7(0) = o The problem starts when heat is first added with the finite boundary at the origin. If heat is added, i.e. H(t) > 0, and if the medium is at a uniform temperature initially,, then the temperature of the slab increases and u(,t) becomes a monotonic; non-increasing function. of x at any instant of time. The temperature range within the conducting medium is from the temperature on the boundary to that at infinity, which rei.ains at the initial temperature because the adding of heat at the left boundary cannot affect the temperature at a point infinitely removed. This condition is nearly met by bodies of finite thickness in which the heating rate is high and the conductivity is lowo Thus physical systems can approximate ~-ery closely the conditions set down here and the semi.-infinite problem does have practical importance. If the heat is added in sufficient quantity and at a sufficiently rapid rate, the heated boundary temperature will rise to the melting point and the boundary will move to the right. It is assumed that the liquid produced by the melting is removed immnedi.ately and none of the heat flux is used to raise the temperature of the fluid. If the heat flux is removed, melting will stop and the temperature throughout the slab will, approach the temperature at infinity. If the heating rate is decreased sufficiently, melting will stop and the problem becomes the usual semi-infinite heating problem with a fixed boundary.

-8In order to prepare this problem for solution on the electronic differential analyzer the following steps are taken: 1. the mathematical model is established; 2. the problem is made non-dimensional; 3. the semi-infinite interval in space is reduced to a finite interval; 4. the moving boundary is removed; 5. the space derivatives are replaced by finite difference approximations; and 6. the problem is scaled for the computer. The Mathematical Model The problem is described by the diffusion Equation (la) which applies throughout the conducting medium; the boundary condition Equations (lb), (lc), (Id), (le), and (if) which apply at the left boundary during the appropriate time intervals; the boundary condition at x = oo, Equation (lg); and the initial conditions, Equations (lh) and (li). - 0A O t-~ t (a) GH Oc d- 0 t)t, (a)b X=0 L ocLtht- (c) X > (e) {I L(X) C, X,& t=O ( )

-9In these equations: = (,p(u)=the density of the conducting medium. c - c(xu) - the specific heat of the solid. k k(x,u) = the conductivity of the solid. L = the latent heat of fusion of the solid. uc - the temperature of fusion of the conducting medium. t1= the time at which melting starts. t2 the time at which. melting stops. g& = the temperature at x - oo q(x,t) = a point source of heat within the conducting medium, It is assumed in this thesis that there are no heat sources within the conducting medium. Therefore, c(xt) - 0. The time at which melting stops, t2, can be infinite for a heating rate which remains large enough to maintain melting once it has begun to melt. Figure 2 shows the (x,t) domain of the problem for such a heating rate. t1 Figure Figure 2. The (xt) Domain.

- 10Equations (1) will now be made non-dimensional to reduce the number of problem parameters and to provide for a more orderly investigation of the problem. The Equations in Non-dimensional Form To make Equations (1) non-dimensional the following dimensionless variables are defined. Lp(QC1 = p( )= i ( A )x~ ~(a) c Co () CO) (b) ib, RD d (X1M) A5 kt ) (C) -A K K OOD C (d) (2),,rK -i-"- I (d) CO () (e) u(t) = L ( UC-U ) (f) In Equations (2a), (2b), and (2c) the quantities po Con0 and ko are dimensional constants equal to the maximum value of p(x,u), c(xju), and k(x,u) respectively. The functions ~py gcn and ~k are dimensionless functions of temperature and position which are characteristic of the material of which the conducting medium is composed.

-11In Equations (2d)? (2e), and (2f), K1, K2, and K3 are conrenient, constant, arbitrary scaling factors. Equations (2d-2f) are derived in Appendix A. The numerical values of the constants in these equations are included for several common materials. When Equations (2) are substituted into Equations (1), the nondimensional Equations (3) are obtained. - 2 [aJ 2~ (a) HtK dt x- ()t (b) K u < =0 o.< t t (d) ~~~dJQE~~~~~~~OJ, (e) x, -6,4-/;-f (f)'co~;z- tC o (g) M-_. - 0 % - o (h) )~ {~O (h) The V 1/o t o (i) The (x\9t) domain for Equations (3) is the same shape as the (xt) domain of Figure 2 on page 9

- 12Two special choices of the scaling factors K1, K2, K3 are of interest: 1. When H( t) is not constant, it is desirable to eliminate these arbitrary constants from Equations (3)o 2. When H is a constant these factors can be chosen to eliminate H. For the first case define the scaling constants as in Equations (4). K/<2 14K (a) (4) K - (b) When Equations (4) are substituted into Equations (3) the results are Equations (5). % {f0 a4] 3s b o (a) H t)ONd i D X=%(t; 6>0 (b).-'-jT. - 0 o —-O t.<t, (c) dW0 dio pt=0;~4(a<t e (e) t = 0e x=2(t);tl - 2 (f) O. cx - -Xo (g) ~~~~~2i~~~~~~o~~ GT ~O~ U = ~(h)

- 13In order to eliminate H when the heating rate is constant, the scaling factors are defined as in Equations (6). K,- = i (a) a-= H' (b) (6) K3= /-qz (c) This transformation is well defined because g, < 0. When Equations (6) are substit-u ted into Equations (3) Equations (7) result. T < 0 % "6.J 04 X t < ~ (a) d _v d >, Oz (7) (c) LK = O Kc -0 O~Utt t2 (f)'- - I C) t=:o1 (g) ", uX -uO (h) Li =qi x) so, t=0 (i)

-14Equations (5) are the fundamental equations of the problem because, as can be seen by comparing Equations (5b) and (7b), the transformation defined in Equations (6) effectively requires a dimensionless heating rate of (1-g). This dimensionless heating rate may be related to any real heating rate and temperature at x = oo by the proper choice of the constant H in Equation (2g) and Equations (6)~ It will be shown below that this heating rate, H = l-g, also will assure that the steady state melting rate obtained from the finite difference approximation approaches the steady state melting rate obtained from the solution of the partial differential equations when the density, specific heat, and conductivity are assumed to be constant. In Equations (5) the problem is described in terms of the following variables: Lo Independent variables, x and t; 2o Dependent variables, U~Xx t ), Qt) 3o Problem parameters: a. Heating rate H(t). b. Initial temperature distribution g(x). 4. Characteristics of the conducting mediuma. Density function, p (xu), b. Specific heat function, Oc(Xu), c. Conductivity function~, Ok(X9u). This is the minimum number of variables to which the problem can be reduced by means of dimensional analysis according to Buckingham's T theorem(2l19) However, when the heating rate is constant one of the parameters, the dimensionless heating rate, can be eliminated as has been done in Equations (7).

The problem can be further simplified by assuming that 1. the conducting medium is homogeneous, and 2. the density, specific heat, and conductivity are independent of the temperature. These two assumptions make the functions Op, c~, and Ok equal to unity. They also make the problem independent of the material of which the conducting medium is composed. When these assumptions are made the problem is described completely in terms of four variables and two parameters. These assumptions are made at this point and apply to the rest of this dissertation. However, the development of Equations (5) to the finite difference form is carried out in Appendix C for the non-homogeneous, temperature dependent, conducting medium. When the heating rate is constant and the initial temperature distribution is uniform in x the problem parameters are reduced to one. It is thus only necessary to obtain solutions for the range of initial. temperatures which is limited (in the dimensional variable) in minimum value by absolute zero and in maximuim value to the melting temperature of the material being considered. The maximum temperature in the non-dimensional form is zero while the minimum value will vary with different materials. For the materials listed in Appendix A the lowest value is -5.01 for stainless steel. Thus the range of initial temperatures over which solutions need be obtained for a complete solution of the problem is limited when the initial temperature distribution is uniform. Of course, if the initial temperature distribution is arbitrary, an infinite number of possibilities exist.

-16N'- N, When Equations (5) are considered individually with p = c=k:=l the problem seems to be linear. However, the problem is nonlinear as shown by Landau(27) because the position of one space boundary, 4(t), depends upon the temperature distribution which in turn depends upon the position of the boundary, as shown below in the transformed Equations (10). The space interval for Equation (5a) is semi-infinite in extent and the left boundary moves with time after melting starts. It is necessary to reduce the interval to a finite one and to remove the moving boundary before the finite difference approximations can be applied and the differential analyzer used. The next transformation is designed to accomplish these results. Reduction of the Problem to a Finite Space Interval and Removal of the Moving Boundary The transformation defined in Equations (8) will reduce the space interval and remove the moving boundary, X= (a) (8) t (b) Thus x is in the interval 1 > x > 0. Let Lt) %K(9)(a) Hu~(x~~~t)=~~~ H ~t) (b) ac0, (x ( ) =:a;;t (e) q Q = g oD (e)

- 17When Equations (5) with 0p -,c = Ok = are transformed by Equations (8) and the relationships of Equations (9) are used, the results are Equations (10). dA u='o x=; t,ItIt (bf) u=q, x~o, to (g, I J -Co;t-O (h) g(x, O<cXvIxI t=o (i) The Jacobian of the transformation (8) vanishes at x - 0. How7ever, this singularity can be tolerated because only u(x,t) need be defined at this point. The derivatives of u(x,t) need be defi.ned only arbitrarily close to x = oo.

-18Figure 3 shows the (x,t) domain of the problem. t Conducting Medium 1 x Figure 3. The (x,t) Domain. Equations (10) describe the problem with boundaries at x = 0 and x = 1. The point x = o has been transformed into x = 0 and x = 4(t) into x = 1. With the boundaries fixed, the effect of the boundary motion dQ d* is now reflected in the term in Equation (10a) which contains dti dt is a function of the heat flux at x = 1 as shown in Equation (lOb). The (x,t) coordinate system is attached to the melting boundary and the heat transfer represented by the term, - x dt -, can be interpreted as the heat transferred by convection due to the relative motion between the conducting medium and the coordinate system. Equations (10) are now in a form suitable for the application of the finite difference approximations to eliminate the space variable derivatives. Before doing this, however, the steady state melting rate will be discussed. Steady State Boundary Velocity for a Constant Heating Rate Equations (10) can be solved, as is done in Appendix B, to obtain the steady state boundary velocity and the temperature distribution for the special case of a constant heating rate. The result of this solution is given in Equations (11).

~d-6__~'9- - X/ )- (b) For convenience here the steady state boundary velocity is de-, de noted by b = d dt Equations (l1) may be transformed to the (x~t) domain to obtain Equations (12). ((a) I =, go (12) (3Ir7 ),LI ~ Zt )-7 J c (b) where g go as defined in Equation (Qe). The temperature distribution is obtained explicitly as a function of x in Equation (11b). However this does not yield the temperature distribution as a function of x, the untransformed distance variable. This is apparent if the boundary position is written as where t > t3 and t3 is the time after which the boundary position may be considered linear in time. Until (t 3) is computed by some other technique, i(t) cannot be obtained for t > t3.

-20The steady state boundary velocity is unity when the heating rate is chosen to be (1-go), as can be inferred from Equations (7b) and (12a). It will be shown below that the steady state solution for the finite difference approximation is also unity for this heating rate. The Finite Difference Approximations The electronic differential analyzer can. integrate with respect to one variable only, (19) but Equations (10) contain derivatives with respect to both t and x. It is, therefore, necessary to eliminate either the derivatives with respect to t or to xo The elimination of the latter is more convenient because the resulting set of ordinary differential. equations are of first order in t. It is required because it is impossible, in practice, to chosen initial conditions, u(Otn), at each station in time such that the n final. conditions, u(l,tn), are satisfied after integrating from x = 0 to x = lo The derivatives with respect to x are eliminated by using finite difference approximations. The approximation used will depend upon the location of the station at which the derivative is to be approximated. Three different types of stations can be identified, (1) the interior stations, (2) the left boundary station, and (3) the right boundary st ationo Equations (10) are written in such a manner that first order derivatives only need be approximated. To effect the finite difference approximations let the slab be divided into equal space increments and let r = the number of finite difference increments or cells in the interval 0 < x < i

-21Then a x r (a), nZ ax $ 0, 1,2,.., (b) (13) u(t)LA (xrt) n - 0 1,2,)...,r (c) The cell size is constant and is uniform in x. However, in the (x<'jt domain )r X,, =r ti- In) r- (a) while in the (xt) domain (14) (x,,= t~ ) - KIcI nL r (b) Thus the station distribution is logarithmic in both the x and x coordinates. The values of xn and xn are functions of time when melting is in progress. Figure 4 on page 21 shows the station position with respect to the moving boundary for the x coordinate. 1 0o.692 r-2 0.405 1.10 r=3 o.J88 0. 92 -.38. j23.1510.1918 1161 r=5 r=6.62 e t05.;92 1.1.. 85 x Figure 4. Values of x, During Warming for Several Values of ro

-22Station number 0 is defined to be the left boundary in x. This corresponds tothe point at x=+o in the (x) domain as shown in Figure 4 for x. Stations 1. through (r-l) are the interior stations, and station r is at the right boundary in x (left boundary in x as shown in Figure 4). The interior station finite difference approximations will be obtained first. Interior Station Finite Differencee Approximations The central difference approxi.mation is used to obtain t+h.e finite difference a.pproximati.on of U at the nth station~ This derivative is approximated by the rqlationship giver. in. Equation (15)o I.....|.Un...~ un-l - 2 ( Un )(15) Figure 5 shows schematically the quantities of Equation (1i5)o Un+1 un Xn- 1 xn xnr+l n=l., 2 r Figure 5. Finite Difference Stations in the Interior of the Slab

-23It is necessary to use the double interval in Equation (13) because both the temperature and the first derivative with respect to x must be evaluated at each stationo The error introduced by this approximation is of the order of (Ax)2 as shown in Appendix C. It is also necessary to approximate a Ex -] and this is done iX Xn again by using the central difference formula. However, a single interval, Ax, can be used for this approximation which is made in accordance with Equation (16). K~ LXJ Do U x(16) However, since LA t n E an = (n l - U)I and Equation (1.6) reduces to Equation (17). nu,, -r2 nr +(n L -7

Left Boundary Approximation At the left boundary, xo, the temperature is fixed. Thus the left boundary equation is Equation (18). lu~~io~t)~=90~ ~(18) Right Boundary Finite Difference Approximation At the heated boundary, Xr, both the gradient and the temperature must be evaluated. The gradient is approximated by a backward difference. As shown in Appendix C, the backward difference approximation introduces an error of the order of Ax when first order differences are used. When second order differences are used with a single interval, the order of magnitude of the error in the first derivative becomes (Ax)2. Both these approximations will be investigated to determine the effect on the time domain solutions. The first backward difference approximation of the first derivative of the temperature at the heated boundary is given by Equation (19). __) Or- Atop (19) When second order differences are used the approximation is that of Equation (20). ( i)< X1 )9 IU/ (2 2M- L 2) (20)

-25The approximations of Equations (19) and (20) are derived in Appendix C. Determination of the Temperature at the Heated Boundary The temperature at the heated boundary must be determined as a function of time during the period of warming before melting starts. This may be done by means of Equation (10a) or (10b). When Equation (10a) is used, the heat capacity lumped at this boundary station is considered in the calculation of ur, and r differential equations must be solved during the warming period. One difficulty with this method, which will subsequently be termed the implicit method, is that in general d (t l) will have a non-zero value at the instant melting starts. This introduces an error into the boundary position because d-y (tl) (3) for the continuous medium~ must be 0 at tI as shown by Citron for the continuous medium The backward difference approximations used at the heated boundary for the implicit determination of ur are. (a) first order differences ()_ sean ore ifrec(21) (b) second order differences l 73 r - r( t(>r - (22) 2

The derivations of these approximations and the assumptions made in these derivations are contained in Appendix C. Equations (21) and (22) were obtained by using the smallest increments in x possible at each step in the derivations. These formulae are not in as common use as Equations (19) and (20) and other formulae could be developed, but they would require using larger increments in x in some point of their derivations. In an effort to eliminate the error in d (tl), the explicit dt method for determining ur is introduced. In this case Equation (10b), which becomes H =a x )ot0tt during the warming interval because d- = 0, is used. This equation allows dt the determination of ur in-.terms of the heating rate and the internal temperatures. This approximation forces di (tl) to be zero. However, it dt neglects the heat capacity lumped at the boundary station and, if H(O) A 0, ur(O) is in error. During the warming period (r-l) differential equations must be solved simultaneously. When ur is determined explicitly during warming, the boundary temperature during melting can be calculated and this value can be used to determine when melting stops if the heating rate falls below that required to sustain melting.

-27It is9 of course, possible to cal.culate ur implicitly and explicitly simultaneously~ Then the implicit calculation can be used to determine the boundary temperature, ur, when there is no melting; the explicit calculation can be used to determine when melting stops. The approximations for the explicit determination of ur are: (a) first order differences L 9, t+ H(t) (23a) (b) second order differences 4 2 r 3 (r- - 3 - +3r H(t) (23b) These approximations follow from EZTations (1.0b), (19), and (20) since d- 0 i.n the interval 0 < t < ti, Determination of the Boundary Velocity and Position During.Melt ing During the melting phase the temperature, uris identically zero and Equation (lOb) is used to determing the boundary velocity and position. For this case, also, both first and second order backward difference approximations for a- are to be inv-estigated. The approximations for dQ are given in Equations (24a) and (24b). dt First order differences~ df \ (24a) Second order differences: d a r 2 -t 2 r l d-t. c~ ~-l~t(24b)

-28These approximations follow directly from Equations (19) and (20) with ur = 0 and Equation (10b). The Finite Difference Equations The space derivatives of Equations (10) are eliminated by means of the finite difference approximations by introducing Equations (15) through (24) into Equations (10). The results are the following set of ordinary differential equations with boundary and initial conditions given in Equations (25). Uot) -=g, t>/O0 (a) du, 2 (2 ) -, L An-, +2n,(r-O (b) 2 Un,,,+ 2(2r )I L4+, LA 2rH~ r(2r-)r r-I)utr (cl)'-:-= 2 r H- r(2 r-I) L,_,- r(r-i) r'Lr= -j — (r)r- +2r(2r-l) -l-j(7r-3) 14tt (c2) Ur = r H (c3) r = 3rH+ 4,- 3 1 2 H (c4) dQ o ott, (d) r(r<~ o 6 t7Stl (e) o. _ r< t,. ta I,.) d -~ H- r H -(r-,) +-2 r U-, (f2) IA, =._ n i b r ~ ()) 1 r

92 9 Equations (25c1) through (25c4) will be used in turn with Equations (25fl) and (25f2) in investigating the effects of first and second order difference approximations. Two combination.s, however, are not worthwhile. These are (25c3) with (25f2) and (25cL4) with (25fl). Equations (25c3) and (25c4) result from the explicit determination of ur. One object of this approximation is to force Ld (tl) to be zero, and this can be done only if the is approximated by means of the same order differences during the warming and during the melting periods. Thus the combinations in which first order differences are used during warming and second order differences during melting should be eliminated. With these two cases eliminated, there are six cases of interest in studying the effects of the various approximat ions o Equations (25) describe the problem in terms of r simultaneous; first order, ordinary differential equations. During warming dt is zero and the equations reduce to a set of linear, ordinary, differential equations which can be solved by conventional methods~ During melting the equations are nonlinear because dt depends upon the temperatures at the interior stations and they in turn depend implicitly upon the boundary vel.ocity. Equations (25) are -the equations to be solved upon the differential. analyzer during both warming and melting. It is reiterated here that the stations at which the temperatures are determined move within the solid during melting and remain a fixed distance from the physical boundary in the (xt) and (xt) domains but are fixed in the (x t) domain.

-30Equations (25) will be analyzed during the warming phase in the next section. Analysis of the Equations During the Warming Phase In order to determine the errors introduced by the finite difference approximations and to isolate these errors from those introduced by the computer. Equations (25) are solved analytically for the warming portion of the problem and this solution is compared with the solution obtained from the partial differential Equations (5). The solution of Equations (25) is facilitated by writing these equations in matrix form and proceeding to the solution in this form. Thus these equations may be written as m =A a (t) t4(t) o<t t, (26) where u(t) is a column vector with r elements for the implicit determination of ur and (r-l) elements for the explicit method; A is the matrix of Equations (25b) and (25c) for the interval. 0 < t < t1, It is an r x r matrix for the implicit method and (r-l)x(r-l.) for the explicit method. f(t) is a column vector with the same number of elements as u. The first element is 1/2 go and describes the boundary condition at x = 0. The last element is a constant multiple of H(t) and describes the heated boundary condition. All other elements are zero.

- 31The solution for Equations (26)'may be obtained by standard matrix methods for the solution of a set of ordinary, constant coefficient differential equations (23) This solution is given in Equation (27). - 3C - f (t) t[ 6'LAC3I- A (27) where [X] = the non-singular diagonal matrix of the eigenvalues of A, e []t - the diagonal matrix with elements ekit B = a matrix whose columns are eigenvectors of A. g a column vector whose el.ements are gl, g2, g3,'.e gr for the implicit case and glp g2 9~2, gri for the explicit case. If it is assumed that the heating rate is constant,, the initial temperature is zero, and all eigenvalues of A are real, distinct, and nonzero then Equation (27) may be reduced to Equation (28). (28) -g tAL +t A f'0) ] - A-' f (t) The first term of the right member of Equation (28) is the transient solution and is composed of the siam of exponentials. The eigenvalues of A calculated for r = 2, 3,.,. 7 are all negative and are listed in Appendix D. It can be shown that A is negative definite. This can be done by showing that all even order principal minors are positive and all odd order

-32principal minors are.negative. Thus the transient solution must approach zero and the complete solution must approach the steady state solution represented by the second term of the right member. The heating rate is constant for this case; therefore f(t) is a constant vector and the steady state solution must be a constant. The solution of Equation (28) can be compared with the solution of the continuous Equations (10) which., for constant heating rate and uniform initial temperature distribution) is given by Carslaw and Jaeger(2) as " 2 H L e —Lr 4c4F]t7 3C efe r " co (29) where erfc (y) = l-erf (y) and erf (y) -' e d This equation is given in the (x~%) domain so that the heated face is at x = 0. Thus x = 0 corresponds to x = 1 or x = Xro The temperature u(:0,t) must be compared with ur(t). The temperatures at each of the stations are calculated from Equations (28) and (29) and are compared graphically in Figures 6 for r = 3 for the several. approximations. In these figures the temperature is normalized with respect to the constant heating rate and is plotted as A- oH

-33Equations (28) were evaluated manually be obtaining (1) the inverse of the matrix A; (2) the eigen value matrix; [X]; (3) the eigenvector matrix B and (4) its inverseB-1 The solutions, un, were written as explicit functions of time. These functions of time were then evaluated for specific values of time with the aid of tables of exponentials and a desk calculator. For Equation (29) the temperature at the heated boundary is 2Ht) -- - -X - qu (30) This temperature varies as'j7. However, it was shown above that for the finite difference approximations the solution at the heated boundary approaches a constant temperature which is Thus the error, u(o,t)- ur', eventually would become unbounded if melting did not occur to stop the rise in surface temperature at t = t The time at which melting starts, tl, is attained when the boundary temperature ur reaches zero, This occurs in. Figures 6 when r =go - -. Thus t1 may be determined from these figures by finding H H the intersection of the line - go with the ur curve. When H l-go, the line -, becomes - The minimum value of this ratio is zero when H go-1 go= 0 while its maximum value is unity because __4f> g9 ~

Finite Difference Approximation 1. 0 -- - - Continuous Medium 0.8 1 i 0.6 U2 0.4 0 0.2 0.4 0.6 0.8 1.0 1.2 Figure 6a. Temperatures During the Warming Phase vs. Time Ur Determination Explicit, 1st Differences.

-351. 2....., Finite Difference Approximation 1. 0 -- - Continuous Medium./l 0. 8 1.0 0 0// Figure 6b. Temperatures During the Warming Phase vs. Time Ur Determination Explicit, 2nd Differences. u2 x00 i~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~1 //00 Ur Determination Explicit, 2nd Differences.

-361.2........ Finite Difference Approximation 1.0 Continuous Medium 0.8 I: 0.6 U2 0.4 0.2 0 0.2 0.4 0.6 0.8 1.0 1.2 Figure 6c. Temperatures During the Warming Phase vs. T:ime ur Determination Implicit, ist Differences.

- 371.2. Finite Difference Approximation 1.0 -- -- - Continuous Medium 0.8 U 0.6 0.4 0. 2 0 0.2 0.4 0.6 0.8 1.0 1.2 Figure 6d. Temperatures During the Warming Phase vs. Time Ur Determination Implicit, 2nd Differences.

-38go go If the temperature ratio, g-1 is in the interval, 0 < —- < 1, the -- go-1 - largest error introduced by the finite difference approximations occurs in the temperature, Urn at the boundary, for each approximation. These results have been obtained with r = 3. It is concluded from the small errors that the logarithmic grouping of the finite difference stations in the (xi) domain allows a very accurate approximation of the gradient of u by means of the finite difference approximations. It is expected that an increase in the number of finite difference increments will improve the accuracy. This effect will be checked by means of the computer solutions for r = 6 EitSat the computer solutions for r = 3 will be checked against these analytical results for the finite difference approximations. The choice of H = 1-go as the heating rate allows the solution during warming to be confined to the interval 0 < ~ < 1. The errors introduced by the finite difference approximation have been shown to be small in this interval after the first 0.1 units of time, t, even when r = 3 when ur is determined explicitly. When ur is determined implicitly, the errors are small over the entire interval. Thus this is a good choice of heating rate for the warming phase. It will now be shown that this is a particularly good choice for the steady state melting phase. Analysis of the Melting Phase The transient melting phase cannot be solved analytically in either the continuous representation or the finite difference representation. Equations (25) can be solved analytically, however, for the steady state

-39melting rate when the dimensionless heating rate, H, is constant. This solution can be compared with the steady state solution obtained for the continuous case in Equations (11). Let H(t) be a constant. Then, in steady state dA L a — b, a constant, and Equations (25) can be written as l +21 (b-)jgo_-2 u+Ii +, I (b-I)]uL=., [I +2 I (b-l1lnlz —2 tn[+ 2 (101 ( L)3 L/1 = n- 2,3,,.(. 1) o = H + r U r-, (I - D,' e fe renceS. ) b - H - r + 2r r (2 nd D,-fe rene) If b - 1, then the above equations become - 2 LA t+ LA = LA 4-2 2c1a L4,1 0O H - -rrY- (I r-'D,+Cfererce- ) H - I + 2 r —2 K ar- (2"d 0,'-Ff Ke nc 5) Thus UL, -(% + LA ) (a) =(Un-l L L1-l) (b) (32)

-40Thus each temperature is the arithmetic mean of the adjacent temperatures and the temperature distribution is linear in x. Since ur = 0 the temperature increment between stations can be expressed as There are r increments so go - r,r-l (33 Therefore for first order differences in the d- equation dt H = - go The relationship for the case when second order differences are di used in the dt equation is obtained by using the equation 2 r - 2 go But r LAr- =2 r tAr-I i29 because r —2 2 r — + Ulr O, Therefore and the heating rate of H= i- qO results in a steady state boundary velocity of unity when second order differences are used, also. From Equations (11), if H = 1-goI

-41then b = 1l Thus when the heating rate is the constant, H = 1-go, the steady state solution for the finite difference approximation is identical to the steady state solution for the continuous case in the boundary velocity. It is emphasized here that this does not mean that the boundary position is exact for the approximation because the boundary position depends upon the value of the boundary velocity during the transient phase as well as during the steady state portion. Since the transient phase solution for the boundary velocity is not exact, the boundary position will also be in error. The above solution for the steady state boundary velocity is not affected by the approximations used during the warming phase at the heated boundary nor by the initial. temperature distribution. The problem is now ready to be solved on the differential analyzer. The computer scaling and circuits are the subjects of the next section. Differential Analyzer Scaling and Circuits Equations (23) are scaled for the computer by means of the transformation defined in Equations (34). USA) =1 l~n~t) Go, a' ge (a) (a) Ai -) (1- ~Z ~(b) = (34) (d) (d)

-42In this scaling the unit for the dependent variables in the machine unit and 1 mu = 100 volts.l The transformation of Equations (25) by means of Equations (34) yields Equations (35). Li (T)C0 ^. (a) L......?J r -—: }I- _fI;4 I (b) ).~ y3 ~~i~- ~ Cl I..., -tI, (I' t' J )'. I i _,1 k'i)'rr (C2). L4 0 -. t t- (c3) Ur -'" + " Ur-! -.4 r, 91 /9U- _.r.. (c4)(35) c. (d) lJr-;' -) "[, S "'FT C?(e) (e) _1 - A Q-i 2' rn _, A = 0 7=>^- (g) -) = D,,,;:; -(,: (C.)o r,~lo(c )2 Ii,, - G,n -o=:, 2 (r5-,)-i -f r r i 3c n)rl C(.+) T= 0 (h) 1A familiarity with the electronic differential analyzer is assumed. For a general treatment of this computer the books by Johnson(22) and Korn and Korn(25) are recommended. The paper by Howe and Haneman(20) is an excellent reference on the use of finite difference techniques in solving partial dirferential equations with the differential analyzer. The report by Howe(19 contains an exhaustive study of the solution of the heat equation with fixed boundaries on the electronic differential analyzer.

-43The scaling constants, al, a2, a3, and a4, are chosen to assure as large a variation in voltage at each amplifier output as possible without saturating any of the amplifiers during the solution. This procedure is simplified for Equations (35) by the fact that the maximum absolute value of the temperatures occurs at t = 0. The computer circuits for solving Equations (35) are presented in Figures 7 through 9. The heating rate,'/I(T), shown in Figures 7 and 8 as the input heating rate may be any arbitrary function of time or temperatures or combination of these. This, of course, is one of the prime advantages of the solution by means of the differential analyzer because these functions can be generated and used with little difficulty. The circuits of Figures 7, 8, and 9 were used to solve Equations (35) for r = 3 and r = 6. A Reeves Model 101 computer with 0.1% components was used to solve these equations.

(2n-) Un -Un-1 - 2a4 -1 -U 1 na2 n- 10A 2n2 n Unth Station quation b dX 1 dX Un+ldT na2 Un+ 1 -Un+1 (2n+i)Un+1dX dX dX Figure 7. Computer Circuit for Determining the Temperature at the nth Station (Equation 34b).

1st 0 2nd U( 24 Ur-2 1 Ur Ur 0 < T < T1 r-2 1st-1L- (2r' - 1)T.1 Ur 0 T1 <T < T2 1st - ( )U a4 UsW - + 100V 0 < T < T1 1 2nd - 2r (2r-1)Ur_ 1 - 100V T1 < < T2 1st + r- (2r-l)Ur U 2nd (7r-3)Ur 1 2 a4 r r Cnr __ r 2 1st - ra3 a1a4 Diff a3 2nd - 3r r 0~a4 (a) Implicit Determination of Ur (Equations 34cl and c2) ist 10 1 Ur-2 2nd Ur-2 1st -Ur-l 2nd- 1 Ur U -U 1-0 — 3 Ur r 1st - a 2 ral 1 - / 2nd 2a3 SW Ae-1,Ct___,~a D.R. (b) Explicit Determination of Ur (Equations 34c3 and c4) Figure 8. Computer Circuits for the Determination of Ur.

1st 0 ral Ur-2 2ndra1 Ur-2 1 ra11 1st - Ura2a4 r-1 2nd - 2ralU dX aair- 1 1 d -U ~~~a2a4 d r- 1 33 1st - a2a4 a3 dX 2nd -a__ D.R.6) a2a4 1 4z-F Servo UrSw (see Figure 8) dX*K lo- 0 0 < T < Tj dT dX* YX dX T T < T 7 dr dT Figure 9. Computer Circuit for the Determination of the Boundary Velocity and Position. (Equations 34fl and f2)

CHAPTER III COMPUTER SOLUTION FOR THE SEMI-INFINITE SLAB Equations (34) were solved on the electronic differential analyzer in order to: 1. Determine the relative accuracy of the finite difference approximations; 2. Determine the effect of increasing the number of finite difference increments; 3. Obtain a set of curves with initial temperature as a parameter which could be used to solve a wide range of physical problems~ To accomplish the first and second goals above, the problem was solved for r = 3 and r = 6 with go = gl = g2 = g3 ='' = gr-1 = gr = -1.128 (where gi = the initial temperature at the ith station) and constant heating rates of H = 1 and H = 4. This value of initial temperature was chosen in order to compare the results obtained in this dissertation with Lundau's solution. (27) First and second difference approximatiornswere used at the heated boundary, and ur was determined both explicitly and implicitly. The results of these solutions are presented in Tables I, II, III,. IV, and V. Table I contains values of t1 (time elapsed when melting starts) and b (steady state boundary velocity for constant H) for the continuous -47

TABLE I tl AND b FOR THE SEVERAL APPROXIMATIONS AND THE CONTINUOUS MEDIIUM go = -1.128 Approximations tl(time melting starts) b(steady state boundary speed) Diff. H Case ur Determination ET Determination Cont. Eqns. Comp. Cont. Diff. Eqns. Comp. r=3 r=3 r=6 r=3 r=6 r=3 r=6 1 A Explicit, lst Differences 1st Differences 1.00 1.07 1.076.99 0.47.343.399.34.39 1 B Explicit,2nd Differences 2nd Differences 1.04 1.032 1.00.416.420.42.42 1 C Implicit 1st Differences 1st Differences.99.99.99.343.399.34.40 1 D Implicit, 1st Differences 2nd Differences.99.99 1.00.416.420.41.42 1 E Implicit, 2nd Differences 1st Differences 1.00 1.00 1.00.343.399.34.40 1 F Implicitj 2nd Differences 2nd Differences 1.00 1.00 1.00.416.420.41.42 4 A Explicit 1st Differences 1st Differences 0.0625 0.00 0.00.034 1.88 2.06 1.96 2.06 1.96 4 B Explicit, 2nd Differences 2nd Differences 0.042.04.061 1.88 1.88 1.88 1.88 4 C Implicit,lst Differences 1st Differences.073.07.065 2.06 1.96 2.06 1.96 4 D Implicit, 1st Differences 2nd Differences.073.07.064 1.88 1.88 1.88 1.88 4 E Implicit, 2nd Differences 1st Differences.o60.058.062 2.06 1.96 2.05 1.96 4 F Implicit,2nd Differences 2nd Differences.o60.059.062 1.88 1.88 1.88 1.88

mediunm Equations (8), the analytical solution of Equations (34) and the computer solution. When H = 1 the time, tl, at which melting starts is very close to that for the continuous medium for all cases. However, for the higher heating rate, H = 4, the explicit determination of ur introduces very large errors into the value of tl when only three stations are used. When six stations are used with first order differences the error in t1 is 45%. When second order differences are used this error drops to 2.5%. The analytic solution of the differential equations indicates that these errors are introduced by the approximations and not by the computer. The errors in tl are smaller when the temperature at the rth station is determined implicitly but the error decreases from 12% for Run C, r = 3,to zero for Run F, r = 6 within the accuracy of the solution. When six finite difference increments are used with second order differences there is little to choose between the implicit and explicit determinations of ur as far as the time at which melting starts is concerned. The values of the steady state melting rate, b, were calculated by means of Equation (12a) for the continuous case. The values for the finite difference approximations were obtained by setting all time derivatires of temperatures in Equations (35) equal to zero, letting t = b and solving the resulting set of algebraic equations simultaneously for b. The solutions are polynomials in b of degree r if r is odd and (r-l) if r is even. The coefficients of powers of b are functions of H and go.

-50These polynomials must reduce to b = 1 if H = (l-go) as shown in Equations (29-31). These polynomials are presented below for r = 1 through r = 6. r=l, 1st differences. b - -.e)' =o r=2, 1st differences. (2-g)b -(H+g3ko r=2, 2nd differences. (I-go)b-H= O r=3, 1st differences. b3- ( H+3,) b2 (23 -12 j)b(2- 9 )=D 2nd differences. b3-(II+6 qo b) + (23-2Zo) I-(23Hf +6o)= r=4, 1st differences. (4 -gD) b3-(4 H+9 )b2+(44 -23gb -4(-4Ht 159) = g 2nd differences. (2l-go)b3-:-aH t o)~l-b2i(22-II, )b-(22 l-+;~o)-= 0 r=5, 1st differences. b-(H+ s5)b4O4(23-8ogo) b3-(23oH+43 Oo 0)b + (/689 -,8o ) b -O( 689 H +5255o)= O 2nd differences. b (H+I og ) b'+ (230- 120c0) K-(g3o H 5H D Og)b2 (35) + (1689-8409 ) b -(/6.9 H +450q) O

r 6, 1st differences. ( - o) b- (6 H+ 25o)b +(8 o 0-230ooJ)b3- (5s bH +s9o 5 ) b + (3254-168990) b -(3254H+9'+-sg) 0 2nd differences -(-go52 b)-(3H+ zo~,,) b4+(2o- /sog~) b3 -(~290o-fs0 (o) be + (1627- o- )k0( - i627H+t42o ) = 0 The values H = 1i and H = 4 were chosen instead of H = 2.128 for the solutions in Table I in order to show the effect of the approximations on b. The results in Table I show this error to be rather large for H = 1 where the minimum error is 10.3% even when second order differences are used with r = 6. For H = — 4 when second order differences are used, the value of b is almost exact. These results may be explained by means of Figure 10. In this figure b is plotted as a function of H for go = -1.128. The values of b as a function of H for the continuous medium are plotted as the dashed line. The curve b vs. H for the finite difference approximation intersects the continuous medium curve at H = 1-go, H - 2(1-go), and H = 3(1-go); and in the interval from H = 1-go to H = 3(1-go) is very close to that for the continuous medium. The two curves diverge outside this interval. Thus the close agreement between the values for b when H = 4 and go = -1.128 are predicted. In all these cases the close agreement between the computer results and the analytical solution of the ordinary differential Equations (25) indicates that the computer solution is accurate and establishes a reasonably high confidence in the electronic differential analyzer solution.

-52-- Continuous Medium -—.- First Order Difference 4 l Approximations Second Order Difference Approximations 3 i/ 0 0 1 2 3 4 5 6 7 8 9 10 1I 12 H Figure 10. Steady State Boundary Velocity vs. Heating Rate.

-53Figures 11 through 14 contain sample solutions from the computer results for r = 3 and r = 6 with the second order difference approximations at the boundary. In Figures 11 and 12, 1/2 inch represents one unit of T; in Figures 13 and 14, one inch represents one unit in rF The transient melting solution will now be compared with the solution obtained by Landau for a constant heating rate and uniform initial temperature distribution. Comparison of the Boundary Velocities Tables II - V contain the values of the boundary velocity for the several approximations used in this thesis and those obtained by Landau. (27) The agreement of the differential analyzer results with Landau's results is very good for the cases where r = 6 and second order difference approximations were used for at the boundary (cases B and F)> The effect of changing the approximations used during the warming phase is to vary t1 and, when ur is determined implicitly, to cause the di di error described on page 25 in at t1. This error in d at tL is very apparent for cases C and E, r - 3,in Table II. In these cases the boundary velocity is negative during the first moments of "melting"~ The use of second order differences in the velocity equation dE steepens the slope of the -t curve during the first part of melting when compared to the first order difference curves. Increasing the number of increments has the same qualitative effect. The values in Tables II - V

t) t tt) 6 1.2 5 1.0 4 0.8!(t) 3 0.6 2 0.4 d I 0. 2 0 0.1 0.2 0.3 t 0.4 0.5 0.6 0.7 0.8 Figure 1la. Boundary Velocity and Position vs. Time r=3; H=4; go= -1.128 Explicit Determination of ur. Second Order Difference Approximation of xU Xr

t 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Us - 0.226 -0.4 51 n Figure l1b. Temperature at the Finite Difference Stations vs. Time axX ~XMr

d f. t(t) 6 1.2 4 0.8 2 0.4 0. 00.05 O. 10 0.15 t 0.20 O. 25 0.30. 35 0.40 Figure 12a. Boundary Velocity and Position vs. Time r=3; H=4; go= -1.128 Implicit Determination of u3. Second Order Difference Approximation of _u ax r

0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 - 0.226 -0.451 U3 u2 -0.902 -0.128 Figure 12b. Temperature at the Finite Difference Stations vs. Time r=3; H=4; go= -1.128 Implicit Determination of u3. Second Order Difference Approximation of 6x. xr

dt 6 1.2 4 0.8 U. t 2 0.4 df dt 0 0 _ _ _ _ 0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Figure 13a. Boundary Velocity and Position vs. Time r=6; H=4; go= -1.128 Explicit Determination of u6. au Second Order Difference Approximation of x axr

0 -5 0.10 0.15 0.20 0.25 0.30 0.35 4 - 0. 226 - 0.4 51 U - 0-P77 0.902 -1,128 ~~~~~~~~~~~~~~U2 Figure 13b. Temperature at the Finite Difference Stations vs. Time r=6; Hft4; g = -1.128 Explicit Determinatiof Second Order Difference Approximation of ax r x

dt 6 1.2 4 0.8 _ 0 2 0.4 dt 0 0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 Figure 14a. Boundary Velocity and Position vs. Time r=6; H=4; go= -1.128 Implicit Determination of u6. Second Order Difference Approximation of ~u Xr

t 0 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 -0.226 - 0.451 -0.677 -0.902 Figre b. Temperature at the Finite Difference Stations vs. Time128 Figure 14b. Temperature at the Finite Difference Stations vs. Time r=6; H=4; g = -1.128 Implicit Determination of u6. Second Order Difference Approximation of -x xr

-62indicate that after approximately 0,1 time units of melting only the approximation for in the d_ equation affects This can ~px i a f o dt dt be seen by comparing the values of Td for cases A: C, and E with each other and cases B. D, and F with each other. General Solution Curves for Constant Heating Rate In Appendix E solution curves are presented for initial uniform. temperatures of go = -0~5, -1.0, -2o0, and -3.0~ In each of these solutions H = 1-go, and ur was determined implicitly. Thus these curves can be used to obtain values of the temperature at any of the six stations and the position and velocity of the boundary at any given time. Equations (2), (6), and (8) are used in conjunction with the curves of Appendix E. The temperature curves of Appendix E indicate that the temperature distribution is linear in x in steady state as predicted above. The steady state boundary velocity in each solution is (1.-go) as predicted~ The relatively few values of go for which solutions are obtained in Appendix E do not allow direct solution of the problem for all possible g.. These curves could be used for values of g, which do not correspond to any value of go contained in Appendix E by graphical interpolation. A better solution could be obtained by solving the problem for small increments of go. Solutions for a Time Varying Heating Rate Solutions for a heat pulse of the form H = h(l- cos2itt), shown in Figure Fl, are presented in Appendix F. In these solutions the initial

TABLE II BOUNDARY VELOCITY FOR THE SEVERAL APPROXIMATIONS AND LANDAU'S SOLUTION go = -1.128, H=1 Approximations t( t dt r Case ur Calculation Calculation 1.00 1.02 1.06 1.10 1.20 1.60 2.00 2.40 2.80 3.20 ur..... A Explicit,lst Differences 1st Differences 0 0 0.03.115.24.28.30.32.32 B Explicit,2nd Differences 2nd Differences 0 0.065.13.22.33.37.40.41.41 C Implicit, lst Differences 1st Differences 0 -.05.00.04.11.24.28.31.32.33 3 D Implicit, 1st Differences 2nd Differences 0.065.125.17.24.33.37.39.40.41 E Implicit, 2nd Differences 1st Differences 0 -.03.01.05.12.24.28.31.32.32 F Implicit, 2nd Differences 2nd Differences 0.05.11.16.23.32.37.39.40.41 A Expliitlst Differences st Differences.025.055.10.125.18.26 31.34 36 A Explicit,lst Differences 1st Differences.025.055.10 ~125.18.26.31.34.36.37 B Explicit, 2nd Differences 2nd Differences 0.07.15.19.26.34.37.39.4o.41 C Implicit, 1st Differences 1st Differences 0.04.08.105.17.26.32.34.36 ~37 6 D Implicit, 1st Differences 2nd Differences 0.09.135.165.21.30.34.37.39.40 E Implicit, 2nd Differences 1st Differences -.04.03.07.105.16.26.31.34.36.37 F Implicit, 2nd Differences 2nd Differences 0.085.12.16.21.30.34.37.38.40 LANDAU'S SOLUTION.00.0712.1289.1600.2103.292.336.362.381. 395

TABLE III BOUNDARY POSITION FOR THE SEVERAL APPROXIMATIONS AND LANDAU'S SOLUTION go = -1.128, H=l Approximations (t) r Case u Calculation Calculation 1.00 1.02 1.06 1.10 1.20 1.60 2.00 2.40 2.80 320 A Explicit,lst Differences 1st Differences 0 0 0 0.008.085.190.307.431.554 B Explicit,2nd Differences 2nd Differences 0 0 0.005.023.138.278.433.595 -- C Implicit, 1st Differences 1st Differences 0 -.002 -.002 -.001.006.080.184.305.430.558 3 D Implicit, 1st Differences 2nd Differences 0.001.002.012.032.147.287.439.597 E Implicit, 2nd Differences 1st Differences 0 -.001 -.002 0.009.087.190.308.434.558 F Implicit, 2nd Differences 2nd Differences 0 0 0.011.030.146.283.436.594 A Explicit, lst Differences 1st Differences 0.002.003.008.023.113.229.368.500 B Explicit, 2nd Differences 2nd Differences 0.001.005.012.035.158.299.451.612 C Implicit, 1st Differences 1st Differences 0 O 0.005.020.108.220.348.438 6 D Implicit, 1st Differences 2nd Differences 0.001.006.011.030.032.256.400.547 E Implicit, 2nd Differences 1st Differences 0 0.001.005.018.104.218.346.484 F Implicit, 2nd Differences 2nd Differences 0 0.004.010.030.132.258.397.543 LANDAU'S SOLUTION -- --.00517.0106.0296.1318.263.400.555.705

TABLE IV BOUNDARY VELOCITY FOR THE SEVERAL APPROXIMATIONS AND LANDAU'S SOLUTION go = -1.128, H=4 Approximations d(t) r Case ur Calculation Calculation 000.05.0625.075.10.15.20.30.40.60 A Explicit,lst Differences 1st Differences.60 1.17 1.28 1.35 1.50 1.68 1.81 1.91 1.97 2.02 B Explicit,2nd Differences 2nd Differences 0.24.47.67.97 1.34 1.53 1.70 1.75 1.80 C Implicit, 1st Differences 1st Differences 0 0 0 1.00 1.40 1.54 1.78 1.92 1.96 -- 3 D Implicit, 1st Differences 2nd Differences 0 0 0.16.61 1.19 1.47 1.69 1.76 E Implicit, 2nd Differences 1st Differences 0 0 1.18 1.40 1.46 1.69 1.81 1.94 2.00 F Implicit, 2nd Differences 2nd Differences 0 0.08.35.78 1.28 1.52 1.72 1.79 A Explicit, lst Differences 1st Differences 0.52.98 1.16 1.36 1.58 1.66 80 1.88 -- B Explicit, 2nd Differences 2nd Differences 0 0.20.72 1.16 1.44 1.60 1.72 1.76 C Implicit, lst Differences 1st Differences 0 0 0.66 1.24 1.54 1.68 1.82 1.86 D6 Implicit, 1st Differences 2nd Differences 0 0 0.48 1.06 1.42 1.56 1.68 1.76 E Implicit, 2nd Differences 1st Differences 0 0.20.96 1.28 1.54 1.68 1.80 1.86 F Implicit, 2nd Differences 2nd Differences 0 0 -.10.62 1.10 1.43 1.57 1.69 1.75 LANDAU'S SOLUTION 0 0 0.84 1.17 1.45 1.57 1.69 1.76

TABLE V BOUNDARY POSITION FOR THE SEVERAL APPROXIMATIONS AND LANDAU'S SOLUTION go = -1.128, H=4 Approximations r Case u Calculation d Calculation r ff 0.00.05.0625.075.10.15.20.30.40 A Explicit,lst Differences 1st Differences 0.046.060.074.110.192.280.470.670 B Explicit,2nd Differences 2nd Differences 0 0.004.010.030.100.162.330.504 C Implicit, lst Differences 1st Differences 0 0 0 0.045.130.200.390.586 3 D Implicit, 1st Differences 2nd Differences 0 0 0 0.008.056.124.286.456 E Implicit, 2nd Differences 1st Differences 0 0.oo6.020.056.136.244.410.600 F Implicit, 2nd Differences 2nd Differences 0 0 0.004.020.072.140.304.476 A Explicit, lst Differences 1st Differences 0.006.014.030.062.140.220.394.572 B Explicit, 2nd Differences 2nd Differences 0 0 0.005.030.o96.174.342.512 C Implicit, 1st Differences 1st Differences 0 0 0.006.034.104.184.354.532 6 D Implicit, 1st Differences 2nd Differences 0 0 0 0.020.038.162.324.498 E Implicit, 2nd Differences 1st Differences 0 0 0.008.038.110.192.370.552 F Implicit, 2nd Differences 2nd Differences 0 0 0.002.026.086.166.328.502 LANDAU' S SOLUTION 0 0 0.0074.0328.0998.176.340.506 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~.3.5.06

temperature distribution was uniform in x. The heated boundary temperature was obtained explicitly in each case. In Figure F2 the heating rate, H = 0.5(1-cos23t),was not sufficiently high to cause the boundary to reach the melting temperature. This corresponds to a heat sink. After the heat input returns to zero the temperatures within the body approach go, the temperature at x = O. or x = oo, because there is an infinite heat capacity in the slab. Solutions are also presented for peak values of h = 1, and 1.5 in Figures F3 and F4. In both these cases melting starts at t1 and eventually stops when the heating rate drops to a sufficiently low value. After the heating stops the temperatures again approach the initial temperature.

CHAPTER IV THE FINITE SLAB The finite thickness slab with one dimensional heat flow can be treated by finite difference techniques also, as has been done by Eyres, et al..(13) Otis(34) and Sunderland,(39) and Murray.(32) In each of these references the finite difference increments used were of uniform size. A non-uniform grouping of stations in which the finite di.l rence increments are smaller where the temperature gradient is larger allows a more accurate approximation of the gradient. Such a method for regrouping the finite difference stations for the heating problem with fixed boundaries was developed by Howe.(19) This method is applied to the melting problem in this chapter after the problem has been made dimensionless by means of Equations (2) above and the moving boundary removed by means of another transformation below. The Physical Problem The physical problem is presented schematically in Figure 15. Conducting - - Insulation Medium c H (t) x=O x x =o x = () Figure 15. Physical Schematic, Finite Slab. The slab extends to + oo in both directions perpendicular to x thereby giving one dimensional heat flow. The heat flux, H, is defined to be positive when heat flow is in the direction of positive x.

-69The Mathematical Model The physical problem is described by Equations (36) in the (x,t) domain of Figure 16o S<' C g = a ii rX 06 Z c R (t);t? t O{ t7 (c) tC e (3d) U = g~x) Oc0~lt) O (g) E=X~~~~~~~~~~~ - S t - o ~(h) O- eO, (e0 0 ~ (i) Here tm = the time at which the slab becomes completely melted. All other symbols are as defined on page 9. It is assumed in Equations (36) that there are no heat sources in the conducting medium. tmt __ A Figure 16. The (x,t) Domain, Finite Slab.

-70Figure 16 is drawn for a heating rate which remains negative sufficiently long to cause complete melting of the slab. If the heating rate becomes zero after tl, melting will stop and the slab will approach an equilibrium temperature. If the heating flux becomes positive the temperature will decrease. Freezing cannot occur because perfect ablation has been assumed and thus the liquid is removed immediately upon formation. Thus de/dt < 0 for all time. Equations (36) are made non-dimensional by means of Equatiors (2). Equations (37) result from this transformation. d O C~gS.<t (c) d~ I <. 0 X-' ~:,5 C.) ( d) TO~ - o ~:~"t~).~,.,'t -,% (f) (37) 3 0'XS-5') t-( (g)'A - (,i:) o~,,g~tY -f, 0 (h).~- s o~ ()) The (x,') domain is the same shape as the (x,t') domain of Figure 16.

-71Choice of Scaling Constants K1, K2, K3 The scaling constants, K1, K2, and K3, can be eliminated from Equations (37) by means of Equations (4) as was done in Chapter II for the semi-infinite slabo When the heating rate is constant, Equations (6) can be used to eliminate H in a manner analogous to that used to obtain Equations (7) from Equations (4), Problem Variables Equations (37) describe the problem in terms of the following variables: 1. independent variables, x and t; 2. dependent variables u (Yxt) and t 3. problem parameters: a. heating rate, H(t); b. initial temperature distribution g (x); c. initial thickness, s; 4. characteristics of the conducting medium; a. density function,,0p(X,U); b. specific heat function, Pc(x,%U) c. conductivity function, k(X) The finite slab problem has one more parameter, namely the initial thickness, s, than the semi-infinite slab. It will now be assumed that the thermal characteristics of the conducting medium are constant both in temperature and in x. Thus Tp mic =bk wl b The moving boundary will. be eli.minated i.n the next section,

-72Removal of the Moving Boundary The moving boundary is removed by applying the transformation defined in Equations (38) to Equations (37) to obtain Equations (39). The transformation is (38) t t Also let h ( t) - Xi The transformed equations are: n 2), _ j_ X 2A A S1 oV6 u ~ /\ v, C o ~ -!; o.-'t- (d) cr ) A,-:\ (39) A~~~dA~~~c.~~ <O %, tt~c~ t ~(e) Awt Q~=o;i~~~:,, O x1t o.)t,$-tm (f) -'" 4t( oXj~ ) <>;t = O (g) t- O (h)

-73The (x,t) domain is shown in Figure 17. This figure emphasizes A A the singularity in the problem at t = tm when the conducting medium disappears. tM A tm ti O _.......... x 0 x=l A A Figure 17o The (x,t) Domain, Finite Slab. If the finite difference approximations of Appendix C were applied at this point, the increments of x would be uniform. In order to improve the accuracy of the finite difference approximations the stations will now be grouped closer together near the heated boundary. A This grouping will be obtained by transforming the x variable nonlinearly. The Nonlinear Space Variable Transformation The transformation used to group the finite difference stations near the heated boundary must satisfy the following two criteria: 1. The region near x = 1 must be expanded. 2. The Jacobian of the transformation must not vanish in A the interval 0 < x < l1 For convenience it is also required that, if x is the transformed variable, its range must be 0 < x < 1.

The nonlinear transformation is defined in functional form to be. Bz -t ( >) O6XCI;, O6 X A (a) ^, ^ (40) t - t o't-tt, (b) The Jacobian of this transformation is -.I A / 1- 2,-f (X) This Jacobian must be positive for all. values of x in the interval 0 < x < 1. This requirement is justified by a consideration of Figure 1.8 and the two criteria above. x=f(x) / A / / // X ~ = O x Figure 18o The Nonlinear Transformation.~ A A A In order to expand x near x = 1, the function x = f(x) mu.st pass throu.gh the origin and through the point (1.1). Therefore f'(x) must be positive somewhere in the interval if f(x) is continuous. Since f' (x) f 0 in the interval 0 < x < 1) then f' (x) > 0 throughout this interval.

-75Equations (39) are now transformed by means of Equations (40) to obtain Equations (41) in which H(t)= i (t) - =.'') - d- -X +'g b ) -'T' ~ 3+ O O~LrtctO (a) H +t) =- fl u ) Lg i m='1; o' t -sty (b), ~-=~~~~O ou~~~o.stt, (c) U4O xAl; 00Q~ t, (d) A- ~~~~~< <~ tt, tct2 (e) (41) LA =g(x; t O (g) 2 p~~~~~ ~s~ t_~tc D (h) k=03L~~~~~~~o ~ X o o t Mb (i) The Equations in Finite Difference Form Equations (41) are in a form suitable for applying the finite difference approximations of Appendix C. When these approximations are introduced into Equations (41), the results are Equations (42).

-76F;n (,Ao r2 I r3Y u' i.t -() t2, r(b) L r - H r — r- 4 34 Fr (t) (t) o s- t (clt, rA ~ f Lr 3J r-f b r, J (e) = V o z t c t, (d) =h-t f r^ I) 4t < tlt2r (e) (42 ) ur.<0 o tt't, (g) IUr - 9O, n = O I (rI t= /0- (h) n= o, 1,.r f C o (C2) t=s t=C (i)

-77Second order differences have been used in Equations (42) at the heated boundary because the results of the semi-infinite slab clearly demonstrate the improved accuracy of the second order over the first order approximations. In Equation (42-C1) the boundary temperature ur, is determined explicitly; in Equation (42-C2) it is determined implicitly. Scaling for the Differential Analyzer Equations (42) are scaled for the differential analyzer by means of Equations (43). U -()r - di LAt) nf=O, r. (a) AT) (Ct) (b) 1V(') (t) ( (43) 1- = CI4 r t (d) Equations (43a, b, and c) are identical with Equations (34a, b, and c). In Equation (43d) r2 is written explicitly in order to eliminate r2 from Equations (42a and b). When Equations (43) are introduced into Equations (42), the results are Equations (44) which are the computer equations.

-78~d r = o0 ( +2 h');r4 fo -I2) 2 a4- 7 X1r) O$I (a) u/unf i (f -r-+f,, / d UU (-, /_ _-' / tJ rr dd A? +U dA Ur o U', - +,\;(-3") (b) Dr=' r _ r, — r — /-( O'', (_) 4-r/' " — I r — (c2) 4- 2 zrc4 2 ) ffrr )a + oyrir, f - < f 4.2 f )4 \> (44) dXC1 il Uf' " Xf-r-)p X e) (d) -F~~~~~~~~7~~~~ =e3~~~~~~~ O <(e) U~ < 0 o??, (f) UrO =T~T ~ (g) U, G - o, I,.,, r-1o ~, (h) n O rli,,r 4or (C x z S.-,

-79The computer circuits for solving Equations (44) are given in Figures (19a-d) on the following pages. A Specific nonlinear x transformations must be selected before solutions can be obtained on the differential analyzer. Two possible choices of this transformation are presented in the next section. Nonlinear Space Variable Transformations Two nonlinear space variable transformations are considered here. The first is a power transformation; the second an exponential transformation. The power transformation is given in Equations (45).;-, - -- J __c) (a) - (b) (45) C and k are arbitrary parameters, In Equations (45) the interval of x is 0 < x < 1 as can be verified by direct substitution. The Jacobian, J = dx/dx of the transformation is >.c (, I C O- = — (l-c) o A s I and f 4 __ ___C _ —) O5 ~

1 1 Un-1 Un- l f1 1f2 a-a. X 1 - uni o-f'n' ax _ __ U rfn dX Un-l dX nl 1.1 dT X) cr 1 1 n I U n UnUn -- 1 n U 1 U NO1'E: 1. Pote ntiometer settings are given for n 1,2,.. 2. For ntO: Pots 1, 2, and 4 are set at 0. Pots 3 and 5 are set at 1/f t(i/f 2 + 1 3. 100V 1 machine unit. Figure l9a. Computer Schematic for the nth Station.

Ur 2 1 1 Ur.-2 ri a2a 4 1 Ur-2 Ur2 1 Ur 7O 1 Ur-V r r 2 4 X' 1 -Ur Ur -2.25 U 4 r- 1 1 U _ r 2f X2 f I f t 1 fft a s l a a r v h2 f'~- f1X 1 rA"S4 U/ 1 a33 331 A rfsa1 2 - Implicit Determination of Ur. Explicit Determination Of Ur NOTE: Ur Ur 0<T<l Figure 19b. Computer Schematic for Forming Ur. U* U U*~~~~; U, r ~~~~ ~~~~D. R. r D.R. NOTE: U w+100V whenU <0 i -l0OV when U,>0 D.R.w Differential Relay. Figure 19c. Switching Circuit for Ure

+Ur-2 1 a U_ 1 _U' 2 aa l dXI' + _ +-* T - 2rfaIa r a2I4 dX1* dX + + ra 4 Servo dC) - 00 < T < T1 dT - dX dX* d- d- T < T< T2 Figure 19d. Computer Schematic for the Formation of dX/dT and X(T).

-83Thus the Jacobian is positive for C in the interval 0 < C < 1 and requirement 2 above is satisfied. Also, if k > 1 and 0 < C < 1, then d2./dx2 = dJ/dx < 0 and requirement 1 is satisfied also because the curve x =[1/l-(l-C)k _l-(l-Cx)k] is convex upward. When Equations (45) are used with k = 2, a uniform finite difference station distribution in x causes the finite difference stations in x to be distributed parabolically. This is shown in Table VI. A second nonlinear transformation is the exponential transformation given in Equations (46). (a) (46) -t =t (b) A Again the interval in x is 0 < x < 1 for x in the interval 0 < x < 1. The Jacobian of this transformation is l_ i- _ Ae Thus J> 0 for all k > 0O Since d2x /dx2 - k2e -kx/1-e k the curve x = 1-e-kx/l-ek is convex upward in the interval 0 < x < 1 ard must lie above the line x = x as required. The station locations for a six station problem are shown in Table VI for the power transformation, Equations (4-5), and the exponential transformation, Equations (46). The values of f' n = dx/dx are also shown at each station. The stations can be grouped more closely near the heated boundary by increasing C and k in Equations (45) or (46)~ While the approximation of the temperature gradients improves in these areas the increased

-84accuracy is not without its price, As the stations are grouped closer together, the temperature difference between stations decreases. In the computing circuits the differences are taken in the summing amplifiers and depend upon the coefficient potentiometer settings on the inputs. Thus as the differences become smaller the errors in the coefficient potentiometer settings become larger by comparison and eventually limit the accuracy of the solution. The same inaccuracies become predominant also when very large numbers of finite difference increments are used to improve the accuracy of the solution, TABLE VI STATION POSITIONS AND dx/dx FOR SIX STATIONS LINEAR, PARABOLIC, AND EXPONENTIAL STATION DISTRIBUTIONS n n 0 1 2 3 4 5 6 xn 0 0.1667 0. 3333.5000.6667.8333 1 0000 xn = - 0 0,2414 0o4530 46346 78653 9081 1.0000 f'= 14(.1-O7x) 1.538 1-346 1.179 1.OC00.821.641.462 n. 0.91. 1-e'x x le 0 0.2)428.4483.6225.7698.8945 1o0000 0.37 -X f, = e 1o 582 1.3392 1,134.9595.8122.6875 0o5820 n 0~37 Computer Solution Results The finite slab problem was sol.ved by means of the same electronic differential analyzer used for the semi-infinite solid solution above with. six finite di.fference incremrents. A linear distribution of stations i.n x

- o85A 1-(1-Cx)2 and a parabolic distribution) x with = 0.7 were used. 1-(L-C)2 Samples of these solutions are presented in Appendix G. In each case u6 was determined explicitly. The time tl, at which melting starts,and the time, tm, at which melting is completed when the heating rate is a constant are shown for the linear and parabolic station groupings in Table VII for three thickness, s = 0.5, 1.09 and 2.0. In all these cases the initial temperature was gn = -1.0 and the heating rate was H = -(l-go) = -2.0. The time, tm; obtained from the computer agrees very closely with that obtained from the continuous case. The value to tm for the continuous case was obtained from the equation t=,[(g-l)/H] ~(o) which is derived in Appendix H. The errors in this value are about one per cent. The value for the time at which melting starts is obtained for the continuous medium by solving Equations (39) during the interval 0 < t < t1 in terms of a Fourier series. The results of this solution are given here for a constant heating rate and a uniform initial temperature. A AA~~~h I Hx, TUMI t +' Z n''Z X Ccu~nL,(47) The value of t1 is obtained from this equation by solving it numerically. The values of t1 obtained from the computer and from Equation (47) are compared in Table VII, When e(O) = 0.5, these values agree

TABLE VII ti AND tm FOR SEVERAL THICKNESSES ANT STATION GRUPINGS H = -2; gn = -1.0 t1 tr A Continuous Computer % Error Continuous Computer lo Error 0.5 A n/r 0.167 o.164 -1.8 0.500 0.494 -1.2 0.5 1 -t, 7n 2 0.167 o.164 -1.8 0.500 0.494 -1.2.91 1.0 4 = n/r 0.19598 0.178 -9.2 1.000 0.986 -1.14 1.0 1-(l-.7xn)2 0.19598 0.199 +1.5 1.000 0.983 -1.7 2.0 AX n/r o.1963 o.16o -i8.4 2.000 1.95 -2.5 2.0 1-(i-.7xn)2 0.1963 0.232 +i8.4 2.000 1.98 -1.0.91 Yn n/r 0.1963 0.196 0

-87within less than two per cent. The regrouping of stations does not affect appreciably this value. When (0O) = 1.0, the percentage error in t1 increases to -9.2%o. Here regrouping the stations lowers this error to 1.5%. When (0o) = 2.0, the error in t1 becomes -18.4%. Regrouping the stations makes the error +18.4% which is still large. Grouping the stations more closely together near the heated boundary might decrease this error still further. However, when the value of t1 from the continuous case (0.1963) is compared to t1 for the continuous semi-infinite case (0.1963), from Equation (29), it is apparent that the slab with e(0) = 2.0 behaves as a semi-infinite body during warm-up. This is further corroborated by Figure G5a and G6a in which u0 does not vary until after melting starts. Thus the semi.-infinite formulation could be used with good accuracy for this portion of the problem at least for initial temperatures go < -1 The decrease in accuracy with an increase in initial thickness was expected because the same number of stations was used in each case. As the initial thickness increases the accuracy of the gradient approximations must decrease because the error in the finite difference approximations varies as the square of the increment. The value of u6 for the continuous case is plotted on Figures G3a and G4ao These values were obtained from Equation (47). The improvement in accuracy with the regrouping of stations is very apparent when these two figures are compared. The maximum percentage error when the stations are linearly distributed in x is 22%; for the parabolic distribution it is 10o35o

-88The effect of the station regrouping is to make the differences between adjacent stations more nearly uniform compared to the linear station grouping in x. This can be seen by comparing Figure Gla with G2a, G3a with G4a, and G5a with G6a. A possible desideratum for station grouping could be that the temperatures be linearly distributed over the stations in the part of the solution of most interest. Figures G7 and G8 are solutions of the melting problem for a heat flux of the form H = 1 - cos 4itt. It can be seen that with this heat input melting starts at t1 and ends at t2. The accuracy of the solution can be checked by considering the energy exchange during the heating. The total heat input is @ - 0 ( I - C OS Z4 T t at - From Figure G7a the equilibrium temperature is -0.048. The rise in temperature was thus 0.952 and the heat required to raise the remaining part of the slab to this temperature was Q - 0,952 ( -o 3) - 0,4+40 From Figure G7b, 0.036 of the slab was melted~ The heat required to raise this much of the slab to 0 was Q 0 o 0 36 The heat required to melt 0.036 of the slab is Q- = O, 03 The total heat required is then Q = O, + Q+ Q3= 051 The error is O -(O, + t +O,) _ O, 06 Q This error could be caused by inaccuracies in reading the curves.

-89When Figure G8 is compared with Figure G7, the effect of the station regrouping is seen to be to spread out the temperature curves and make the temperature increment between stations more nearly equal as noted above. In this chapter the finite slab has been analyzed and solutions obtained on the electronic differential analyzer. The results are not as striking as those obtained with the semi-infinite solid, but it has been shown that the accuracy of the solutions can be improved appreciably by regrouping the finite difference stations and that, for slabs of such thickness that the finite thickness formulation gives unacceptable errors, the semi-infinite formulation can be used to reduce the errors during the warming phase and the first part of melting.

CHAPTER V CONCLUSION The problem of the diffusion equation with a free boundary has been analyzed in this dissertation for both semi-infinite and finite thickness conducting media in which diffusion proceeds in one dimension. The object of the analysis was to render the problem into a form which could be readily solved on the differential analyzer and which would yield accurate solutions with a modest amount of computing equipment. It is believed that these objectives have been attained. Description of the Solution The problem equations were first made non-dimensional. The particular dimensionless ratios for the time variable, the space variable and the heating rate which are derived in this dissertation are believed to be original. In the process of eliminating the dimensions, scale factors were introduced which allow physical temperature, time and space variables, and physical heating rates of any range to be related to any convenient non-dimensional range. It has also been shown that, for a constant heating rate, these scale factors can be chosen to eliminate the heating rate. This reduces the semi-infinite problem to one parameter and the finite slab problem to two. The effect of varying these scale factors is to expand or contract the space and time variables in the non-dimensional variable domain. For the semi-infinite solid the infinite space variable was reduced to the unit interval. This permitted the approximation of the derivatives of the temperature with respect to the space variable by means of the finite differences over the entire interval -90

The finite difference appr,'_oxi.mati.on reduced'the problem, to a set of ordinary differential equations which could be solved on tihe differential analyzer. It was shown for a ut'niform initial temperature distribution and a constant heating rate that during warminng, the trari. — sient solution of these differential equl.ations is a suli of decayi ng exponent'ials plus a constant at each station), while the so.lu.tion of the continuous case proceeds as the squ.a.re root of time at+ the he, <ated sur, - face. Thus the error between the conti.nuous and fi.ni'te difference solu-, tions at the heated boundary could eventuall~y become unrbonded. However it was also shown how the dimrlensilonless heat"ing:ra't.e could he chosen to keep this error very smrall until Tnelting started, F.l-,r-ther i t was sJto'wn that this particular choice of constartit heai:,_,-Lrg rate caused. tIhe s,eady' state melting rate o&tained from.a the fin -it.e dIftter'ence eqilations to be exactly the same as that f'or the CenJtj rlv.o,.s case,(: Two methods fo)-r dete-:i ri. ntg'the K'>.a,u:td hon ulrir(-.a../ triTpeb. r-tare during the warming phase were dev-e.lopedo I['The imp.lci.tc met-ehod has th'Ie advantages of (1) allowLng the i.rn:Ltia.l a1bo)ur",da.r>t te.rrperaatture to ue arl:i_trary and (2) being more accurate than the e-xpl,.i,Lci.t method~ The explicit method (1.) forces the tjoundary v'elzorlty to be zero w.ler. m?:lting star'ts and (2) allows the heat balance -at th Kir:,,:j:A.dlary to,-, tbe d.etermined dunl rrilg melting in order to determine whenr the heatirng rate decreases sulfficientl.y to stop melt.ingo The semi-infinite slab prtoh-lemr; was solved orn an electt;ro:nJic differential. analyzer for three and for six.fini te di:.fferelnce incremrenrtso The results show that reasonable accu.racy can. te obtained'wir;t1th six i.ncrementso The one-parameter problem for constant heating rate and

92uniform initial conditions was solved for several initial temperatures~ These solutions are presented in Appendix Eo The problem. was also solvedfora heat pulse of the form (l-cos25tt) to demonstrate the ability of the electronic differential analyzer to solve the problem with arbitrary heating rates. Samples of these soluti.ons are included in Appendix F. The finite thickness slab solution was mnade more accurate by grouping the stations closer together near the heated boundary'. It was shown further that when the slab became too thi.ck to obtai n accurate results with only six finite difference iLncrements, the semi-infini.te slab formulation could be used during the warmi ng phase to obtain accurate results for a constant heating rat;e. This problem also was solved for a heat pulse of the form (1.-cos4,Tt) and sample solution curves are presented in Appendix Go Extensions of the Soluti.tion The problem considered in. thi.s dissertation has been described in. terms of heat flow and the melting of a solid~ The solution obtained i.s directly applicable, in its non-dimensional form, -to any problem described by the diffusion equastion with a free boundary of this type. Thus it is directly extended to problems of mass diffusi on, crystal, growth., sublimation, and neutron diffusion and the control. of nuclear reactors. In each case, of course, one dimensional di.ffusion must occur and the equivalent of perfect ablation must occ')ur. Also, suitable non-dimensional ratios would have to be derived to relate the non-dimensional solu.tion to the dimensional. probl em O

-93The solution obtained in this dissertation is limited to the case of perfect ablation. It can be extended to the two phase problem in which the liquid and solid phase remain in contact at the free boundary by writing the equations for each phase. The space variable must be normalized with respect to the liquid phase thickness in the equations for this phase and with. respect to the solid phase thickness in the solid phase equations. The finite difference approximations could then be applied. Care would have to be exercised in approximating the gradient at the free boundary, where it is discontinuous. This difficulty could be overcome by using forward and backward differences at this point. Care would also have to be exercised during the warming period and just after melting starts, because the liquid thickness would be zero during these times and the equations have a singularity when the thickness is zero. This could be overcome by starting melting with a small, non-zero liquid thickness and assuming that the melting rate would be linear just after melting starts. In order to solve the two-phase problem it must be assumed that there are no convection currents in the liquid. This is a very restrictive assumption and is unrealistic after the liquid thickness becomes large enough to allow free motion of the fluid. This restricts the validity of the solutiono The one dimensional heat flow problem has been described in this dissertation in a Cartesian reference system. It can also be described as radial heat flow in cylindrical or spherical coordinates. The same methods could be used in the latter two coordinate systems.

The extension of the results of this dissertation to two and three dimensional heat flow is complicated by having to determine what occurs at the "corners". No attempt has been made in this dissertation to solve this "corner" problem~ The solution obtained here does provide a set of useful. dimensionless ratios for the melting problem; provides for the selection of an optimum heating rate for the semi-infinite slab when the heating rate is constant; demonstrates inaccuracies inherent in the finite difference approach and developes means for avoiding large errors due to these inaccuracies; provides a means for solving the melting problem which yields very good accuracies with only a nominal amount of computing equipment; refines the solution of the finite slab by regrouping the stations and shows that the semi-infinite solution may be employed for slabs sufficiently thick to cause large errors when using only a few finite difference increments.

APPENDIX A DIMENSIONAL ANALYSIS Equations (la-i) are expressed in terms of p, c, k, L, x, t, H, and u. These physical quantities must have units of mass (M), length (D), time (T), and temperature (G) or combinations of these four basic quantities. The eight physical quantities and the four basic dimensions then allow, according to Buckingham's T theorem of dimensional analysis,' the formation of four dimensionless ratios. The problem can thus be expressed in terms of four dimensionless variables. The problem will be analyzed in terms of x, t, u, and H. That these form a dimensionally independent set may be shown by forming the product m M (atb [ ije) MVi O -0 0 o Now d; m () - T 3 oi,' (Fi)= VITThus D ITbM/V T-de - M/ )~'T o 1 A complete discussion of Buckinghal's theorem can be found in the books by Durand(9) and by Huntley. (21) -95

-96The following set of linear equations is obtained by setting the exponents of like factors equal. = o b- d-o C =O d- O Clearly the only solution is the trivial one thereby proving the dimensional independence of the chosen set. The dimensionless ratios are determined by forming the product of p, c, k, and L, each raised to arbitrary power, with the first power of x, u, t, and H. Thus four dimensionless ratios Tt1 A22 r3, and A4 are formed. The first of these is derived here. Let iT -' I b)-c Le Then dim (r ) MD= M n O 0 - d e-) aCT 3 2b 2-rb9-bl~g027 0 D T Equating exponents of like factors yields: a+o2- O ()) -3aQ+2 +d2e+ I O (0) -2b-3 -2 - o () The solution of these equations yields Er,=PcL iS

-97Similarly E L U In order to make Equations (la-i) dimensionless, the following dimensionless variables are defined, K-K At ~K~aC. In these definitions K; is a non-dimensional constant of arbitrary value. Its introduction here allows the non-dimensional variables to be chosen in any convenient range and still be related to any given set of dimensional variables. The numerical values values of the non-dimensional variables are presented in Table A-II for several common materials whose physical properties are listed in Table A-Io In all cases the physical properties are assumed to be independent of temperature. In many ablation problems heating rates on the order of 102 to 103 are of interest. For these values of H, a convenient value for K would be 10-4.

-98TABLE Al PHYSICAL CONSTANTS FOR SEVERAL COMMON MATERIALS1 -gm - cal cal cal C Material p -m k _c a C L uc cm3 ggm sec Oc gm C gm Aluminum 2.7 1.01 0.25 76.8 658 Copper 8.9 0.858 o.o 096 42.0 1083 Ice 0.92 o.oo004 0.50 79.7 0 Magnesium 1.8 0 376 0.28 70.1 615 Silver 10.3 0. 992 0.07 21 o.1 961 Stainless Steel 7.8 0o107 0.12 48(iron) 1475 Water 1.0 0.0014 1 0 79.7 0 1All data are from the Handbook of Chemistry and Physics, 36th Edition, Chemical Rubber Publishing Co., Cleveland, Ohio, (1954).

-99TABLE A2 DIMENSIONLESS VARIABLES FOR SEVERAL COMMON MATERIALS Material xx 10-4 t x 10-9 u x 10 Aluminum 3.78 x 2.14 t 3.26(u-uc) 8.50 H Copper 4.17 x 1.75 t 228(u-u) 6.37 H Ice o.663 x 382 t 6.27(u-uc) 23.7 H Magne s ium 7.25 x 3.94 t 3. 99(u-uc) 14o6 H Silver 21.6 x 6.41 t 3.32(u-u —) 155 H Stainless Steel 39.2 x 17.6 t 2.87(u-u ) 5.97 H In this table x is measured in centimeters t is measured in seconds u is measured in ~ Centigrade H is measured in gram calories/square centimeter-se cond K1 = K2 = K3 = 1

APPENDIX B STEADY STATE SOLUTION FOR THE SEMI-INFINITE SLAB WITH A CONSTANT HEATING RATE In steady state the temperature distribution in the (x,t) domain becomes a function of the space variable, x, only when the heating rate, H, is constant. This provides a means for solving Equations (8) for the steady state boundary velocity and the temperature distribution. Since the temperatures are functions of x only O.; = 0 (Bl) The boundary velocity, de/dt, must be a constant, also, because dQ = H - e ~= I (B2) Here H is constant by hypothesis and -4- is constant because it'DX X I is evaluated at x=l. Therefore de/dt = b, a constant. The steady state problem now may be statedXv, 7[-7-b b m XO O2cl (a) u A - (b) (B3) LA o (c) Expand (B3a) to obtain a" -..Xa -1 b00 d I, -100

-101This equation may be written x *d l+- -b)t The second factor is now equated to zero to obtain x 01'm (I - b) = o (B4) Equation (B4) may be integrated once to obtain a first order differential equation. Thus: a — _U -( -b)( - C, or x LA-b L = C (B5) where C1 is a constant of integration. The solution of (B5) is C rC;x S d y C, lr u xb But Thus - = C -b + C, b { 4 d A C2 +a, (B6) 1/1=C3_Xb t- -

- 102The constants are evaluated by means of the boundary conditions (B3b) and (B3c). Thus and C2 -go The solution, with constants evaluated, is LA (I- Xb) go (B7) Equation (B2) provides a means for computing b: (B8) Equations (B7) and (B8) indicate that if H = 1-go then the temperature will be linearly distributed in the x coordinate. This property will be used to demonstrate that the finite difference approximation solution in steady state is identical. to the continuous mediim steady state solution for this heating rate.

APPENDIX C FINITE DIFFERENCE APPROXIMATIONS AND ERRORS When Equations (5) of Chapter II are transformed by means of Equations (8), the results are Equations (C1). q <M )i G(;)tA );*=i,t, <) i'()s < X, LA 4A ( tA o S ~ octst, (c) LA ~o ~_i; o~ctct, (d) (C1) X~~~~~> o ~tctbt (e) iX =z O X;V', ta't (f) 1a-gr =9~O~ t> (g) LA- g on s, = (h) =~ too (i) In these equations the derivatives with respect to x are to be replaced by finite difference approximations. The equations are written in such a manner that only the finite difference approximation of the first derivative is required. These approximation formulae and the errors introduced by their use will be derived for a general function, G(x,t), and then the formulae will be applied to Equations (C1)o -103

Central Difference Approximations To obtain the finite difference approximations the conducting medium is divided into r increments as in Figure C1. ~earJ i I I 11 U,,I I I lI n 0 1 2 n-l n n+l r-2 r-1 r Figure C1. Finite Difference Increments. The function Q(x,t) is first expanded in a Taylor series in x about the nth station. Thus i9 (p<,-Y- +-en, t) =r (t) t |) 2G w X/c< lXn/ Ym! +'(C2) By means of Equation (C2), G(xn+l,t) can be evaluatedwith Thus &4, =+ & ~ 3 t) 0n _t t Icn 2 L' DXlxr 31 and o,-l=@7i -0* ~-/~ti-K4- 9 X14 3!

-105on and i-| I and all even order derivative terms are eliminated from these two equations by subtracting the second from the first. Thus Anal ZX + b =,nx+2 -3 + "' The approximation for is obtained by solving for 219 from this equation and neglecting all terms of higher order than the first in (Ax). The required approximation is D 0 0.- ct I-ep I (c3) ~__ )_e e, - 6n_, and the error introduced by neglecting the higher order terms is - Bl. X3!'a x,5lx.5!'3 If Ax is small and all derivatives exist as required by the Taylor series expansion, the error will then vary as (Ax). Equation (C3) is the central difference approximation of the first derivative of a with respect to x evaluated at the nth finite difference station. It requires that the values of G at each adjacent station be known. At the right hand boundary, the rth station, it is impossible to evaluate G at a station to the right. Thus backward differences must be employed to evaluate X r a Backward Difference Approximations The function, G(x,t), is now expanded about the rth station (refer to Figure C1). Thus.. {,.,, A i,2-~ f. a)_1

Then l-~ may be approximated by means of the values of G at the rth and (r-l)st stations. Solving for - from these two equations yields i —Xr __ X I - - 1 I (C4) Equation (C4) is the backward difference approximation of when all terms of degree greater than unity if (Ax) are neglected. The error introduced by this approximation is -~,' z>1 K s' In this case the error varies as (Ax) for small Ax. In Equation (C4) it is required that the value of the function 9 be known at two pointso This equation employs first order differences onlyo A more accurate approximation is obtained by using second order differences in which the values of Q at the rth, r-lst, and r-2nd stations are required. These values are A;)) I It'~Yu SXY 41 I -3-1 @l -( - ___ Gr-'9~~ X ~ i, 258t<( - < r m} <

-107When 60r and ~~*~ are eliminated from the right members of these Xr three equations the result is 3 ir2 IL9- = _Dr3_ _Lx( ) This equation is solved for r to obtain dl <<r=> Z Lo -2- +2 sr- + T3 rx() - 2 The approximation of, is thus LZ [r2 89r29~~ r-2] (c5) The error introduced by this approximation is This error varies as (Ax)2 and thus is of the same order of magnitude as the central difference error of Equation (C3). Equations (C3), (c4), and (C5) will now be used to eliminate the x variable from Equations (C1). The Finite Difference Approximations of the Problem Equations The space variable, x, is eliminated from Equation (Cla) for the stations n=l through n=r-l by means of the central difference approximation of Equation (C3). Equation (Cla) is X. LA _..-. I D L..A

-1o8The first term of the right member is considered first. Let Then, by means of Equation (C3) But bX = 7, X M -I rtx _ 2Thus The second term is approximated by JdCrL~ Y l+l - Uri) since Thus Equation (Cl.a) is approximated by Equation (C6)o ~hG s k~s-),a t~o (ah (ntei te. t4 n= -2d t it +K( +i) n+ (C6) This is Equation (25b) of the text when ) 1

-109When the heated boundary temperature is determined implicitly during warming, Equation (Cla) is used to obtain dur/dt. During warming de/dt = 0 and the second term of the right member of Equation (Cla) is zero. Also, from Equation (Clb) Let Then -~ I_3,r r J ~],a.,r / -' C'_) X2i' ( ] Here Ax = 1/2r and it is possible to use the half interval to improve the accuracy of the approximation. Thus Equation ( reduces to Equation (5cl) of the text when Equation (07) reduces to Equation (a5cl) of the text when Jo = =c =k = 1. When second order backward differences are used to determine the boundary temperature implicitly let Again H = ~)l) }&r In this case the following approximations are used >A | ~ LasArr- and c Ir, - " a/r

When these four relationships are substituted into Equation (Claj with Equation (C5) in. which Ax = 1/2r the result is Equation (C8). 3 r _ __ (I) (r- I h r- l uD r )- 14(2r-l r') (c8) Equation (C8) reduces to Equation (25C2) of the texto Equation (Clb) is approximated with first order differences by means of Equation (C4) to obtain ur explicitly. Thus Sit(t) ~=~~ j 4sr d ( r (c9) Equation (C9) reduces to Equation (25C3) during warming when 2/= 0 if I It reduces to Equation (25fl) during melting when Ur = 0 When second order backward difference approximations are used in Equation (Clb) the result is Equation (ClO). H(t)=,,Y d r( LA ae lA + r-2) (C10) Equation (Cl.0) reduces to Equation (25C4) during warming and to Equation (25f2) during melting if )op = = I o Additional Approximations for the Finite Slab To approximate the gradient at the insulated boundary, x = 0, a virtual station is placed at x = -nxo The temperature at this virtual

-111station is required to be such that Thus LA = -l (Cil) This boundary condition approximation could have been obtained also by ending the conducting medium at x = 1/2Atx and then forcing U0 = U1 The approximation of dur/dt during warming to obtain the implicit formulation. of ur requires that Equation (40a) be written as _ -0 I I __ID - - - ~S; f L6ffi)'O' a ) -ffi-, During warming, of course, d,/dt is zero. The quantity, is the heat flux at any point in. the slab. Thus if F represents the flux and if a second order backward difference is used to approximate aF/6x at the heated boundary with a finite difference interval of Ax/2, the result is Butc= X[r 3 F +2 E iF - Fr. But F -H ___ r- wr — _r

-112Thus - ~ 7- — 3 r H -4 + L-/- ___ rr~r _iEL. Lp r (C12) L'fi r-~ {f r- l L t) Only second order differences at the heated boundary are used for the finite slab.

APPENDIX D EIGEN VALUES OF THE MATRIX A The eigen values for the matrix, A, of Equation (26) are listed below for r = 3 through r = 7 for the explicit determination of Ur and for r = 2 through r = 6 for the implicit determination of uro o. Explicit method, First differences ur = ur.l +(l/r)H n r = 3 r = 4 r 5 r = 6 r = 7 1 -0.32056 - 0.25155 - 0.21413 - 0.19008 - 0.17310 2 -4.67945 - 3.17196 - 2.42348 - 2.15232 - 1.90842 3 - -14.07558 - 9.91283 - 7.97425 - 6.82072 4 - - -29. 34803 -21. 31184 -17o 34784 5 - -50.87021 -37.86498 6..- -78.88755 2. Explicit method, Second differences ur = 4/3 ur_ 1/3 u,,_2 + 2/3r H n r = 3 r 4 r r =6 r = 7 1 -0.215370 - 0.20967 - 0179190 - 0.176217 - 0.16326 2 -3.11963 - 2.54390 - 2.19925 - 1.95634 - 1.77689 3 - -11.12466 - 8.48754 - 7.13172 - 6.26940 4 - - -25.12152 -18.90817 -15.80432 5.. -45.16184 -34.35316 6..-71.63703 -113

-1143. Implicit method, First differences dur/dt = r(2r-1)url - r(2r-1.)ur + 2rH Xn n r =2 r = 3 r =4 r = 5 r = 6 1 -0.39445 - 0.28273 - 0.231.74 - 0.20160 - 0o18132 2 -7.60555 - 3.80059 - 2.81687 - 2.32581 - 2.02428 3 - -20.91580 -11.67256 - 8.85285 - 7.35941. 4 - -41 27602 -24.47488 -19.15244 5 - - - 68~86556 -43.47948 6 - -103 79791 4. Implicit method, Second differences dur/dt = -r/2(r-1)ur_2 + 2r(2r-1)ur1i - r/2(7r-3)ur + 3rH Xn n r =2 r = 3 r = 4 r 5 r - 1 -0.31534 - 0.27128 - 0~22794 - 0~19985 - 0.17792 2 -6.68466 - 5.67350 - 2.76862 - 2.30338 - 1.97764 3 -33.05158 -.11.69403 - 8.80312 - 7.16368 4 - - 63.30421 -25.16050 -18.64387 5 -1103.52709 -42.78379 6 - -119 o24 925 All the eigen values that have been calculated for the matrix A, negative,

APPENDIX E GENERAL SOLUTION CURVES FOR CONSTANT HEATING RATE AND UNIFORM INITIAL TEMPERATURE The Semi-Infinite Solid -115

00.1 t 0 0.~~~~~02 0 q0..637. 0 0.3~~~~~~~~~~040.. -0. I U6 u5i -0.2 Un I I -0.3 U3 -1CSI~ -0.4 --- 8.5~~~~~~~~~~~~~~~~~~~U Figure El. Temperature Vs* Time gn=-O,o5- r=6

0. 1 0.2 0.3 I0.4 0.0 0.6 07. - ~ ~ ~ ~ ~ Fgue2. 0.2. -0.4 Un -— ~~~~~~~~~~~~~~~~~~~~~~~ —--- - 0.6 0.8 U n -t.0 j Figure E2. Temperature vs. Time gn — -1.0; r= 6

0.1 0.2 t0.3O. 0.5 0.6. 0 0.30. 0 - 0.4 m U6 - - 0.8 ----— c U5 Un - - 1.2 -1.6 - Zo0 Figure E3* Temperature Vs* Time. gn -2,oO; r =

t 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 0 - 0.6 _____ ____ "b -1.2 Un 3 Wm" —-Oo Now-OM I t -1.8 H9 u2 -2.4 -3.0 Figure E4., Temperature vs. Time gn= -3.0; r=6

d C(t) dt 1.5 0.6 0.5 0.2 1040IIiV0 dt o o 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 t Figure E5o Boundary Velocity and Position vs. Time gn -0o 5; r - 6

de (t) C dt 1.5 0.6 1.0 0.4 r df dt 0.5 0.2 O 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 t Figure E6. Boundary Velocity and Position vs. Time gn —l,0; r =6

dP dt 3.0 0.6 2.5 0.5 2.0 0.4 1.5 0.3 1.0 0.2 d t dt 0.5 0.1 0 0 _ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 t Figure E7. boundary Position and Velocity vs. Time. gn= -2,C; r=6

T Q(t) d t 3.0 0.6 2.5 0.5 2.0 0.4 (t) 1.5 0.3 1.0 0.2 dE dt 0.5 0.1 0 0 0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 Figure E8. Boundary Velocity and Position vs. Time. gn= —-3.'; r=6

APPENDIX F SOLUTIONS FOR HEATING RATES OF THE FORM H = h(1-cos2At) Semi -Infini te Slab -124-_

1.0 - 0.8 0.6 H (t) 0.4 0.2; I i i I I 0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 t Figure Fl. The Heating Rate Pulse. H(t) = l-cos2Tt

t 0 0.4 0.8 i.2 1.6 2.0 2.4 2.8 3.2 0 -0.I -0.2 Un - 0.3 H - 0. 4 _ _ _ _ _ __ _ FT Ul Figure F2. Temperature vs. Time H =C.1 (1-cos2-fft); 9, z-C-5; r= 6

t 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 0 -0.1 uS.0.2 Un - 0.4 U' -0.4 - 0.5 Figure F3. Temperature vs. Time H=l.C(1-cos27Tt); gn=-C.5; r=6

d f [M(t) dt 1.5 0.6 df 1.0 0.4 ____ __ Go 0.5 0.2 0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 Figure F4. Boundary Velocity and Rosition vs. Time H =1.C(1-cos2 7t); gn =-0.5; r =6

t 0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 0 u6 -0.1 0.2 Un -0.3 U3 -0.4 u2 05.5 Figure F5S. Temperature vs. Time H=1.5(1-cos27Tt); gn=-C.5; r- 6

dl d t 2.0 0.8 d f d t 1.5 0.6 1.0 0.4 0.5 0.2 0 0 0 0.4 0.8 1.2 1.6 2.0 2.4 2.8 3.2 I Figure F6. Boundary Velocity and Position vs. Time HI l1.5 1-cos2Tr t); gn= -C,*5; r;

APPENDIX G SOLUTION CURVES FOR THE FINITE SLAB -131

t 0 0.0625 0.1250 0.1875 0.2500 0.3125 0.3750 0.4375 0.5000 0 - 0.2.... I -2 -0.2 u6 - 0.4 Un -0.6 1u -0.8 -1.0 Fiure Gla. Temperature vs. Time 4C)C-.5; gn=-1.C; II=-2.C; r=6 Linear Station Distribution in x.

0.4 0. 3 I (t) 0.2 0. I d I -d -2.5 d t dt 0 0.0625 0.1250 0. 1875 0.2500 0.3125 0.3750 0.4375 0.5006 t Figure Gib. Boundary Velocity and Position vs. Time AC) =C.5; g -1l.C; H=2.C; r=6 Linear 3Station Distribution in x.

t 0 0.0625 0.1250 0.1875 0.2500 0.3125 0.3750 0.4375 0.5000 0 -0.22 u4 - 0.4 Un - 0.6 - 0.8 - I.O ~i'IO Figure G2a. Temperature vs. Time Z6) =C-5; 9,= 1.C,; H=2.CQ; r=6 Parabolic Station Distribution in xA,

0.5 0.4 0.3 e (t) F0.2 0.1 50 d e 50 d t0 0.0625 0.1250 0.1875 0.2500 0.3125 0.3725 0.4750 0.5000 Figure G2b. Boundary Velocity and Position vs. Time 2 (C) =0.5; gn= -1.0; H=2.C; r= 6 Parabolic Station Distribution in x.

t 0 0.1250 0.2500 0.3750 0.500 0.625 0.750 0.875 1.000 0 -0.2 U5~~ -0.4 Un -0. 6 - 0. 8 — I,0 Figure G3a. Temperature vs. Time. i~c) -i. c g" = -iC; 11 = 2 r: A Linear Station Distribution in x,

0.8 0.6.500 0.625 0.750 0.875.000 t[ ~ ~1 (C)=.C; g=-1.C; 2; r 6 0.4 0.2 dl t0 0.125 0.250 0.375 0.500 0.625 0.750 0.875 I.000 Figure G3b. Bonmdary Velocitytand Position vs. Time Linear Station Distribution in x.

t 0 0_125 0.250 0.375 0.500 0.625 0.750 0 875 0 -0.2~~~~~~~~~~~~~~~~~~~~~~~~~WS u5 0. 2 A ~~~~U3 -0.4 Un U0, Oor ~ ~ ~ U -0. 8 Figure G4a. Temperature vs. Time L (C)=.o;; g =-1 CP H =2.C; r 6, Parabolic Sta'ion Distribution in x.

0.8 0.6 (It) 0.4 0.2 -2.5 dt 0 0.125 0.250 0.375 0. 500 0. 625 0.750 0.875 1.000 Figure GLb. Boundary Velocity and Position vs. Time. Parabolic Station Distribution in x.

t 0 0.25 0.50 0.75 1.00 1.27 1.50 1.75 2.00 0 U6 - 0.2 u u - ~~~~~~~~~~~~0.8 6 Figure G5a. Temperature vs. Time Li(n) =a2.Ci = -l.C; tIb 2.Ci r 6 Linear Statig; Distribution in x

2,0 1.6 1.2 e~~~~~~~~~~~ (I) t } < l~ 0.8 0.4 -d 2.0 d t l0 0.25 0.50 0.75 1.00 1. 25 1.50 1.75 2.00 t Figure G5b. Boundary Velocity and Position vs. Time. (c) -=2. C; g -1.C; HI = 2 C; r6 Linear Station Distribution in x.

0 0.246 0.492 0.738 0.984 1.230 1.476 1.722 1.968 0 -0.2 - 0.4 U3 Un Ui -0.6 _______ -0. 8 I 0M Figure G6a. Temperature vs. Time ~(C>= 2.C; -g.Cc; H 2.C;. r 6. A Parabolic Distribution rf Stations in x.

2.0 1.6 1.2 e(t) 0.8 | 0.4 -2.0 dt o0 0.246 0.492 0.738 0.984 1.230 1.476 1.722 1.968 t Figure G6b. Baundary Velocity and Position vs. Time. (C) =.C; g=-1i.C; H=2.C; r= 6. Parabolic Distribution of Stations in x.

t 0 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000 0 "6 U5 -0.2 - 0.4 Un -0.6 -0.8 Figure G7a. Temperature vs. Time. -1.01 0.5; g =-1.0; H=l-coi~fft; r-6 Linear Station Distribution in x*

H(t)!(t) -1.0 05 -0.8 0.4 -0.6 0.3 H (t) -0.4 0.2 -0.2 0.1 0 0 -2.5 dt- 0 0.125 0.250 0.375 0. 500 0.625 0.750 0.875 1.000 t Figure GTh. Boundary Velocity and Position and Heat Pulse vs. Time. L(O)=C.5; gn=-1.c; H=l-cos4lTt; r=6 Linear Station Distribution in i.'

t 0 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000 0 - - 0.2 - 0.4 Un - 0.6 - 0.8 Figure G8a. Temperature va. Time. ~(o),o0.5; gn=-1.C; H=l-cos41Tt. r-6 Parabolic Station Distribution in

H(t) 1(t) -1.0 o 0. H(t) - 0.4 0.2 - 0.2 0.1 0 0 i -2.5 d t 0.125 0.250 0.375 0.500 0.625 0.750 0875 1.000 Figure G~b. Bpundary Velocity and Position and Heat Pulse vs. Time. J (C)C.5; gn -1.C; HIf=l-cos41Trt, r:=6. Parabolic Station Distribution in x.

APPENDIX H THE EXACT VALUE OF tm1 FINITE SLAB An equation for the value of the time required to melt the slab, tm, is derived by Landau in Reference 27 by means of a contour integration about the (xt) domain of Figure 1.5, Chapter IV. The same equation can be obtained, again for a constant heating rate, H, and a uniform initial temperature, g, in the interval 0 < x < e(0) by considering the energy required for warming the slab and for melting it. Under these conditions and with unit specific heat, density, and latent heat of fusion as in the nondimensional problem with temperature independent thermal characteristics the pertinent energy quantities per unit surface are arez -tmH = total energy supplied during the interval 0 < t < tmn -gie(O) = the energy required to raise the entire slab to zero temperature without melting;.(0) = the energy required to melt the slab without any change in temperature. The total energy supplied is the sum of the last two quantities. Thus -t H -- ('o) + (0o). This equation yields the following expression for tm', =F-~ 2,(0). -148

BIBLIOGRAPHY 1 Brillouin, M. "Sur queques problems non resolus de la physique mathematique classique propogation de la fusi.on." Annais de l'In.stitute Henri Poincare, 1_, (1931(, 285-307. 2. Carslaw, H. S., and Jaeger, J. C. Conduction of Heat in Solidsn Oxford Press, 1948. 3. Citron, Stephen J. Heat Conduction. in a Melting Slabo Institute of Aeronautical Sciences Report No, 59-68, 1959, 4. Crank, J. The Mathematics of Diffusion- Oxford' The Clarendon Press, 1956o 5. Crank, J. "Two Methods for the Numerical Solution of Moving-Boundary Problems in D;iffusion and Heat FlowT " Quarterly Journal of Mechanics and Applied Mathematics, 10, (1957), 220-231. 6. Datzeff, Ao "Sur le problems lineaire de Stefan II'. Annuaire de L'University de Sofia (Godishnik) Faculte des Sciences, Livre 1 Mathematiques et Physique, 46, Sofia, (1950), 271-335. 7. Douglas, J., Jro "A Uniquenress Theorem for the Solution of a Stefan. Problem." Proceedings of the American Mathematical Society, 8, (1.957), 402-408. 8. Douglas, J, Jro, and Galli.e, Ji. M., Jr,' "On the Numerical Integration of a Parabolic Differential Equation. Subject to a M4ovi.ng Boun.dary Condition'" Duke Mat+hematical Journal. 22, (1955, 557-571 9. Durand, WO F. Aerodynamic Theory,. Durand Reprin.t.;i.ng Comml.ttee, 1, Chapter 4, California Instit-ute of Technology, Pasadena, California. 1943. 10. Ehrlich, L. WO "A Numerical Method of Solving a Heat Flow Problem w'ith Moving Boundaryo" Journal of the Association. of Computing Machinery, 5, No. 2, (April, 1958), 161-176. 11. Evans, Go W., II. "A Note on. the Exis-tence of a Solution -to a Problem of Stefan." Quarterly of Applied Mathematics, 9- (1951), 185-193. 12. Evans, G. W., II, Isaacsorn, E., and MacDcr.ald; J. K. L, "Stefan-Like Problems." Quarterly of Applied Mathemati.cs, 8, (1950), 312-319.0 -1.49

-15013. Eyres, N. R., Hartree, D. R, Ingham5., Jacksorn, R., Sargan.t, R. J., and Wagstaff, J. Bo "The Calculation ctf Vartable Heat Flow in Solids." Philosophical. Transaction.s of the Royal Society of London, Series A, Mat-hematical ar.d Physical Sciences, 240) London, (1948)) 1-58. 14. Franklin., J. N. "N-umerical Stlability in- Digital and Analog Computation for Diffusion Problems." Journ.a:l of' Mathema.ics arid Physics, 37, (1959), 305-315o 15. Friedman, A. "Free Boundary Problems for Parabolic Equations Io Melting of Solids." Journal of Mathematics and Mechanics) 8) No. 4, (1959). 16. Goodman, T. R. The Ablation o f Melting Bodies with. Heat Penetration into the Solid, OSRTM 58-789,, August, 1958 17. Goodman. T. R. and Shea9 J. J. The Melting of Finite Slabs. AFOSRTN 58-824, August, 1958. 18. Hildebrand, F. B. Methods of Applied Methematics. New Jersey-: PrenticeHall, Inc., Englewood Cliffs, 1956,. 19. Howe, Ro M. Application of Difference Techniques of Heat Flow Problems Using the Electronic Different l.l Analy zer, Report AIR-10, Engineering Research Instiitute, Univers i ty of Michigan., A.rn Arbor, Michigan, May) 1954. 20. Howe, R. M., and Harnemarn? V S,, "So;Lutionr of Partial Differential Equations by Difference Methodcs Ursi:ig the Electronic Different~ial. Analyzer " Proceed:igs of the Insti,1,.ute of' Radio Engineers9 41, (October, 1953).)9 1497 21. Huuntley, H,. E Dimensior- a.i sis. New York: Rinehart And Co,, 1951. 22. Johnson, C. L. Analog Comput;er Techniq-ues.e New York" McGraw-Hill Book Co., O 1956 23. Kaplan, W. Ordlnary Differer*.tialI Equatio.ns0., Massachusett s: AdisonWesley Publishing Co,, 1958,, 24. Kolodner, i. I.o "Free Boundary Problem for the Heat; Equat-ion with Applications to Problems of Change of Phase " Comnmunlications on Pure and Applied Mat.hematics, 9, (1956), 1-31. 25. Korn, G. A. and Korn, T. M,. Electronic Analog Comput;ers. New York. McGraw-Hill Book Co., 1952.

-15126. Kyner, W. T. "An Existence and Uniqueness Theorem for a-Nonlinear Stefan Problem." Notices of the American Mathematical Society, Abstract No. 546-40, 5, No. 2, September, 1958, 27. Landau, H. Go "Heat Conduction in a Melting Solid." Quarterly of Applied Mathematics, 8, (1950), 81-94. 28. Longwell, P. A. "A Graphical Method for Solution of Freezing Problems." American Institute of Chemical Engineers, 4, No. 1, (March, 1958), 53-57. 29. Martini, W. R. Heat Transfer Aspects of Flame Propagation. Unpublished Report at the University of Michigan, 1955. 30. Miranker, W. L. "A Free Boundary Value Problem for the Heat Equation." Quarterly of Applied Mathematics, 16, No. 2, (July, 1958), 121-130o 31. Miranker, W. L., and Frish, H. L. "Analysis of the Nonlinear Stefan Problem." Notices of the American Mathematical Society, 5, No. 4, Abstract No. 548-96, August, 1958. 32. Murray, W. D. Numerical and Machine Solutions of Problems in Transient Heat Conduction. Ph. D. Thesis, New York University, April, 1958. 33. Murray, W. D. and Landis, F.'Numerical and Machine Solutions of Transient Heat-Conduction Problems Involving Melting or Freezingo I. Method of Analysis and Sample Solutions)" Transactions of the American Society of Mechanical Engineers, Series C-HT, (May, 1959), 106-112. 34. Otis, D. R. "Solving the Melting Problem Using the Electric Analogy to Heat Conduction." Heat Transfer and Fluid Mechanics Institute, (1956), 29. 35. Roberts, L. "On the Melting of a Semi-Infinite Body of Ice Placed in a Hot Stream of Air." Journal of Fluid Mechanics, 4, Part 5, (September, 1958), 505-528. 36. Rubinstein, L. "On the Solution of Stefan's Problem." Izvestia Akadamie Nauk, SSSR, 2, (1947), 27-53. 37. Scarborough, J. B. Numerical Mathematical Analysis, 3rd Edition, Maryland: The Johns Hopkins Press, (1955), Chapter 7. 38. Stefan, J. "Uber die Theorie der Eisbildung, inbesonders uber die Eisbildung im Polarmere." Annalen der Physik und Chemie, 43, (1891), 269. 39. Sunderland, J. E. The Transient Temperature Distribution in SemiInfinite Solids with Melting. Ph. D. Thesis, Purdue University, 1958.

901 0:3466.2059 -15240. Sutton, G. W. "The Hydrodynamics and Heat Conduction of a Melting Surface." Journal of the Aeronautical Sciences, 25, No, 1, (January, 1958), 29-3-2. 41. Aldrich, H. P. and Paynter, Ho M. Analytical. Studies of Freezing and Thawing of Soils-First Interim Report, Arctic Construction and Frost Effects Laboratory-New England Division, Boston, Massachuetts, June, 1953. 42. Soodak) H. Effects of Heat Transfer Between Gases and Solids. Ph. D. Thesis, Duke University, Durham, North Carolina, 1943o