APPROXIMATION TECHNIQUES FOR THE STEFAN PROBLEM by Hart Katz A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Department of Mathematics We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA 1974 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f the r e q u i r e m e n t s f o r an a d v a n c e d d e g r e e at t h e U n i v e r s i t y o f - B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by the Head o f my Department o r by h i s r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . Department o f The U n i v e r s i t y o f B r i t i s h C o l u m b i a V a n c o u v e r 8, Canada i ABSTRACT The macroscopic d e s c r i p t i o n of matter undergoing a phase change (the S te fan Problem) can be fo rmula ted as a se t o f coup l ed , non- l i nea r p a r t i a l d i f f e r e n t i a l equa t i ons . For the case of one space d imens ion , the t h e s i s deve lops th ree approx imat ion methods to s o l v e these equa t i ons . ) (a) Asymptot i c Expansions With the Green 's f u n c t i o n s to s u i t the g i ven boundary c o n d i t i o n s , the system can be t ransformed i n t o a se t of i n t e g r a l equa t i ons . For the case where the i n i t i a l phase grows w i thou t l i m i t as x 0 0 , the l a r g e T expansions f o r the i n t e g r a l s are c a l c u l a t e d and the l a r g e x behav iour of the i n t e rphase boundary found. (b) P e r t u r b a t i o n Expansions 1. When the l a t e n t heat of f u s i o n of a m a t e r i a l i s l a r g e r e l a t i v e to the heat content o f tha t m a t e r i a l , a u n i f o r m l y v a l i d p e r t u r b a t i o n expans ion i n a parameter r e l a t e d to the r a t i o of these heats i s p o s s i b l e . The f i r s t few terms of the expansions f o r the temperature d i s t r i b u t i o n and the p o s i t i o n of the i n t e rphase boundary are c a l c u l a t e d and found to be i n good agreement w i t h the .few:known exact s o l u t i o n s and n u m e r i c a l l y c a l c u l a t e d s o l u t i o n s . (c) Numer ica l Techniques Rather than use a t r a d i t i o n a l f i n i t e d i f f e r e n c e scheme, on l y t ime i s d i s c r e t i z e d and an a n a l y t i c e x p r e s s i o n f o r an . approximate temperature i i i s found f or each time step. This method gives good r e s u l t s f o r the temperatures over the time i n t e r v a l s required by the physics of the problem. F i n a l l y , the method i s applied to describe the freezing of a shallow lake. i i i ACKNOWLEDGEMENTS • My thanks go to the members of the Institute of Applied Mathematics; to Dr. G. Bluman for suggesting the problem, to Dr. L. Halabisky for helping me to see i t through, and to Dr. J. Varah for his patience with my complaints. The support from my wife during these Ph.D. studies i s a story i n i t s e l f . Quite simply, I could not have done i t without her. TABLE OF CONTENTS Introduction CHAPTER I Asymptotic Methods 1. The Semi-Infinite Bar 2. The Finite Bar CHAPTER II Perturbation Methods 1. The Case e 0 2. The Case £->-«> CHAPTER III Numerical Methods 1. Heat Flux Specified At x = 0 2. Temperature Specified At x = 0 3. Description of a Freezing Lake APPENDIX I APPENDIX II TABLES . . '" . BIBLIOGRAPHY V LIST OF TABLES TABLE I.a Comparison of perturbation solution with numerical sol-ution 1 1 4 TABLE I.b 1 1 6 TABLE II 1 1 8 TABLE III Comparison of numerical solution with exact sol-ution 1 1 9 TABLE IV.a 1 2 0 TABLE IV.b 1 2 1 I! v i LIST OF FIGURES FIGURE 1 Temperature dependence of c p . FIGURE 2 The interphase boundary FIGURE 3 Numerical scheme for the Stefan problem 82 FIGURE 4 Cross-section of a shallow., lake 94 FIGURE 5 Calculated and observed temperature profiles for Seneca Lake 122 - 1 - • INTRODUCTION There are a number of problems in heat conduction in which matter, in addition to supporting thermal diffusion, also undergoes a phase change. The problem is usually treated macroscopically; no attempt is made to detail the underlying molecular processes which govern the phase change. Further, complications due to mass motion (convection) are also ignored. Within these restrictions the system can then be characterized by a temperature f i e l d and the boundary between the phases. The evolution of the temperature f i e l d i s governed by the heat equation while the position of the interphase boundary is determined by an energy balance across this interface. The position of this free boundary is not known in advance but must be determined as part of the problem.. For these reasons the study has been known as the free boundary problem for the heat equation. The problem f i r s t arose over a century ago i n the work of K. Stefan [1] who studied the thickening of polar ice. In fact, the free boundary problem was originally known as, and is s t i l l often called, the Stefan problem. The methods he used were extended to describe other phase change problems such as the freezing of a lake, the recrystallization of metals, or the evaporation and condensation of water. The free boundary problem has been studied analytically in a few important special cases. It i s instructive to start with the general description of the problem and see what restrictions lead to these special cases. - 2 -In f u l l generality the system can be modelled as follows JR . R is considered to be a heat bath; i t is a source of heat to R but i s unaltered by R . (b) R contains a medium with temperature dependent con-ductivity k(4>) and enthalpy/unit volume H(<j>) . It is this temperature dependence which characterizes the presence of the two phases. In the absence of such "transitions both phases, are considered to have temperature independent properties. (c) The evolution of the temperature ty within R i s > determined by the partial differential equation (d) If 6) is the volume fraction of the i n i t i a l phase, then (a) Let R be a three dimensional region with boundary (at) where c^ i s the specific heat at constant pressure L i s the latent heat; assumed temperature independent. 0 - 3 -(e) There i s a temperature interval in which the medium undergoes most of i t s phase change. Let the midpoint of the interval be $ ^ , the length £ . Assume that Cp and k can be written (Figure 1) where Ag i s a phenomenological function such that for £ —» O where 0 is the Heaviside function. Figure 1: Temperature dependence of c - 4 -Here "f JT are the assumed temperature independent properties of phase I in the absence of phase II; analogously for £ k jE ) Since most of the phase change occurs in the interval •I $ M ~~ f ^ - f r ) i 1 : follows that 1* 1 ( 0 - 4 ) f 1^ The volume fraction can then be assumed to have the form The actual form i s not important; i t is required only that -*$tb-'$tty for 6-^0 so that (0.4) holds for e-> O It then follows from (0.2, 0.3, 0.5) that H can be written in the form where ^ i f * ^ ) i s continuous. (f) Finally, the boundary conditions that w i l l be considered are (i) Dirichlet conditions; temperature specified on 2 R ( i i ) Neumann conditions; heat flux specified on BR Consider now the t—O limit of the preceding system. At a fixed time t , let the surface SQ(t) be defined by (Figure 2) In the cases of interest, S q divides R into two regions: R y for '4* * ®H a n d f o r 4 " c $ H " E N C L O S E e b y a cylindrical p i l l box of linear dimensions e . The position of the p i l l box wil l vary in time; for v the velocity of y ^ w e have Figure 2. The interphase boundary - 6 -Now for the enthalpy H , where >/ = V-lr» for n the outward normal. Integrate (0.1) over V Apply the divergence theorem to the left-hand side and (0.6) to the right hand side to get H i s bounded, hence 1 U <JV- 0 ( t ? ) uniformly in t . The l e f t hand side i s O C £ Z ) Now l e t £—* O . The sides of the p i l l box do not contribute as there, —> O Only the discontinuities of k and H contribute at X 0 to give 7 -where the plus or minus sign i s taken i f v_ points into or respectively. The sign cannot be determined a p r i o r i ; that i s an additional input to the problem. The next restriction i s imposed by symmetry. Assume translation inyariance of the system i n two space dimensions; one space co-ordinate survives and the interphase boundary s(t) satisfies . <j?C£Cfc' , ~t}^ Q Further, focus on one of the phases. This can be done when the second phase is a heat'bath (its evolution would be known and i t s contribution to the f i r s t phase at the interface assumed given) or the second phase i s removed as soon as i t i s formed. . , The equations then reduce to 4 . where at . x = 5 6 f ) Here, H(t) can be considered a known heat input from the second phase. Under these restrictions there are two types of problems (a) The semi-infinite bar - 8 -At +•= O the f i r s t phase occupies the semi-infinite interval O^ X. <. Co . The partial differential equation i s satisfied on S 0 M < yc<. Co . An i n i t i a l condition 4 K«,o) = on (>£ ic c t>" i s needed to render the problem well posed. We consider the cases where — > d?t> for "X — > t» ; i t follows that 4 > ^ % » ' t ^ - J > 4>0 for X — * t*> . Non-dimensionalize as follows: fe where _& is a reference length and o~7~ • Then, letting — ^ — - *^9-M and - — - j — — ifly >^ C*} and suppressing bars, we get ~ 4 > T <r ^ x ^ Oo . . to.M Co 6 i ) where £s= ____ is the Stefan number. - 9 -For d e f i n i t e n e s s , assume and take the plus sign i n (0.8c). We r e s t r i c t the study to the cases where ~r~ "> o > these cases describe the melting of a slab. For c e r t a i n boundary conditions, closed solutions can be found a f t e r an appropriate change of v a r i a b l e s . With the transformation M _ , (0.8a) becomes CO where f— - ^ c f a , while the boundary conditions transform For l - ^ — 7 = a s o l u t i o n i s J x - — A — const; hence & — / 2 X T and ©c? A-«* 4 e~ "5 *1 The temperature d i s t r i b u t i o n s a t i s f i e s , i n the o r i g i n a l co-ordinates,, the i n i t i a l c o n d i t i o n 0 Jt >o - 1 0 -and X satisfies the transcendental equation This i s a Neumann solution [2] presented by Carslaw and Jaeger [3]. There i s a minimal H q required to support melting i.e. for which ( 0 . 9 ) has a real positive solution. We shall find this H q i n Chapter I. Another simple solution arises from the transformation v|_ <T » ' Then (0.8a) becomes while the boundary conditions transform to • city* For W l t ) = M 0 assume cn= | A 0 T . Then i f <j>„i.,)_» jg we get i 2- — M 11 -where U_=r ^ l O , 1 l ) This i s a solution f i r s t found by Stefan [ 1 ] . (b) The f i n i t e bar For $C«»1» O we look at the problem complimentary to the semi-infinite bar; rather than examine the diminishing semi-infinite phase, we treat the problem of a growing f i n i t e phase (which at \— Q i s not present). For SC**} > O the f i r s t phase occupies the interval •0£.v>c at t = 0 . Here s may increase or decrease. In both cases the part i a l differential equation i s satisfied on "Sc <• • Non-dimensionalize as follows where Si<^ i f S t * 1 > o , otherwise is a reference length. Letting _f »"4?*Tr* a n c* s u P P r e s s i n § bars, we get If <T(o)— o no in i t ia l condition need be specified. If <J"C"1 > o (hence &{ft\= 1 ) an in i t ia l condition is required: 4>(*.~w .o^**.* lowly For the condition at x = 0 either the temperature is specified, or the heat flux is specified. For definiteness we choose the minus sign in (0.12c). Let the problem (0.12a-0.12e) be called £>XC<*1«">, M ") and (0.12a-0.12d, 0.12f) be called ^ ( ^ 1 , H " ) . To get a closed solution, again let to give with the boundary conditions and i f « r (.<=>}=-- ^ For ^Ul=s 'cffw-J-t- , the problem *>) has a s o l u t i o n of Neumann type: 4Vl£**!>)« ^& 5 hence Here i s a parameter which i s the s o l u t i o n to the transcendental equation e f ^ ^ e . ( C M * . ) ) - 14 -The boundary conditions at X=«r(x) in both (0.8) and (0.12) are non-linear in <J~ ; hence the free boundary problem as a whole i s non-linear. As in most nOh-linear problems i t has been profitable, not to search for more closed solutions, but to develop effective approximation techniques. The principal results of this thesis are conclusions about the effectiveness of these approximation methods. (1) Asymptotic Methods In many practical cases only the behaviour of the free bound-ary <3*(?ty is sought. For example, as a lake i s freezing over, i t i s usually more important to know the thickness of the ice rather than the temperature distribution of that ice. In particular, knowledge of the i n i t i a l CT—^O) response of' <y to the boundary conditions may be sufficient. Such small time expansions have been calculated for <y in system (0.8) by Boley [4], Landau [5] and Evans et a l . [6]. In each case the authors transformed (0.8) into an integral equation for CT . The free boundary was then expanded in an asymptotic power series in J~X and the coefficients CF^t were found by matching l i k e powers of Jx "X- • Although i t has not been done, the same approach w i l l work for system (0.12). - 15 -There has been no d i r e c t attempt to f i n d the X * — 0 0 response of <T to the boundary c o n d i t i o n s . However, p r o v i d i n g cr—^oo as X—*r 00 , the work of Cannon and Denson-H i l l [7] can be used i n an i n ve r s e problem where the l a r g e x* behav iour o f a r i s s p e c i f i e d and the boundary c o n d i t i o n which g i ves r i s e to t h i s behaviour i s e s t ima ted . They show tha t f o r any C T C T Y > the temperature d i s t r i b u t i o n i n s a t i s f i e s 4**.— 4 V a n d t * i e f o l l o w i n g c o n d i t i o n s at x = <J" I t f o l l o w s tha t An asymptot i c expans ion f o r cr^ can be s u b s t i t u t e d and the behav iour o f *|y{«=>,-j") e s t ima ted . For example, i f <rv- | * X * f o r | * > « » j «»•> JL , an easy es t imate g i ves the bounds I t i s shown i n Chapter I t h a t , p r o v i d i n g tr~^oo as X — > « o , the l e a d i n g behav io r o f tr depends on l y on the boundary c o n d i t i o n s as - 16 -"X co and not the i n i t i a l c o n d i t i o n . Hence, es t imates l i k e (0.15) do not depend on the p a r t i c u l a r i n i t i a l c o n d i t i o n which Ay s a t i s f i e s . Gannon and Densbn-H i l l used the above es t imates to s tudy the monotone dependence of the f r ee boundary on the boundary c o n d i t i o n s , and not p r i m a r i l y to o b t a i n an asymptot i c expans ion . Thus t h e i r method- i s o f l i m i t e d va lue i n c a l c u l a t i n g h ighe r order terms i n the expans ion f o r CJ" • In p a r t i c u l a r , the e f f e c t of an a r b i t r a r y i n i t i a l c o n d i t i o n cannot be c a l c u l a t e d . Chapter I develops a method f o r c a l c u l a t i n g to a l l o rde rs the l a r g e "T behav iour o f cr f o r both systems (0.8) and ( 0 . 12 ) . Here the i n i t i a l c o n d i t i o n s present no d i f f i c u l t y i n the c a l c u l a t i o n . (a) The s e m i - i n f i n i t e bar For the system (0.8) we f o l l o w Boley [4] and extend the space domain to o& tc <. <c» a t the cos t o f i n t r o d u c i n g another unknown 4** l<\-c1» F l T l • The equat ions can then be so l ved i n c l o s e d form on t h i s domain. When the boundary c o n d i t i o n s at x = cr a re imposed on the s o l u t i o n a p a i r of coupled i n t e g r a l equat ions f o r CT and F r e s u l t . The a n a l y s i s breaks up i n t o th ree separate cases a c co rd ing to A c r * whether f\— approaches z e r o , goes l i k e a cons t an t , or i n c r eases w i thout bound as T~ y r p • F i r s t , the asymptot ic expansion of the i n t e g r a l — A M - ^ for T—^r 0 0 must be calculated. - j For A - ^ O C ^ ) t^ i e starting point i s the work of Handlesman and Lew [8] who use the Parseval theorem for Mellin transforms to derive an asymptotic expansion of . , . ij T C * * } C du ° for A—^ O . I n the case their method raises the second problem of calculating the analytic continuation in z of the asymptotic expansion of , for T".—*r ao . These expansions are carried out in the Appendix to Chapter I. For A > O ( l ) the integral (0.16) can be expanded by standard steepest descent techniques. It i s then straightforward to apply these expansion techniques to the integral equations of Boley. The main results are the following. [1] To support melting ( > o ) the heat flux H cannot Ax decrease faster than pr -t-*» ("p.) , X~^r 0 0 where [2] When H has the form jffff * "^j^ , x 0 0 t h e n A - ^ o - I n particular, when H - fe 4- ^ + * ' M A > ^ > 1 then <r = O C X 1 " 5 ) for S-c i and O0»»x") for. S - 1 • The heat flux contribution dominates the i n i t i a l condition contribution to the order calculated. [3] When H is of the form M= j=f 4- "^j^} , T , - > « » where H,> then A= O^ ") • I n t h e c a s e where H =. ~!= + * the i n i t i a l condition dominates at the second order in the expansion of the integral equations. It follows that <r, J T -f- -^^ -j.--- where A „ i s the solution to (0.9) and A p . can be expressed e x p l i c i t l y in terms of and j f d+ . . . [4] When the heat flux satisfies M > Ot^) , then A—?- 0 0 For the particular case H = O C T * * ) > b >— - L i t follows that 0"=r O^ -l****) • I n general, when A is unbounded the i n i t i a l condition contribution i s exponentially small relative to the heat flux contribution. For this reason a perturbation procedure applied to the partial differential equation gives the same results as the analysis of the f u l l integral equations. - 19 -(b) The f i n i t e bar For the ana l y s i s of (0.12) we follow Friedman [9] and Sherman [10] to generate an i n t e g r o - d i f f e r e n t i a l equation f o r CT The techniques developed i n the study of the s e m i - i n f i n i t e bar are used to analyze the dependence of CT on the boundary condition at the fix e d [x~0] end. Here the major i n t e g r a l to be expanded i s , 5 t Nevertheless, the ana l y s i s s t i l l depends on the behaviour of A"=" as *XT • *^ oo • The problem %>T.W%°\ i s analysed i n the cases ^ 0 ( 1 " ) , the problem (ii**) when A > 0(^1 . The main r e s u l t s are the following. [1] With 4-to ,-c1= ~ $<ST""%V o ^ v ^ l follows that A--2^ - o and the free boundary cf has the leading behaviour • ^ j- • l - v . 12] I f 5 g W - $ § J > t h e n A = 0(41 and <?^ JT\Jv where X^ , s a t i s f i e s (0.14a) [3] I f _ 4>*C«vO«-Hs&;l ' A > - x then A*-**eo . The behaviour of a~ for t h i s f i n i t e bar case i s - 20 -compared with the behaviour for the semi-infinite bar undergoing the same heating, but at x = cr . The free boundary <j" grows slower in the f i n i t e bar case. (2) Perturbation Methods The f i r s t suggestion of a perturbation procedure appears in a paper by Landau [5] who studied in a particular case the £ — > o and co limits of the semi-inf inite bar. Sherman [11] did a detailed analytic study of the case but the results are not suited to numerical calculation. In Chapter II we study the conditions under which the problem (J- j-C* 5) admits an expansion of the form e f c T c w o - « S e^crv . ( a 17) Even though these expansions are in fact asymptotic expansions ordered with respect to the parameter £ , the methods used to derive them are distinct from those of the previous section. To make that distinction, the expansions (0.17) will.be called perturbation expansions. The main results are the following. [1] When x \ i~ Q ( i l f° r T " — ^ °° > a perturbation expansion for and Yz=.J^<s'*~ , carried out in the fixed boundary co-ordinates Cy=Jr"jT) w i l l be uniform on o ^ y ^ 1 , T > o In most physical problems two terms of the expansion give good results. In the case that 4>t.<'«x^ ^ s oscillatory the regular procedure produces secular terms. These terms can be summed, however, - 21 -to produce a uniform expansion. The specific case <d> (~, =. 4- S'm .ax-i s worked out. Finally, even i f £ i s not small, the perturbation expansions can serve as asymptotic expansions for T — - • s > « « [2] If ^ f o T j ^ O f r O for X — t h e regular perturbation expansion f a i l s to be uniformly valid. In fact, after the independent and dependent variables are rescaled, the f i r s t order system incorporates the leading growth for <5~Cx7 but i s as d i f f i c u l t ' to solve as the original system. [3] A regular perturbation expansion for *Tct> also exists. However, by a rescaling of the variables, the f i r s t order system can be shown to be of the same d i f f i c u l t y as the original free boundary problem. (3) Numerical Methods There have been many types of numerical calculations for the free-boundary problem. One approach i s to solve the general problem (0.1) where the phase change occurs over a f i n i t e temperature interval. Although the enthalpy and conductivity are then functions of temperature, the domain - 22 -f o r the p a r t i a l d i f f e r e n t i a l equa t ion i s f i x e d . F u r t h e r , the method can be formula ted i n any number o f d imens ions ; the two d imens iona l case has been c a l c u l a t e d by Hashemi and S l i e p c e v i c h [12 ] , the th ree d imens iona l case by Meyer [13 ] . For the one d imens iona l c ase , when on l y the boundary c T i s r e q u i r e d , the i n t e g r a l f o r m u l a t i o n o f the f r e e boundary problem can be used . Chuang and Szeke ley [14] d i s c r e t i z e d the t ime i n the i n t e g r a l equat ions o f Friedman and so l ved f o r <T" a t the d i s c r e t e t imes . F i n a l l y , i n one d imens ion , when both the temperature and the boundary a re to be c a l c u l a t e d , f i n i t e d i f f e r e n c e schemes are used. A c l a s s i c paper on t h i s method i s tha t o f Douglas and G a l l i e [15 ] . In Chapter I I I we i n t r oduce a s em i-ana l y t i c approach f o r the s o l u t i o n o f the f i n i t e ba r . Only the t ime i s d i s c r e t i z e d and the boundary <r approximated by a s t r a i g h t l i n e i n each t ime s t e p . In each such t ime i n t e r v a l the s o l u t i o n to system ( 0 . 1 2 ) , s a t i s f y i n g a l l the boundary c o n d i t i o n s except ( 0 . 1 2 c ) , i s c a l c u l a t e d . Th i s s o l u t i o n , t i l which s t a r t e d at the beg inn ing of the n t ime step i s then made to s a t i s f y the f l u x c o n d i t i o n (0.12c) a t the end of the n ^ s tep to s t c a l c u l a t e the s lope of <T f o r the (n+1) s t e p . Th i s method was f i r s t suggested by G. W. Bluman [16] who used the theory o f L i e Groups to show tha t such an a n a l y t i c i n t e r p o l a t i v e f u n c t i o n i s a s i m i l a r i t y s o l u t i o n to the i n ve r s e S te fan problem w i t h a l i n e a r f r e e boundary. We de r i v e the same s o l u t i o n by a s i m p l i f i e d method which has the advantage of more d i r e c t l y i n c l u d i n g the case of - 23 -inhomogeneous boundary conditions. The s o l u t i o n i t s e l f has two possible s e r i e s representations; which representation i s chosen depends upon the boundary conditions being considered. The cases where f i r s t the heat f l u x and then the temperature i s s p e c i f i e d at " x = 0 are solved. These cases demonstrate where each representation of the a n a l y t i c s o l u t i o n i s more e f f e c t i v e . F i n a l l y , the method i s applied to the problem of a shallow lake which i s f r e e z i n g . We take advantage of the d i s p a r i t y between the ' d i f f u s i v i t i e s of i c e and water to reduce the problem to an equivalent one phase problem with an e f f e c t i v e l a t e n t heat and e f f e c t i v e heat -source. The r e s u l t s of the numerical c a l c u l a t i o n are compared with the data for a shallow lake i n northern Michigan and are i n good agreement for both water temperature and i c e thickness. - 24 -CHAPTER I: ASYMPTOTIC METHODS The o b j e c t i v e of t h i s chapter i s to determine the l a r g e X" response o f the boundary c r to the i n i t i a l temperature d i s t r i b u t i o n and the s p e c i f i e d boundary temperatures o r f l u x e s . These boundary temperatures o r f l u x e s are to be s p e c i f i e d a t x = cr f o r the semi-i n f i n i t e bar and at x = 0 f o r the f i n i t e b a r . The a n a l y s i s i s r e s t r i c t e d to the cases f o r which CF_2».«*D as V « o . I t f o l l o w s tha t the boundary c o n d i t i o n s p rov ide the l e a d i n g c o n t r i b u t i o n to the growth o f <r w h i l e the i n i t i a l c o n d i t i o n c o n t r i b u t i o n en te rs o n l y at h ighe r o rde rs i n the l a r g e X " expans ion f o r cr • The cases of the s e m i - i n f i n i t e and f i n i t e bar a re s t u d i e d s e p a r a t e l y . We beg in w i t h the s e m i - i n f i n i t e b a r ; the equat ions are e a s i e r to ana lyse and g ive d i r e c t i o n to the d i s c u s s i o n o f the f i n i t e b a r . 1. THE SEMI-INFINITE BAR The system to study i s ( 0 . 8 ) . S p e c i f i c a l l y , choose the p l u s s i g n i n (0.8c) and take $ M > 0 • T h i s s i t u a t i o n d e s c r i b e s the m e l t i n g p r o c e s s . F o l l o w i n g Boley [ 4 ] , extend the space domain o f (0.8a) from (<5", cp) t o •i c't o c ) . T h i s i s done at the cos t o f i n t r o d u c i n g another unknown which Bo ley chooses to be the f l u x at x = 0 25 -* £ O . T 1 — F ix) Then solve instead the problem 4***.—. o£- x • • .n-\+) The solution to (1.1) can be written in terms of the functions tyo* •% both of which satisfy (1.1a) a n d 1*** o 3 J > I * « , T I=- o <fcl*v 3ac % This gives ± i r e . ± -*ct--^ i •» where Gr U , T ; J - . f W - p i - J are the even and odd Green's functions for the half-plane. Finally, from the Duhamel theorem [3] the solution to (1.1) can be written - 26 -Note that the conditions (0.8b) and (0.8c) have yet to be satisfied. It is the application of these conditions that determines the equations for O-fT) . Before applying these boundary conditions we point out that the bar may have had a heating history before melting commenced. It is this history that determines the temperature distribution at the onset of melting i .e . ( p ^ t x l of (1.1c). Temporarily, let X== 0 A. mark the beginning of this history and <J7 the temperature distribution during this in i t ia l phase. . Then, since there is no melting, the true space domain is £ c , c » ^ and system describes the physics exactly. In this case F(x) can be taken to be a known flux. If, for example, the bar starts at a uniform zero temperature and is subjected to a constant heat flux input HQ , then set FCT} = — ^"o , > ° a n d 4^ ix") '= O i n t o get from (1.2) Melting wil l then commence at a time T^j such that ^(".T^ra. $ ' which can be calculated [4] as - 27 -With the onset of melting reinstated as T==* o and F no longer known, -system (1.1) i s used with Now, using (1.2), apply the conditions (0.8b) and (0.8c) at x = «*" to get To expedite the asymptotic analysis, make the transformation i n the f i r s t integrals of (1.3a) and (1.3b) to get - 28 -where l F C < - , X , K l and A l r ^ — ^ ' 4 T This i s a set of i n t e g r a l equations f o r F and cr . Although only the large X expansion for <T~ i s required, i t w i l l be necessary to f i n d the large X expansion f o r F i n the process. The asymptotic expansions of F and <T f o r X - ^ « w a r e calculated as follows: (a) Assume an asymptotic s e r i e s f o r F and <T" P —- H P f c Ai , M (1.**) where tAtf and are asymptotic sequences f o r X — a n ^ (TT %. it- Is c 1 - $ M ( x f e * - 2 9 -(2^ C^ tfc - a r e c o e f f i c i e n t s yet to be determined. For the types of heat fluxes being considered, i t i s s u f f i c i e n t to assume f o r F the asymptotic expansion where o° and V ^ - o , C L ^ O . In pr a c t i c e i t i s convenient to leave the asymptotic expansion f o r cf unspecified f o r the moment. (b) Substitute the expansions (1.5a,b) into (1.4a,b) and generate large X" asymptotic expansions f o r the i n t e g r a l s . (c) When the asymptotic sequences -^ f^c-T and L ^ *«-"^ have been chosen c o r r e c t l y i t i s poss i b l e to match to determine the c o e f f i c i e n t s (T^ . i n p r a c t i c e , once the form (1.6) f o r F has been assumed, the corresponding asymptotic sequence f o r O " i s determined by the requirement that matching be po s s i b l e . Step (b) i s the most d i f f i c u l t . The behaviour of A as X"—>cw i s c r i t i c a l . The ana l y s i s d i v i d e s up, accordingly, into three separate cases: A^OCO A^OCO and A>OtiV as "T** ^ oo . In t h i s chapter, the Stefan number S i s a parameter that does not enter into the c a l c u l a t i o n s , so without l o s s of gen e r a l i t y , set £.— 1 i n C1.4b). - 30 -CASE 1: J W o for T " — ^ Oo To complete step (a) i t i s consistent to assume that V = 0 in (1.6). In step (b) the f i r s t integral to be expanded is IPC^X^A^ . This i s a complicated calculation and i s done in Appendix I. For reference, the result i s 4- U ^ , ^ F U ^ o A % The f i r s t sum i s called the asymptotic contribution and for T — ^ Q° > The second sum i s called the domain contribut ion and for . TT* ^ vo 4- U x L- v. fe 4- L- I w. i r ^ i * Here, the ^ for which M 4- i i s called , and Q is the sequence I . If Q is empty the second sum in (1.7b) is absent. The M{c are unknown constants which - 31 -can be determined in the matching process to come. The next integral to consider i s , from (1.4a), where a 6* For simplicity, consider the case where 4 ^ * 1 ^ |*>0 , X " "y ©° • This i s not an unusual condition; for example, i t is met by the i n i t i a l condition (1.2a) which arises when the melting state i s approached in a "natural" way as was described in the introduction to this chapter. The exponential in (1.9a) can be expanded to get i H t e n — f OC^) CVio) The assumption on the i n i t i a l condition i s not c r i t i c a l . For more general i n i t i a l conditions the method developed in step (a) of Appendix I can be applied to (1.9b) to generate the required expansion for T — • - 32 -F i n a l l y , u s i n g (1.10) i n (1.8) we have The l a s t i n t e g r a l to cons ide r i s from (1.4b) •» The methods j u s t desc r i bed y i e l d the expansion T < M T . M « . ^ J;4>-15^57.+ OH*) - evil) The r e s u l t s ( 1 . 7 ) , (1.11) and (1.12) can now be used to w r i t e the i n t e g r a l equat ions (1 .4a ,b ) i n an expanded fo rm. For (1.4a) we have .C*=C£.si>-K**i\ ti,i) j- D l f j } air) - - tuft f I** J E £ + - W 4 5 + Oli.) w h i l e f o r (1.4b) we have - 33 -IT £ where from (1.6) I - — ^ } t - ^ c o • Examine f i r s t the l e a d i n g behav iour o f ( 1 . 13 ) . The l e a d i n g o rder on the R.H.S. i s 0C |^) , on the L .H .S . the domain c o n t r i b u t i o n dominates. To match, s e t ST then from (1.7b) i t f o l l o w s tha t t>Ft$ ,"* , - ! ) — ^ ,4.") • S ince "IT > matching then-gives Next cons ide r the l e a d i n g o rder o f ( 1 . 14 ) . The asymptot i c c o n t r i b u t i o n dominates on the L .H .S . From (1.7a) and (1.15) t h i s c o n t r i b u t i o n i s A F-( . | ; ,T t «>)- ' s o tha t matching to f i r s t o rder g i ves o r , u s i n g (1.16) and A— 4T V - 34 -7 c - T W H VWi^ 4^ V > ^ / We have assumed that A-> o , hence 0 ^ ) . For 11 - jrrr 4~ *> V.J» ) t h i s c o n d i t i o n i s f u l f i l l e d . This behaviour i s a turning point: f o r M < : melting cannot be supported as then ~r- <. & , whereas for H> we have A ^ O t f ) . The minimal <%x ' — H described i n (0.9) i s thus JJ} % 0 ft The matching has been c a r r i e d up to ^ ^ J ^ ^ . The i n i t i a l c o n d i t i o n c o n t r i b u t i o n does not enter u n t i l and thus has not contributed to the leading order i n the expansion. It i s yet to be determined how Cf-^co . This behaviour can be found by matching the expansion at second order. On p h y s i c a l grounds, the heat f l u x must be the agent to drive <r to i n f i n i t y ; the i n i t i a l c o ndition c o n t r i b u t i o n cannot enter even at t h i s second order. This i s demonstrated i n two cases. ' (i ) 00 a l g e b r a i c a l l y Here the heat f l u x must be of the form 1 H I + — T + where H^ > 0 and -^jCici. • F has a s i m i l a r expansion F, - 35 -Assume also that h^Sjvi^ , ftl^ > a , A0*> o . Then Now from (1.7a,b) So for (1.13), omitting already matched terms The R.H.S. i s of higher order. Matching to lowest order gives and 5- -Note that for we have o<: > so from (1.17) <*"-jy <x» algebraically. Again from (1.7a,b) J X T V - 36 -So f o r ( 1 . 1 4 ) , n o t i n g tha t & 4£.x~\,\ = o and o m i t t i n g already-matched terms The f i r s t term on the L .H.S . i s o f h igher o r d e r . Us ing ( 1 . 1 8 ) , the other terms can be matched to get so tha t f i n a l l y , f o r (1.17) T M-SV4-Note tha t the i n i t i a l c o n d i t i o n does not enter i n t o t h i s e x p r e s s i o n , ( i i ) <r~vO-» l o g a r i t h m i c a l l y In t h i s case the. heat f l u x must decrease l i k e H - / i r e . 4 " , x- * F must have the same expans ion For (1.13) use - 37 -to get, omitting already matched terms The R.H.S. i s of higher order. Assume that A ^ A„ v — 7 + - - - 0.19* T and match coefficients to get the ratio In (1.14) use so that using (1.19) and omitting already matched terms X N/A* X* Ci.aoy _ x {*ZF* _ it! 1 4 . . .' - . • • The lowest order i s OCj^ ~Uvr^ ' Match coefficients and use (1.20) to get Therefore 2\kv— and f i n a l l y Again the i n i t i a l condition makes no contribution to this leading behaviour. - 38 -CASE I I A=OCO f o r X — V o o Th i s case a r i s e s f o r a heat f l u x o f the form where H Q > . Set A= A,,-*- Ag.1^ where A t f i s a p o s i t i v e constant and Ag,t<-1^ yO f o r T—=J To complete s tep (a) i t i s aga in c o n s i s t e n t to se t V = o i n ( 1 . 6 ) . For s tep (6) the f i r s t i n t e g r a l to be expanded i s I F (j'V.T, A") Th i s c a l c u l a t i o n i s done i n Appendix I and the r e s u l t i s I F C r , x , A ) - - C t l > t > P C n x , t » * 0 A " (A 1*7 where C F^ T -PCJH+I-*^.) X r C i w + i - , r J » i H ^ ^ r , A.,") • 4 - ? i> -«- £ W • 0 i 1 ^ Here Q i s as de f i ned i n (1.7b) o r Appendix I,XTCA.V*,-fc) i s the con f luen t hypogeometric f u n c t i o n o f the second k i n d , and * - r " t n t » . A _°" r- |g-The next i n t e g r a l s are the i n i t i a l va lue c o n t r i b u t i o n s i n ( 1 . 4 a , b ) . To f i r s t o r d e r , these i n t e g r a l s have the same expansions as g i ven by (1.11) and (1 .12 ) . - 39 -With the r e s u l t s ( 1 . 1 1 ) , (1.12) and (1.21) the i n t e g r a l equat ions (1 .4a ,b ) can now be w r i t t e n i n expanded fo rm. For t h i s c a l c u l a t i o n i t i s p o s s i b l e to choose a heat f l u x such tha t the i n i t i a l c o n d i t i o n s c o n t r i b u t e at the second o r d e r . Such a f l u x s a t i s f i e s A p ** Assume tha t Ag,'-'^=,+ and put J^==. J" . Then (1.4a) has the expansion -ft *• w h i l e f o r (1.4b) ^ r t i ) IT C l , 4, O + 3%. v 4 F i r s t , matching to O C ) g i ves F p t ^ C t i . A > - f „ t v** . t i n * ) f v T j r C i . l , A . )= " fe) e*V C t ^ ) Nex t , matching to and us ing Ckc = j . - 40 - -H o € A * P. V> \, A „ J o ( 1 2 ^ ) F q can now be e l i m i n a t e d between (1 .24a,b) to de r i v e a t r anscenden ta l equa t ion f o r A*, ^ Compare t h i s r e s u l t w i t h ( 0 . 9 ) . With the A o f (0.9) i t f o l l o w s tha t « X '•ef1 l, and tha t equat ion (0.9) can be r e w r i t t e n as However, w i t h the Kummer t r a n s f o r m a t i o n ^CA«Kl"l= 1 T*(i4a-*»> 1~*?> > t h i s equa t ion can be t ransformed i n t o ( 1 .26 ) . Thus, to f i r s t o r d e r , the asymptot i c behav iour o f A w i t h H — ' -pr agrees w i t h the exact , = ^ = . T h i s i s to be expected on p h y s i c a l grounds and i s a check on the v a l i d i t y o f the expans ion t echn ique . U s i n g , say ( 1 .24a ) , F q can be e l i m i n a t e d from (1 .25a,b) and the s o l u t i o n f o r found - 41 -where Note that i t was possible to eliminate M o Finally, the expansion for «r" is Here the in i t ia l conditions enter at second order. CASE III A-> o o for T _ j y o » This case arises for a heat flux satisfying For definiteness, choose a flux of the form H — V 4 0 x 1 ' T — > ° -* t»>— . Then to complete step (a) it wil l be necessary to choose V> o- in (1.6). To begin step (b), examine the role of the in i t ia l condition - 42 -i n t e g r a l s i n ( 1 . 4 a , b ) . In (1.4a) we cons ide r The f i r s t term a r i s i n g from G + i s L a p l a c e ' s method then g i ves the es t imate As i n Case I and I I cons ide r those 4^ f o r which 4^,0) ^ O C ^ ) I* y c? , X — - j . < X > • Then and s i n ce <T > A f o r X — i t f o l l o w s tha t T, 4^ i s e x p o n e n t i a l l y s m a l l f o r T~=y p « The second term a r i s i n g from G + i s 1 e -Here the boundedness o f 4^ , and the Lap lace method y i e l d the es t imate where - 43 -Together w i t h (1.27) t h i s es t imate i m p l i e s tha t X4^ ^ ^ ^ X A ^ ' S i m i l a r l y , i n (1.4'b) , f o r we have tha t Xcf^ ^ O Cy^ ~ ) ' The l e a d i n g term to be matched i n (1.4a) i s $ w ^ • The i n i t i a l c o n d i t i o n c o n t r i b u t i o n i s t he r e fo r e e x p o n e n t i a l l y s m a l l r e l a t i v e to t h i s l e a d i n g o rder and so does not p a r t i c i p a t e . i n any order o f the match ing . The same comments app ly to (1 .4b ) . Combine s teps (b) and (c) to ana l yze the f i n a l i n t e g r a l With V^ e> T ^ o r XT—^ O C T , assume tha t Ac,x**' • Then F c o n t r i b u t e s to the exponen t i a l i n the in tegrand of IF to get where ^CtO = * (z~tzY^~ A .^ " I f f has an i n t e r i o r maximum at then I F s C I C ' "XZ J so to match the exponen t i a l on the R.H.S. of (1 .4a,b) we r e q u i r e - 44 -" f t V ^ l s 3 • U s i n 8 -f'l**»fll=0 i t follows that U M = and -£ ( U H ) so that -f-C^tt) is indeed an interior maximum. The Laplace method can be used to expand (1.28). When the expansions are used in (1.4a.,b) and matching carried out to first order the leading behaviour of cr is found to be ^ ^ — T b * * \\,ud) For the case b = 0 the result reduces to (0.11). Again, the asymptotic behaviour of CT for H—» H 0 , oo corresponds to the exact solution for H = H o Higher order terms cart be calculated in a straightforward manner. However, a faster procedure to generate an asymptotic . expansion for <T is to use the partial differential equation (0.10) directly. Take, for example, the case \Xl where Vl^ lt : )—V o f ° r =^ ©° . Then rewrite (0.10), introducing the art i f ic ia l parameter \ , as follows J/r-where =- j j ^ . • Now expand in powers of A - 4 5 -* X Ar- XT* + I f Kfc,^ " Cr*6=) A N D — "•CT^) for T — , then with A 1 2 3 1 the s e r i e s (1.29) are asymptotic expansions for the s o l u t i o n s . The f i r s t order system has the governing equation For t h i s equation an a r b i t r a r y i n i t i a l c ondition cannot be s p e c i f i e d . But by the previous d i s c u s s i o n the i n i t i a l c o n d i t i o n makes only an exponentially small c o n t r i b u t i o n to the expansion for ^A. • Consider, for example, the case where U & - » o a l g e b r a i c a l l y . Then a simple . c a l c u l a t i o n gives f o r p. ( a f t e r s e t t i n g \=*-\ ) Ji=- + i i * . U „ $ H JW„ + exponentially small terms . This i s an asymptotic s e r i e s f o r X—^oo • to f i r s t order i t agrees with r e s u l t (1.28a) obtained through the rigorous expansion of the i n t e g r a l equations. - 46 -2. THE FINITE BAR For the f i n i t e bar, the large TT response of the boundary <jr(x") to a heat source at x = 0 i s examined. Here, as opposed to the d i r e c t heating case of the s e m i - i n f i n i t e bar, the heat source i s remote from the free boundary and must d i f f u s e through the material before a f f e c t i n g the phase t r a n s i t i o n process. As a r e s u l t , the boundary cr responds d i f f e r e n t l y to the same heating applied i n these two d i f f e r e n t s i t u a t i o n s . The relevant system to study i s (0.12). S p e c i f i c a l l y c h o o s e the minus sign i n (0.12c) and take H = 0 . For the analysis an i n t e g r a l equation f o r cr must f i r s t be formulated. Consider f i r s t the problem as defined i n the Introduction. Following Friedman [9 ] , suppose that <*"(*cj * form a s o l u t i o n to the problem. Then the Green's i d e n t i t y can be integrated over the domain o < X- <• *»"trl ,: o<. £<• T < C T - £ to get, upon l e t t i n g o Now set )£=<"'Cxi to get an i n t e g r o - d i f f e r e n t i a l equation f o r <rCxl : - 4 7 -For -^0,01"1 a similar analysis gives r1 + — of-As in the semi-inf i n i t e bar, the behaviour of A = ^ . j - for X"^ =?-oo> breaks the analysis up into three distinct cases. > Case I t \ < o l O For this case analyse the system fJjOf, «»1 where for "H^>o, V> c . As in the semi-inf inite bar, in order to carry out the matching, the integrals in (1.30) must f i r s t be expanded. The f i r s t integral to consider i s \ ^ Cr ( c r t x j . x ^ t x l . x l d T . In what follows l e t - 48 -Cx.x \ » 4 ( X - T 1 The f i r s t term of Xtf" » is straightforward; since — 0 ( A ) uniformly on o ^ X - ^ X * the exponential can be expanded to get > The m*"*1 term of this expansion i s ( O C A ! * * } • For the second term 4-the argument of the exponential has a singularity at X"_ -XT Consider the transformation then U ^ o for T - * o and for X — ^ x . If this transformation i s applied to ~Xxr the analysis in case I of Appendix I can be used. The result i s - 49 -0 0 In the nomenclature of Appendix I, the f i r s t sum i s the domain con t r i b u t i o n . Here T*<r Cx,"^ ) i s the a n a l y t i c continuation i n £ °f The second sum i s the asymptotic c o n t r i b u t i o n . Here, f o r X — . X (.X the inverse transformation of ( 1 . 3 3 ) , the function l S ^ U , - ) ) 1 has the expansion f o r r . f O S ' , The and l~ are yet to be determined. *»» X, h i H i The m = 0 term of the domain c o n t r i b u t i o n cancels the m = 0 term of (1 .32)- , so consider the next o r d e r — t h e m = 0 . term of the asymptotic c o n t r i b u t i o n . -For 5*"(Xw)assume the expansion o 0 & (X" w1 ^ X — T + . . . u—>• 0 0 ta-Then expand (1.33) about X— X and substitute in the above expression to get The leading behaviour-of i s then so that the asymptotic contribution to (1.34) has the leading behaviour o • • i 1 The results (1.32) and (1.34) can then be combined to give T - <J<T The second integral to expand is This i s in a form covered by Appendix I , and to f i r s t order the result i s - 51 -The last integral is Here, since the domain of integration is finite and A—^ o , the exponentials in can be expanded to get Now assume that cr,- J*X<*" , <*.-^ -0 » T-*y <» , and collect the results (1.35), (1.36) and (1.37) ~ ^ v u r n where For 0^y a , (1.35') dominates (1.37'). We expect this result. In the phase transition process the in i t ia l heat contained within the bar contributes to the latent heat of fusion to cause <y to grow. But then no finite in i t ia l amount of heat could provide the infinite amount of heat required to cause O-^oo . Matching the powers in (1.35') and (1.36*) implies that - 52 -Ct=- '•» s o °-y° f° r 4. . This i s also an intuitively clear result; i f ^ i s not integrable for X-^ -cw , then cy~--^- oo . Finally, matching the coefficients gives ga^ —^3^ so that V I . -CASE II A- Q(0 For the system this case arises when 1x1 ^ I X " ) where <J?£ is a constant and o for X ^ C ° . Then A—- A^,4-Ap, lTl where A j z . - ^ o for X—<v> • As has just been shown i n Case I, the i n i t i a l configuration does not affect the leading behaviour of <r for X—Y«» . In particular, the solution to B-xCljo) approaches the solution to ^>j-l"t °1 providing <£^lx) is the same for both problems. For the same reason, the solution to approaches the exact solution (0.14), so i t follows that & 2. Jh0T ^ T — o o where A,«> satisfies the transcendental equation J -e-a. C A S E I I I A > 'oiA ) For this last case examine 8j£ with a heat flux satisfying - 53 -where U-0 > o i s a constant and ' 8 4 ^ — ^ 0 for X—"froe» . The i n t e g r a l s i n (1.31) must be expanded f o r X — * * » Consider f i r s t the i n t e g r a l In the notation of Case I, f o r the f i r s t term I f<r we have Tf^*—OCA.) for c ^ v ^ T " • Define the transformation then j[ f o r T - ^ o and u - ^ o f o r X - ^ - X . Since ©° , only the u - ^ o behaviour of the transformation i s needed. A simple o c a l c u l a t i o n gives Therefore, the f i r s t term has the leading behaviour The second term i s exponentially small r e l a t i v e to the f i r s t so that Use the change of v a r i a b l e - ~ ~ X rj-«. t*-J t o transform t—x x the second i n t e g r a l 1 ^ 1 x 1 = - / M R m C^t<Hic\x-, * / c l <Tc t o ^ i v e - 54 -to get for t"J»*' the leading behaviour A It is easy to show that the last integral r l _ ' •"" I f l i^ l x W J <$?(>-) CV C<y£x),T', la-jo) J x - has the leading behaviour i — A r so that I ^ W = 0(4 ) l 8 4 & C T l Assume that IWp is of lower order than 14-^ • Then, matching (1.39) and (1.40) gives This is a transcendental expression in A and cannot be simplified further. However, i t follows that h— *> (lr»x") so that ~p ^ & A and X4> I s indeed of higher order. Finally, this result for the heating of the finite bar can be compared with the results for the same, but direct, heating of the semi-infinite bar. From case III of the semi-inf inite bar, c?~— Q (x/) » so that A= Olx1 as X__^-c»o f whereas for the remote heating - 55 -of the f i n i t e bar we have j u s t shown that A— o C l^t} SUMMARY The behaviour of crCcI for X—^<*» has been c a l c u l a t e d f o r a v a r i e t y of boundary conditions. Although, i n most cases, only the leading behaviour has been c a l c u l a t e d i t i s c l e a r that the determination of the higher orders i s straightforward. Further, the asymptotic forms f o r the boundary conditions were chosen f o r demonstration purposes only; the methods presented i n t h i s chapter do not depend upon the forms chosen. - 56 -CHAPTER II : PERTURBATION METHODS The free boundary problem encompasses two distinct physical mechanisms: (a) Heat transfer through diffusion. This mechanism is represented in the equation which is valid throughout the material (b) Phase transition. This i s represented by the flux condition at x = a r c - t ^ Here £.= j i s the non-dimensional Stefan number which specifies the coupling between these two physical modes. For the £—*f-o case, the problem &j R O» C ' ) i s examined. In this case i t is possible to calculate a perturbation expansion of the form C *. ___ tr* - • • The arguments of 1 have deliberately been l e f t unspecified. In fact, the original co-ordinates ( X , T ) are. not suitable to generate the expansions. Problems of uniform validity of the expansions for o ^ x ^ t o arise only in the particular case t X — ^ c for f__^.«w> . In this - 57 -s i t u a t i o n , however, uniform v a l i d i t y does hold i f and only i f i s bounded f o r X—> i » . When the approximation i s uniform there are the following a d d i t i o n a l features: (a) In ( 0 . 1 2 ) l e t K~ H^XW^tUJJ \ Then by r e s c a l i n g the temperature • • the system remains inv a r i a n t except that now l$g|<<£ 1 , I<4*,-!^ -4.. and i n ( 0 . 1 2 c ) there i s a new expansion parameter For many p h y s i c a l problems £,H«, i s small; f o r example, i f i c e i s being melted by water at 10°C we have e M 0 ~ „ 1-4.This s c a l i n g argument then ' guarantees that f or £ tta small, even i f the expansion i s c a r r i e d out i n the o r i g i n a l non-dimensional co-ordinates, only a few terms of the expansions ( 2 . 1 ) are needed to give good r e s u l t s . (b) The expansions ( 2 . 1 ) are asymptotic f o r X — > o o . Thus even i f c r i s i s not small, the perturbation method i s a d i r e c t way to generate a large X expansion. For the £_^e>o case the s e m i - i n f i n i t e bar i s examined. I t t i s shown that a uniform expansion i n powers of £ i s p o s s i b l e , but that the equations to determine the f i r s t order system are as d i f f i c u l t to solve as the o r i g i n a l problem. 1. THE CASE £ Consider the problem & c(' l ,«>} where <r_^-©o for T ; — V * 3 0 Set H = o in (0.12c). It i s easy to show that a perturbation expansion of the form (2.1) in the original co-ordinates w i l l not work. Indeed, from (0.12c) i t follows that — =• o so that o^ix^ cfl^)— \ • T n e boundary condition at x - <r can now be expanded: (a) From (^.^ "O^ c we have (b) From — t ^ t ^ - c ) — we get Note that, since <jr_^.©» , these are already non-uniform expansions. Consider, in particular, the case where The 0(1) systems i s , for T—yoo The i n i t i a l condition need not be considered as i t does not contribute to the leading behaviour of the solution as T-^©^ . Now for large -r; i t follows that This expansion can be used in (2.2) to give Therefore Cf%-= Ol"C) • But we have shown in Chapter I, case II of the f i n i t e bar that cr— O C / T ; ) • The non-uniformity of the expansion is then apparent. To avoid the non-uniformities presented by the expansions of the boundary conditions at x = cr , use instead the fixed boundary representation (0.13) to generate the perturbation expansion. . system CO.13) can be rewritten as together with the flux condition The analysis breaks up into three cases according to whether ^-Ofa'J « ^ O O V - o r ^ > G ^ , as T - * c ~ . - 60 -CASE I Q{0 Expand as follows and rather than cf , expand R To analyse the systems generated by substituting these expansions into (2.4) the following results are needed. • Let <j» satisfy the diffusion equation with non-homogeneous term £f«*,T) * * n ~ S i ts , with the boundary conditions Then ^ can be written Where 4^c arises solely from the i n i t i a l condition, 4 ^ from the boundary condition, and ^ from the source term. See Appendix I I . -61 -One can show that (j? L*it%}~ £HC ) • Since we examine only the large time behaviour of the solutions this term can be ignored. Further, assume that for T*** a n c * l*b> i ^ < ? ° The Bromwich contour integrals in Appendix II can then be expanded for large X to give *7T sr. n v ^ n T -t-CH* ) where C^I%*1== Specifically, to ensure that CTL-^-co and yet —>c» > choose ° ^ •yfef< ± • The f i r s t order problem i s then - 62 For XT—$»- oo the leading behaviour is then, from (2.5a) ley So from (2.4a) we have ~ X ° , hence This completes the first order calculation . The second order can now be calculated as * The 0(1) system can be written in the original co-ordinates as This is precisely the equation one would write down to describe diffusion in a bar with one end stretching at a rate . An element at x would then experience a velocity " I T C J I V S . ^ , giving rise to the convection „ . ^ accounts for the change of thermal conductivity of the bar due to the redistribution of matter caused by the stretching. - 63 -From (2.5b) the leading behaviour is then So from (2.4a) we have, providing -V T • * Clearly and T*=^(T~) for T — V 0 0 For k -^ o the general problem is A tr iv ia l induction on this system, using (2.5b) repeatedly, demonstrates that E . ^ = c ( P ) e ^ and T ^ <> (.T ) as T oo . In fact, i t is easy to show that T t = O C x - ' k v - ) "x, then (2.4a) can be integrated to give - 6 4 -Thus, even i f e is not small, but s t i l l 0 ( 1 ) , and the expansions break down as perturbation 4>~ 21 £ * - T C b > and 71 e* expansions, they are, for £ fixed, asymptotic expansions for <*> -Finally, note that R q is the largest term in the expansion for R . Thus from (2.3) and (2.7) the leading behaviour of of can be calculated as This agrees with the asymptotic result (1.38) with £=- \ To determine its accuracy, the perturbation solution can be compared with a numerical calculation. Since, when <r is .unbounded., asymptotic agreement as " I T — h a s already been shown, i t is sufficient to compare the solutions in a bounded case. For this case, moreover, the numerical scheme is most accurate and comparison most significant. The system (2.10) can be used to generate such a perturbation solution. Take for simplicity for J 7 T | % taTT, Y\= o, \ - - , . This satisfies the 0(1) equation and can be used to generate higher order terms. A simple calculation gives for the next order - 65 -. 1 To second order, therefore, the solution for the temperature i s This solution satisfies the 'boundary conditions and has the corresponding free boundary where A C i - eT^ % ) The numerical scheme, as outlined in Chapter III, can now be used with the given boundary conditions. The computed temperature and free boundary can then be compared with the perturbation results ( 2 . 1 1 ) and ( 2 . 1 2 ) respectively. - 66 -Choose "g = 1 . The calculation i s then carried out up to £ at which time the boundary <$~ has essentially stopped growing. For £ — only the O(c) correction i s needed to achieve two •figure accuracy for both the temperature and the free boundary. For two figure accuracy i s maintained when the Q contribution i s added. .See Tables l a , and Ib. CASE II J X t i ^ OCil For the case •$^lx\J}¥„+'%-C'+.-. Y± y o » and X—^- o& » the results of case I can be applied with V = o From (2.6) and (2.8) ? ! while from (2.7) and (2.9) ' For X—>©e» these expansions can be compared with the exact solution CO. 14) with <f?ft=. T-J'^ . F i r s t , using (0.14d), expand A f H i ) in powers of £. (215} - 67 -and note from CO.14b) that lT)=^ • Now, since R—> 0 0 , i t follows from (2.3) that f / - c P for X—*©o . Therefore, from (2.14) and (2.15) we have r i x l — t I a l T l , x — yoo ^ Next, substitute (2.15) into (0.14a) and expand 4^1^-.^ ^ C 4 - ^ ) + + Q U M % Comparing (2.13) and (2.16) we have It is clear that these results are also valid for solutions to £-yl.o,o) with the same boundary condition. Thus, for both cases, the time-independent contribution to (2.13) and (2.14) can be summed and the perturbation expansion written as where "V , W —> 0 for X — a n d <T(0) takes on the values 0 or 1 . Providing (j^txl i s bounded and approaches a constant as -r; ^ co > we have just shown that the perturbation solution has the 68 correct asymptotic behaviour. However, i t can happen that C x i is bounded but that $>fglx\ does not exist. For these oscillatory boundary conditions the regular perturbation method breaks down. Consider, for example, the simplest case f ^ l T l - 1 ^ 4 . \LS\*JIT > ( 2 , 1 0 ) Here, take ^ 5 < 3 > 0 and 8 . If the expansion procedure as outlined in case I were used, the terms 2.2-4^ and *j <£\^ in (2.4) would generate secular terms of the form T ^ i t i A t , <io£JTT tv -^l *l Therefore, the perturbation expansion would not be uniformly valid for XT—> 00. ' ~ To remedy the problem, the expansions (2.17) for r and 4* can be used. Here, 4*» i 9 ) a n d AC*^) are assumed to be derived from the time independent part of IPoix'S *-n (2.18); The important feature is that even though and \ depend on £ , they are given by (0.14) and thus can be calculated prior to the perturbation expansion. Put A 0 - A T n e n f o r either problem & TGt,«0 o r % i& t>"\ > substitute (2.17) in (0.13) to give with the boundary conditions - 69 -and, i f c F - i c I — J . >. the i n i t i a l c ondition F i n a l l y , the f l u x condition becomes simply Now proceed as i n case I and assume the expansions E L ^ W f c l T - l „ We show that the secular terms are no longer present. In p a r t i c u l a r , we derive a large X * expansion f o r WQ and show that i t i s bounded fo r X—-> oo The 0(1) system i s The behaviours of f^-C^, o* and %j-\oto\ f o r X — > c*» are - 70 -the same but the case <jr£<e>l= O is simpler to analyse. For this reason,' choose o t v l s o in (2.20). Then no i n i t i a l condition need be specified and the solution can be found first by taking the 'Mellin transform of (2*20): fco where The solution is where, with M(a,b,-Z) the confluent hypergeometric function, To get the solution to (2.20), invert (2.21) where A is a Bromwich contour in the strip —A^ <. . The left limit of the strip is dictated by the right most pole of f*C$ft $'tt\J£ j , • The right limit arises from the convergence requirements at y = 0 as follows. For §>= £ > c , using S t i r l ing ' s approximation for' ^ * C | > 1 we have for f X e m p j — > • Because of the osci l lat ions of X ^ along A- , convergence of (2.21) for y = 0 then holds for ^ Providing the integral (2.22) continues to converge, the contour can be displaced to the right to derive a large X expansion f o r 1^" . To -show this convergence, start with the estimates [17] •ft* _ l t % / I H ^ 6 " C2-24-) /if so that Hence with = arg p , for !Xfr»»j,|—along any Bromwich contour. F i n a l l y , using"(2.23), (2.25) and that t^ ^^ .|->^ for f 3Vnf>8-*<*' , the integrand of (2.22) damps l ike -- 72 -for (I»«|> j - sy-eo a l o n g J > = . For &j > O this estimate allows the contour fc. to be shifted to the right and the contribution at l l c^p\— oo to be ignored. From (2.19b), <M%^ i^t^ , s o that and from (2.22) we have, therefore Displace the contour £ to the right to pick up poles at p = 1 , and on the set The set i s a positive, s t r i c t l y increasing divergent sequence. Further, for a l l n, O C * 5 T „ } [17]. Since X 0 ~ O Cf) even c ^ , , i s large, and the asymptotic estimates (2.24) can be used. Hence we have -J _ »• Tf ' and 2 * ' - 73 -where for a Bromwich contour, in the strip <1 f ^ | » < *^H+J. ' F R O M the previous estimates it is easy to see that O C T " ) and so (2.26) is an asymptotic expansion for T—y co . i t is not , convergent. as (xl l -»«» for X fixed, N—>co . Note that for. T - V > ° , WQ is bounded and is 0(E.) uniformly. Thus, in the expansion for r r= X^T-+ eW^4-2 the term C^g is 0 ( £ ) and is uniformly smaller than the first order term. The expansion for <y follows simply from (2.3). The expansion (2.26) for WQ holds as well for the problem g Ctj* 7! * t n * s case, however, it is worthwhile to compare the perturbation solution with a numerical calculation to determine the accuracy for X. moderate. - 74 -Put CTlc4= in (2.20); then an i n i t i a l condition i s required. For simplicity, choose So that from (2.19a), ^ ^ > e J » o . The problem (2.20) can then be solved by an integral transform in y [18]. The solution i s where for ^ H ^ C+M*. f ( 4 4 - v v . J l ^ 7-F°r A«, small and <;< <sto , the following approximations are valid: 75 -and °^ fr,'~ ~\~ ' These approximations can be used in (2,28) to give 5. t o r ^Tr<4) This representation is valid away from y = 0 and can be summed to Hence, from (2.19b), -g^ —- |* —^—- , so that Finally, from (2.17) The solution for the free boundary cr of f>^£4,©) with boundary condition (2.18) and i n i t i a l condition (2.27) can now be calculated numerically and compared with (2.24). Rather than f i x £ and solve the transcendental equation (0.14b) for A- , choose - 76 -\ o and calculate the corresponding £ . For A^=--i we get £=.-W34 • Then with the correction term WQ in (2.2a) two figure accuracy i s achieved over the range O ^ T - ^ -io • See .Table I I . CASE III > a ^ For this case, the regular perturbation procedure w i l l not yield a uniformly valid series. Attempts at multi-scale expansions w i l l also f a i l . A set of scale transformations can verify these statements. Assume, for def initeness, that ^ , ^ 1 ^ 0^X^ > sf y ^ T — ^ t * * Make the following transformations: £ 4? Hence T - a . £ , x and ~~7Z-=• & • Finally, put Then <*r — £x and *~ yx~ *~ ^ . The system (0.13) becomes ^ « * ^ ± and - ? M U , l W * Now we can write <J^ f-s.1=r t 4- 3£p. l where for some k > 0 > f ^ lx\ £ 0CV*"*) . It follows that $ Assume a regular expansion of the form I T , 12=0 The f i r s t order system i s therefore (dropping the bars) and — T Thus the f i r s t order system incorporates the leading growth behaviour. Any multiscale expansion procedure would have t h i s system as i t s f i r s t order problem. But t h i s system i s as d i f f i c u l t to solve as the o r i g i n a l one; hence a perturbation expansion i s not us e f u l here. 2 . THE CASE g_y ©o Consider the s e m i - i n f i n i t e bar. F i r s t note i n the s p e c i a l - 78 -c"<r case (0.10) and (0.11) that admits a regular perturbation expansion m powers of £. Hence, we do not expect the expansion in the general case to be of a singular nature. Let &-*r £• in (0.10), to give the system Examine the melting case where ^ ><? • T n e n > t o support melting, Het) > o Now assume a regular expansion for 4* a n d <y • Then the 0(1) system is The flux condition at y = 0 no longer contains < J * e x p l i c i t l y . The dependence can be made explicit, however, by using (2.31) at y = 0 to get - 79 -' O This system i s no simpler than the original one, hence again a perturbation expansion i s not useful. Similar comments hold for the f i n i t e bar. SUMMARY • For the one-dimensional free boundary problem we have shown, for a selected set of boundary conditions, when a feasible perturbation expansion i s possible. For other situations, eg. H ^ 0 in (0.12c), , or the flux rather than the temperature specified at x = 0 in (0.12b), the approach i s the same and the details can be worked out. Finally, for the case f£—^.Q , the expansion techniques of this chapter can be applied to an n-dimensional free-boundary problem providing the system i s rotationally invariant in n-dimensions. Then the equations are, with H(r)~. 0 . With the change of variables «4 = *— the equations become - 80 -where and CF l&\=^ 4- • This problem has the same structure as those already studied and a s i m i l a r a n a l y s i s can be c a r r i e d out. CHAPTER I I I : NUMERICAL METHODS The ob j e c t i v e of any numerical c a l c u l a t i o n f o r the free boundary problem i s f i r s t to locate the boundary CFXXl at time X and then to represent the s o l u t i o n <fyte-tx} at that time. For the method presented i n t h i s chapter we accomplish t h i s objective by the following steps . ( S e e Figure 3 ) (a) Given the requirement to solve the problem oh [0,T] , choose a p a r t i t i o n X^T^ .^-- ^ T^ where T , = Q and X^ = RY* . Let (b) P i c k an R-parameter family of C lc>^«,fl,) curves ~ • Then represent cr cm C«,»*TP J by , where < X H i s the c h a r a c t e r i s t i c function with support on I and a^ th. th i s the value of the i parameter i n the j i n t e r v a l . The method fo r determining these values i s yet to be s p e c i f i e d . Let (c) For each I c a l c u l a t e a s o l u t i o n efy*%Lx,t X — X ^ / ) T O the heat equation on the domain - 82 -- 83 -Subject to the c o n d i t i o n s This task i s a p a r t i c u l a r example of the s o - c a l l e d i n v e r s e Stefan problem. The s o l u t i o n *pn i s the a n a l y t i c i n t e r p o l a t i o n f u n c t i o n to be used on the i n t e r v a l I . Note t h a t the f u n c t i o n has been de f i n e d n on the i n t e r v a l t«',T w — ~J ; thus f o r numerical purposes, <g?n need be c a l c u l a t e d only f o r small times. -(d) Determine the parameters '"CJ i by r e q u i r i n g At*) t h a t CT be c t o / P J . Here, use the r e s u l t s of (c) together w i t h the f l u x c o n d i t i o n at x = <5* to determine these v a l u e s . (More d e r a i l s l a t e r . ) To demonstrate the method we do numerical c a l c u l a t i o n s f o r both £.j-C4,o-l and %>^- <o~\ 1. HEAT FLUX SPECIFIED AT x = 0 ' The d i f f i c u l t step i n t h i s scheme i s step (c) s i n c e , i n ge n e r a l , no a n a l y t i c s o l u t i o n s to the i n v e r s e problem have been found. However, such i n t e r p o l a t i n g s o l u t i o n s can be c a l c u l a t e d f o r a s p e c i a l three parameter f a m i l y of curves. To f i n d these s p e c i a l s o l u t i o n s i t i s convenient to use the f i x e d boundary r e p r e s e n t a t i o n (0.13) and then so l v e the system subject to a l l the boundary conditions except (0.13c). Even though the domain i s rectangular, the equation cannot be solved by i n t e g r a l transforms i n y or X " because the c o e f f i c i e n t of " i f ^ i s a function of both y and T . To eliminate t h i s d e r i v a t i v e put ^ -XTe**' * •'•nen t n e equation becomes where J^cr 1 " . The d e r i v a t i v e can be eliminated by s e t t i n g %V — — — cj , so that W—— ~ ti"2- PCX'S where F i s an a r b i t r a r y function. The c o e f f i c i e n t of V then becomes F i r s t set - p . . The T-dependence can then be eliminated provided D i f f e r e n t i a t e t h i s expression with respect to T* to get ' Y ~ ~ Q Integrating three times, we have then J . The transformation of the dependent v a r i a b l e can now be written - 85 -where "VT s a t i s f i e s For the numerical c a l c u l a t i o n we choose a l i n e a r boundary cr— c*+ h>x , which i s contained i n the three parameter f a m i l y j u s t d e r i v e d . For t h i s case i t f o l l o w s from (3.1) that Cg— O and the equation f o r "IT reduces t o In .the k*"*1 time step I T then s a t i s f i e s the boundary c o n d i t i o n s v ci,x)^ o where Here a^ i s determined by c o n t i n u i t y - 86 -and b^ i s determined by matching the fl u x condition at " £ " ^ - 1 ' Since J^> has already been ca l c u l a t e d from the previous step and does not involve b^ , the equation i s as e x p l i c i t as i t appears. Now, omit reference to the time step. Return to (3.3) and put Then (3.3) becomes simply = • This can be solved by a Fourier transform i n y . With (ps-XyiT, h=. 1- 2 • - - a n d -\T /^ireoSV,,^ we have which has the s o l u t i o n Vhere ^ I c , ^ £ $ ^ 1 < ^ * V H . j . « ^ ^ ^ and X H H O - ..//GCPI e V j t ^ F I J F / i^i) - 87 - ; The i n i t i a l value integral (3.7) can be approximated by the simple quadrature where O e - O * * , h = £ Z £ and $ f e = ^ 0 { J . This is a simplified Filon quadrature which was found to be adequate for the present calculations. For N = 10 the quadrature was accurate to four significant figures over the range O ^ V > h ^ 2 5 * As for the boundary integral (3.8), note that the calculated solution i s to be used only for small time increments, hence small Since tT has already been approximated by a linear function in these time steps, there i s no significant added error in linearizing the boundary condition at x = 0 as well. Hence, represent G(£0 by the linear approximation The integration can then be carried out to give Finally, the solution to (3.3) has the expansion So from (3.2) the interpolating function (J* can be written in the original co-ordinates as 88 -This expression can be d i f f e r e n t i a t e d term by term to f i n d 4 ^ £ x l as r e q u i r e d i n (3.5). .fie With the above r e p r e s e n t a t i o n f o r the (p , the method can be checked against a Neumann s o l u t i o n s i m i l a r to (0.14) where H=. f ~ and X s a t i s f i e s the tr a n s c e n d e n t a l equation I n i t i a l i z e time at T = 1 and choose A— -ST and H = 1 . o Then £ = .6420 and 4^ i s a s o l u t i o n of the f r e e boundary problem w i t h the heat f l u x and the i n i t i a l c o n d i t i o n * -V 3C The corresponding f r e e boundary i s J . <rtxl= 0+ f ) *~ . For the numerical c a l c u l a t i o n equal time steps are chosen - 89 -and the temperature i s computed at f i x e d y = i n t e r v a l s . I n t h i s c o -ordinate the exact s o l u t i o n i s time independent and comparison w i t h the exact s o l u t i o n i s e a s i e r . See Table I I I . With a step s i z e of AX = .1 the c a l c u l a t e d temperature was accurate to w i t h i n 1% except near x = 0 . The e r r o r i s l a r g e s t at i t h i s boundary because the F o u r i e r s e r i e s r e p r e s e n t a t i o n f o r the tem-perature converges most s l o w l y . The e r r o r i n the f r e e boundary &~ i s l a r g e r ; i t grows as .03X. 2. TEMPERATURE SPECIFIED AT x = 0 For t h i s case, t o determine the i n t e r p o l a t i n g f u n c t i o n we must again s o l v e (3.3) f o r V . The boundary c o n d i t i o n s (3.4) are the same except that at y = 0 . V * U . - r i - ^ C T ^ x ^ l ) c r t t T i A - P M T I % H e r e a f t e r , omit reference to the time step. S o l v i n g the equation by F o u r i e r transforming i n y. leads to a s i n e s e r i e s r e p r e s e n t a t i o n f o r £ . Th i s s e r i e s , however, i s not s u i t a b l e f o r numerical c a l c u l a t i o n . At x = 0 i t does not converge to J^(o^-c') and converges o n l y very s l o w l y f o r x near zero. F u r t h e r , the s e r i e s r e p r e s e n t a t i o n f o r 4*^ at x = a + bx does not converge at a l l . These problems can be e l i m i n a t e d when (3.3) i s solved i n s t e a d by a Laplace transform i n the v a r i a b l e - 90 -A-lot-9- feu) The s o l u t i o n can then be w r i t t e n where and X^>v are defined i n Appendix I I . Here i d e n t i f y $p l t f> of the appendix w i t h P ( r ) . Consider f i r s t the boundary terra "^g v • Expand the integrand of (II.3) = r r -Sv t * f e * . o _ _ / p T and use r L . J P<40e e t o get the r e p r e s e n t a t i o n (3 ^) - 91 -As i n the previous problem the boundary condition at x = 0 can be l i n e a r i z e d . Here, using the approximation u e .leads to the problem of evaluating i n t e g r a l s of the form 1 3 f o r a = 0, 1 ; b = - — ,•-- . ' Rather than evaluate the i n t e g r a l s d i r e c t l y , use the i d e n t i t y A power se r i e s representation for "^J can then be used when \ j£ 5* . For A- •* 5" . a one term asymptotic expansion f o r J^^if. oo gives an accurate r e s u l t . ( See [17]. ) For the i n i t i a l value term a s i m i l a r c a l c u l a t i o n y i e l d s the representation 4 ^ j ^ ^ - — . L_ j IM04 t y(i- ! i)ie 4 r 4_ c j 1 4 e J i & . 4 - e ~~ J j To c a l c u l a t e » i n t e g r a l s of the type must to evaluated. For a l a r g e range of £ the F i l o n quadrature g i v e s good r e s u l t s . Here h = , and (pC^0 . The expansions and (3-1c>) together give an expansion f o r I T . The expansion f o r can then be found from C ? . 2 ) •:. With t h i s method of r e p r e s e n t a t i o n the numerical c a l c u l a t i o n was checked against two exact s o l u t i o n s . See Tables IVa, and IVb. (a) A Stefan s o l u t i o n w i t h <. 0 • Acr (b) A Neumann s o l u t i o n w i t h —.—* > o In (0.14) i n i t i a l i z e time at T = r and take the p a r t i c u l a r case - 93 -The c a l c u l a t i o n was done w i t h ..step .sizes £ > T .b.e.t.w.e„en .1 and .5 . Then w i t h at most two terms i n the s e r i e s (3.9) and (3.10) the accuracy f o r the temperature was .5% w h i l e f o r the f r e e boundary <T* an even higher accuracy of .1% was achieved. ' 3. DESCRIPTION OF A FREEZING LAKE The numerical technique can be a p p l i e d to a s p e c i a l two phase problem - the f r e e z i n g of a shallow i a k e . ( F i g . 4 ) A -lake i s considered shallow when during the p e r i o d of f r e e z i n g the e f f e c t of the l a k e bottom must be taken i n t o account. Let r\l> be the p e r i o d i n which f r e e z i n g takes p l a c e . Then there are two p o s s i b i l i t i e s : . • (a) I f SCT) O > c l e a r l y the bottom must be in c l u d e d i n the d e s c r i p t i o n . For the l a k e we study, t h i s i s not the case. (b) Even i f •SCT )^-— H the heat f l u x from the bottom of the l a k e may a f f e c t , through d i f f u s i o n , the water temperature near the top of the l a k e . - 94 -Figure 4 Cross-section o f a. shallow lake -- - 95 -Consider the second e f f e c t . In the l i m i t ^CT^)^ - H the water phase can be approximated by the fi x e d boundary problem a t K w • — -7-7 £)< X < H with the boundary conditions where 1*^ 6* i s the d i f f u s i v i t y and i s the conductivity of water. The phys i c a l heat fluxes are approximated by t h e i r time averages over the fr e e z i n g period. Let t~ = T T " =• <£ctr»£T-The s o l u t i o n can be wr i t t e n as where a ^ ^ - { ( 7 ) L i t 6 + 0 ••.V^ "* The e f f e c t of the bottom w i l l become s i g n i f i c a n t f o r a depth H such that k 1 • This w i l l hold only i f i t holds f o r the fe= terms i n the sums for 4 V . r and 4. . F o r " 7 ^ - £ t follows that - 96 s i r • « • * ~ ^ which can be r e w r i t t e n as For a given l a k e t h i s equation can be solved f o r C and the c r i t i c a l H determined. I f the l a k e i s shallower than t h i s H, i t can be considered a "shallow" l a k e . We apply the numerical method to an a n a l y t i c model of Seneca Lake i n northern Michigan [.19] . For t h i s l a k e , \. so from (3.11), C~*fl-S^ which i m p l i e s a c r i t i c a l depth of 10 meters. Seneca Lake i s 2 meters deep and so can be considered shallow. • To set up the model, f i r s t non-dimensionalize w i t h the water parameters M £r f°r both water and i c e * T * _ ^SZ. M. •• temperatures. For the i c e , the equations are 6* - 97 -When t h e d a t a for Seneca Lake was accumulated the i c e was covered by 20 cm. o f snow. As a r e s u l t the temperature $g at the top of the ic e remained c o n s t a n t throughout the free z i n g period. . Since e r i ^ ) - = - -ji , no i n i t i a l c ondition need be s p e c i f i e d f o r the i c e . ; . -• For the water the equations are: = >T ***** For Seneca Lake the stored summer heat provided a heat f l u x M K t x | ~ i So ^ i / ^ m ^ - c e e ' T h i s f l u x overshadows the ever--A m . Tit £ present geothermal gradient- of Jo /cwO-— ^e.<» * An i n i t i a l c ondition i s required for the water F i n a l l y , the f l u x condition at the interphase boundary CT ( .x1 i s - 98 -fey; 1 ^ 6 £ w where # j A two-phase . c a l c u l a t i o n where both the water and i c e temperatures are c a l c u l a t e d n u m e r i c a l l y leads to d i f f i c u l t i e s . With JA"^ - -J. i n (3.12) the i c e temperature propagates almost ten times f a s t e r than the water temperatures and so makes the c a l c u l a t i o n awkward. This d i f f i c u l t y can be turned to an advantage, however, f o r the i c e e q u i l l i b r a t e s ten times f a s t e r as w e l l . Thus, i t i s a good approximation to represent the i c e - s o l u t i o n by a p e r t u r b a t i o n -series i n and match i t to the numerical s o l u t i o n f o r the water at every time step. This r e s u l t s i n an equivalent one phase problem f o r the water w i t h a modified l a t e n t heat and e x t r a heating terms. Assume a p e r t u r b a t i o n expansion f o r (3.12) % - + • • • to get CT* * Therefore This expression can be used i n (3.13) to give the f l u x c o n d i t i o n § <5$<r ~f <4TW - 99 -where i s the e f f e c t i v e l a t e n t heat and i s the added heat source. This f l u x condition can now be used i n the numerical s o l u t i o n for the water temperature. temperature above the i c e , the f a s t e r the lake freezes over. * To generate higher orders i n K' , note that the expressions are to be used i n conjunction with a l i n e a r boundary f o r a given time step. Thus i n the perturbation c a l c u l a t i o n there i s no s i g n i f i c a n t added error i n s e t t i n g For the f r e e z i n g case , T J , <C O , hence Lj^ ^- A and O f o r small . Both of these e f f e c t s enter the f l u x condition so as to make more negative as «T* becomes more negative. In p h y s i c a l terms, the lower the c& T - 100 -APPLICATION TO SENECA LAKE [ l 9 l On December 2nd, after 25 days of alternate freezing and thawing, a permanent ice cover of approximately 15 cm. formed on Seneca Lake. The numerical calculation was i n i t i a l i z e d to this day. As the ice grew during the month of December, the temperature of the water near the bottom rose from 1°C (Dec. 2) to near 4°C (Dec 26). At 4°C, water is most dense so that any further heating from the bottom cf the lake results in convection rather than an increase in water temperature. This convection i s beyond the scope of the model; hence the calculation was terminated on December 26th. Further aspects of the calculation are the following: (a) The heat flux at the bottom was determined by s o i l temperature profiles beneath the lake. In the period* Dec 2 - 2 6 the total heat released to the water was 303 cal/cm^ to give an average heat flux of -4 2 1.45 x 10 cal/cm -sec. In dimensionless units, (b) Despite the large fluctuations in air temperature, a 20 cm. snow cover kept the temperature at the top of the ice very nearly a constant -1°C. In dimensionless units "T^cS — , (c) The calculation was started when the ice was 15cm. thick or in dimensionless units CTtol ^ - 1 2 5T - 101 -At this time the water temperature was linear in the depth; for the calculation, the in i t ia l temperature was taken to be ~r w 1 9 , ~ ) ^ . 0 4 g . c• 2 5 - tj) , (d) For the temperatures entering the calculation, the O f^O term in Lj* represented less than a 1% correction and so was dropped. In Fig. 5 the numerical calculation is compared with the experimental obervations. At depths of 50 cm. and 150 cm. the difference between the calculated temperature and the observed temperature was within 20% of the observed temperature. Only this rough agreement can be expected because the instantaneous heat flux was not available and a heat flupc averaged over the freezing period was used. Finally, Bilello reported that by Dec. 26 the ice was 20-25 cm. thick. The calculated value of 22 cm. is within this range. - 102 -APPENDIX I ' fo r r ~ J - s | - T ^ «*» The function FtT? i s assumed to be l o c a l l y i htegrable on (of o»7 and to have an asymptotic expansion of the form P ( T I ~ TL F ^ T ' * * d . l ) fo r T - ^ o o , % CASE 1 A * ° C ° f o r T - > < * > ; . . The expansion i s c a r r i e d out i n ,two steps. Ca) F i r s t expand the i n t e g r a l f o r f i x e d T , A *r o where \ Q(. ( f O | i s an asymptotic sequence f o r A — > o Cb) Then expand the I f ^ f.«yT - l f o r X — V Oo where 'CH-. lT l l i i s a n asymptotic sequence f o r X—> «** , Step Ca) With X f i x e d at t h i s step, put -4"(*«):=.—• „ P C 1 ( 1 4 - * *•*.•'/ I t i s then required to c a l c u l a t e the expansion of the Laplace transform - 1 0 3 -w ( 1 . 2 ) f o r s m all A . We i n d i c a t e how the work of Handlesman and Lew [17] can be used to complete t h i s step. Brealc up 4*- at i*= x* a n ^ put -fuo 9U—*)-4- ^ 9 t « - 0 Then H 4. \.-%~\ — f U.^'* (W"J Jt4 defines a f u n c t i o n a n a l y t i c f o r Rt ^Jr y (X- f o r some o <• ex. ^ 4. . Even i n the case' 0.= , the ' • m ^ "fj Cit3 e x i s t s l-e- M 4* t^-^ does not have a pole at — \ Now f o r u oa , " T l * * V has the asymptotic expansion * H = c where l-T} fe-*". ^ ( M (1.3) and frw + V" » n o t e t n a t no i s an i n t e g e r , I t f o l l o w s that ti£%ttJ= Xj. ^ttO<*** d e f i n e s a f u n c t i o n a n a l y t i c f o r <. S0 which can be a n a l y t i c a l l y continued to a meromorphic f u n c t i o n w i t h simple poles at x.= $ and residues 0%^ . - 104 -Define HJ = H ^ l * } 4- ' H + ^ t O . I f G. S e , then the i n t e g r a l representations of -^f^ and are defined on the common v e r t i c a l s t r i p S for which The M e l l i n transform then e x i s t s i n t h i s s t r i p S and so that H^ ,4/ . i s the a n a l y t i c continuation of If • S^.'C <v , so that 5 i s empty , the function s t i l l plays the r o l e of the M e l l i n transform. Hereafter, the d i s t i n c t i o n between H^-y and M f i s suppressed and the notation H-f i s used.. . With t h i s d e f i n i t i o n of the M e l l i n transform, the Parseval theorem f o r M e l l i n transforms i s v a l i d : JF ^ ) 9 t a v a t i = / A HI fri H 9 r i - * 3 ^ . : . . Choose . 5 1 M U e 7 ; then P.Ct-V) A and the Bromwich contour A i s to be taken i n the s t r i p \ < ^ The contour A can then be displaced to the r i g h t to pick up the poles of H-f CxT and P {4-%} ', since no i s an integer these poles are simple. This gives an asymptotic expansion f o r (1.4) - 105 -where the and S m are given by (1.3). Even i n the l i m i t 0L—> dL the r e s u l t i s v a l i d . The f i r s t sum of (1.4), which depends only on the asymptotic behavior of •£ w i l l be c a l l e d the ASYMPTOTIC c o n t r i b u t i o n : the second sum w i l l be c a l l e d the DOMAIN c o n t r i b u t i o n . F i n a l l y , using the d e f i n i t i o n o the r e s u l t s (1.3) and (1.4) can be c o l l e c t e d to give «»^«, **** (I.5a) where A F C r , x , t*»1— *— I r txl (I.5b) and ^ F ^ x , ^ 1 i s the a n a l y t i c c o n t i n u a t i o n i n of X-1 (I.5c) This completes step ( a ) . - 106 -Step (b) The asymptotic expansion of AP (_*"",T, ^ ~) f o r T~—*f 0 0 f o l l o w s simply by s u b s t i t u t i n g (1.1) i n t o (I.5b) (1.6) The expansion of PP ( r , x t t w + i*} f o r X—>PO i s more d i f f i c u l t . For convenience, transform the i n t e g r a l (I.5c) to the i n t e r v a l (0,1) w i t h the t r a n s f o r m a t i o n 1*4^ ' *"**" This i n t e g r a l i s of the general form I J(xi= i 3U0 F ( x t*i a t « (1.8) where <j has a convergent power s e r i e s expansion f o r | u l < ±_ 0a . - 1 ^ 1 . (1.9) Now the M e l l i n transform of ty •> u 3 it*i an i s a meromorphic function with poles at - £ = — te. > fesi o , 1, 1 • • - • Tn f a c t , use the s e r i e s representation (1.9) i n the transform; since the expansion i s absolutely convergent, i n t e g r a t i o n and summation can be interchanged to get M 9 L V U E T (i.io) Define then, by s u b s t i t u t i n g (1.9) into (1.8) we have JCxl= H 9 f c M F C T , fe+i) ( i . i 2 ) while from (I.11) we get 7 r l T MKT, tvO J= T Fui # ( t . 1 3 ) Now use the expansion (1.1). I f there i s a ^ such that J^-s 4 X denote i t by ; then from (1.13) MFU, M i ) . T ^ _ ! _ -.108 -with the understanding that i f there i s no such the term p —FR1 ^ i s absent. Here Hg^ i s an i n t e g r a t i o n constant. Let Q — £<9-gj.} _ Then from (1.10), (1.12) and (1.14) we get the f i n a l form ^ T « I _ T- % M 9 v» + n % ^ * + C 9 F C where i f no r~ i s an integer then G$. i s empty and the second sum *-i s absent. This r e s u l t can now be applied to (1.7) with the i d e n t i f i c a t i o n 1 4 - " / u so that where g ^ w-j ^ H l x l P f j u l i s the Beta function. The a of (1.9) can be read o f f from the expansion te=t> - 109 -to give f i n a l l y 6 , . . F ^ 1 Id-") *• n rr 1) The expressions (1.6) and (1.15) are the r e q u i r e d expansions to complete step (b). These expansions can now be used i n (I.5a) to complete case 1. CASE 2 A - O C 1 ) . f o r X—><*> W r i t e A~- A„4- ApCxI where A*, i s a p o s i t i v e constant and Ap, (xi— -> o f o r X—>- oo • As i n case 1. the expansion i s c a r r i e d out i n two steps: S t e P Q > With X f i x e d , put 4lMl= C | - ^ - ^ r - F C l T T l ) * I t i s then r e q u i r e d to c a l c u l a t e the expansion of the Laplace transform f o r Ap.—> © . As i n case 1, H 4^ and H-f be d e f i n e d . < X can 110 -Then, because -§r i f e x p o n e n t i a l l y decreasing f o r U ^ fy\ i s e n t i r e and M - f - H T l + i s a n a l y t i c f o r ^ 5. • The contour fe> i n the P a r s e v a l theorem can then be s h i f t e d to the r i g h t . There i s no asymptotic c o n t r i b u t i o n ; the r e s u l t i s simply L f t A ^ l ~ H H | • ( I 1 7 ) t*»l Note that the same r e s u l t could have been obtained simply by expanding the exp o n e n t i a l i n (1.16). F i n a l l y , from the d e f i n i t i o n of -f. and (1.17) i t f o l l o w s t h a t rx» I** i F l r . T . J O - II t>FCr , x , ^ 1 ) ( I 1 8 a ) where f>P K T , » » • 1 , = j ^ ^ — ^ . T ^ ^ " ) e, (I-18b) Step (b) The expansion of Pf= (rf TT, t*vf-0 f o r T—^-«*> i s c a l c u l a t e d as i n case 1, step (b). Here, a f t e r transforming the i n t e g r a l (I.18b) to the i n t e r v a l C<», i } we get * * ' . ' . so that - I l l -where i s the second confluent hypergeometric f u n c t i o n . The c a l c u l a t i o n of the power s e r i e s expansion 9 * L a i - C S^,,*, u t e i s s t r a i g h t forward, but t e d i o u s . To the order which the matching of the i n t e g r a l equations i s c a r r i e d out, we need o n l y that v= d. F i n a l l y , the r e s u l t analogous to (1.15) i s This completes step (b). When (1.19) i s s u b s t i t u t e d i n t o (1.18a) case 2 i s completed. - 112 -APPENDIX I I The s o l u t i o n to the system can be found by using the Laplace transform. The s o l u t i o n has three co n t r i b u t i o n s : (a) The i n i t i a l value term which s a t i s f i e s Transform i n X" and in v e r t to give where and P i s a Bromwich contour f or which |> > c? * (b) The source term which s a t i s f i e s . _ ' . <H - ,X j = 4>(*,Xl*=o Then with £ jT 5 t«j ,x> & <1 X i t follows that - 113 -(c) The boundary term which s a t i s f i e s w i t h $ ^ t p ) - J " ^ g l x l e rftr the s o l u t i o n i s ^ v ^ . x » % i r . J r $ l P ) — — — ^ d p ( I I 3 ) With the s o l u t i o n given i n i n t e g r a l form by ( I I . 1 , I I . 2 , II. 3 ) both the l a r g e X and sm a l l X expansions f o r can be found. - 114 -TABLE I.a P e r t u r b a t i o n S o l u t i o n to a Stefan Problem w i t h boundary c o n d i t i o n s __ <J>(0,T) •» e~ T s i n ( l ) , . . s i n ( x ) " s i n ( / 2 x) , •Kx.O) = sin (1-x) + E [ ^ T J J - S I N ( ;/2 ) J ' e » .1 Comparison w i t h a numerical s o l u t i o n ; temperatures compared at f i x e d y=x/a i n t e r v a l s P e r t u r b a t i o n Soln. Numerical Soln. Free Boundary x=.25 1.0219 1.0218 Temperature .0864 .0811 y=.22 .3349 .3294 y=.55 •4816 .4779 y=.7.7 Pree Boundary | x=.50 1.0386 1.0384 Temperature .1337 .1254 y=.22 .2608 .2579 y=.55 .4256 .4248 y=.77 Free Boundary x=.75 1.0514 1.0513 Temperature .1041 .0990 y=.22 .2491 .2487 y=.55 .3315 .3318 y=.77 - 115 -TABLE I.a (CONTINUED) P e r t u r b a t i o n S o l n . N u m e r i c a l S o l n . F ree Boundary T = 1 . 0 1 . 0 6 1 3 1 . 0 6 1 3 • Temperature . 0 8 1 1 . 0 7 7 9 y = . 2 2 . 1 9 4 0 . 1 9 4 5 y= . 5 5 . 2 5 8 I . 2 5 9 0 y = . 7 7 F r e e Boundary T = 2 . 0 1 . 0 8 3 0 1 . 0 8 3 2 Temperature . 0 2 9 8 . 0 2 8 5 y= . 2 2 . 0 7 1 4 . 0 7 1 3 y = . 5 5 . 0 9 5 0 . 0 9 5 2 y = . 7 7 F r e e Boundary x = 3 . 0 1 . 0 9 0 7 1 . 0 9 I . O Temperature . 0 1 7 1 . 0 1 6 1 y= . 2 2 . 0 2 7 6 . 0 2 7 7 y= . 5 5 . 0 3 6 7 . 0 3 6 9 y = . 7 7 - 116 -TABLE I.b Perturbation Solution to a Stefan Problem with boundary conditions _. <J>(0,T) - e ~ T s i n ( l ) , * r . j /1 \ i r sin(x) sin ( / 2 x) , • <x,0> = s i n ( l - x ) + e[ ~[^5J - s l n ( h ) ] e « .5 Comparison with a numerical s o l u t i o n ; temperatures compared at f i x e d y=x/a i n t e r v a l s P e r t u r b a t i o n Soln. Numerical Soln. Free Boundary *=.25 1.0883 1.0914 Temperature .1574 .1496 y=.22 .3967 .3927 y=.55 .5390 .5376 y=.77 Free Boundary t x=. 5 1.1581 1.1617 Temperature .1301 ,1227 y=.22 .3203 .3131 y=.55 .4286 .4253 y=.77 Free Boundary x=.75 1.2129 1.2165 - • .1059 .1004 y=.22 .2563 .2505 y=.55 .3387 .3357 y=.77 - 117 -TABLE I.b (CONTINUED) P e r t u r b a t i o n Soln. Numerical Soln. > Free Boundary T = l . 1.2559 1.2596 Temperature .0852 .0815 y=.22 .2038 .1998 y=.55 .2667 • .2646 y=.T7 Free Boundary i T= 2 . 1.3523 1.3562 '. Temperature .0336 .0321 y=.22 .0784 .0774 y=.55 .1006 .0996 y=.77 Free Boundary x=3. 1.3869 1.3915 Temperature .0133 .0132 y=. 22 .0308 .0311 y=.55 .0392 .0392 y=.77 - 118 -TABLE I I Perturbation Solution to a Stefan Problem with boundary conditions <K0,T) - 1. + .25 sin ( T ) <Xx,0) = 1. - ••erf(/OT2)_x) erfOTxTT) ) i e =» .1034 X = .1 Comparison with numerical s o l u t i o n f o r O ( T ) Time X Free Boundary O(TT) P e r t u r b a t i o n Soln. Free Boundary O(-B) Numerical Soln. 1.0 1.1060 1.0995 • 2.0 1.2132 1.2035 3.0 1.3042 1.3000 4.0 1.3726 1.3777 5.0 1.4270 1.4375 1 6.0 1.4839 1.4912 7.0 1.5532 1.5523 8.0 1.6304 1.6242 9.0 1.7020 1.6991 10.0 1.7588 1.7656 . - 119 -TABLE I I I Numerical S o l u t i o n to a Stefan Problem w i t h Boundary Conditions < o , t ) = - - 1 * ( x , o ) = s1 e - ' 2 ; 7 4 at X - /l +. T X I n t e r p o l a t i n g temperature found by FOURIER TRANSFORMATION Time Step .1 Number of terms i n F o u r i e r Representation 4 Number of nodes i n F i l n n Quadrature 10 Temperature compared w i t h exact s o l u t i o n at f i x e d y=x/o i n t e r v a l s Numerical Soln. Exact Soln. Free Boundary x=.5 1.239 1.225 Temperature .178 .182 y=.22 .483 .485 y=.55, .697 .701 y=.77 Free Boundary x=1.0 1.438 1.414 Temperature .178 .182 y=.22 .483 .485 y=.55 .697 .701 y=. 77 Free Boundary x=1.5 1.613 1.581 .178 .182 y=.22 .482 .485 y=.55 .696 .701 . y=.77 - 120 -TABLE IVa Numerical S o l u t i o n to a Stefan Problem w i t h Boundary Conditions • ( 0 , T ) . - . 1 . - e T - 1 <Kx,0) - 1.- e * - 1 1 I n t e r p o l a t i n g temperature found by LAPLACE TRASFORM Time step .25 Number of terms i n Laplace Representation 2 Number of nodes i n F i l o n Quadrature 10 Temperature compared w i t h exact s o l u t i o n at f i x e d y=x/c i n t e r v a l s Numerical Soln . Exact S o l n . Free Boundary x=.25 .7500 ,7500 Temperature .1554 .1476 y=.22 .2905 .2835 y=.55, .3623 .3588 y=.77 Free Boundary x=.50 .4958 .5000 Temperature .1016 .1052 y=.22 .2502 .2425 y=.55 . 3279 .3222 y=,77 Free Boundary •c=.75 .2506 .2500 Temperature .0517 .0540 y=.22 .1223 .1297 y=.55 .1636 .1767 y=.77 - 1 2 1 -TABLE IVb Numerical Solution to a Stefan Problem with Boundary Conditions • <xf0) - j f e ~ ^ / 4 d 5 1 - C2/A *(0 , x ) = / e .« Interpolating Temperature found by LAPLACE TRANSFORM Time Step .5 Number of Terms i n Laplace Representation 2 Number of Nodes i n F i l o n Quadrature 10 Time X Free Boundary o(x) Numerical Soln. Free Boundary O(T) Exact Soln. 1.0 1.0958 1.0954 2.0 1.1829 1.1832 3.0 1.2647 1.2649 4.0 1.3418 I . 3 4 I 6 5.0 1.4150 1.4142 6.0 1.4848 1.4832 i 7.0 i l 1.5516 1.5492 8.0 1.6157 1.6125 9.0 1.6775 1.6733 - 122 -WATER TEMPERATURE 5 h OBSERVATION CALCULATION 150 cm Depth 12 15 18 21 24 (DAIS) F i g u r e 5 C a l c u l a t e d and Observed Temperature P r o f i l e s f o r Seneca Lake; December -2 to December 26 - 123 -BIBLIOGRAPHY [ I ] K. STEFAN, Ann. Phys. u. Chem. 42 (1891) 269 [2] F. NEUMANN, Lectures c. 1860. See Die P a r t i e l l e n D i f f . Der Math-. Physik 2 (1912) 121 -[3] H. CARSLAW & J . JAEGER, Conduction of Heat i n S o l i d s ; Clarendon P r e s s , Oxford 1969 [4] B. BOLEY, Heat Conduction Problems, J . of Math & Phys., 40 (1961) 301 [5] H. LANDAU, Heat Conduction i n a M e l t i n g S o l i d , Quart, of A p p l i e d Math., £ (1950) 81 [6] G. EVANS, H. ISAACSON & T. MACDONALD, Stefan L i k e Problems, Quart, of A p p l i e d Math., 8_ (1950) 312 [7] J . CANNON & C. DENSON HILL, E x i s t e n c e , Uniqueness, S t a b i l i t y and Monotone Dependence i n a Stefan Problem f o r the Heat Equation, J . of Math. & Mech., 17 (1967) 1 [8] R. HANDELSMAN & J . LEW, Laplace Transforms i n the Neighbourhood ' of the O r i g i n , Siam J . of Appl. Math., _1 (1970) 118 [9] A. FRIEDMAN, Free Boundary Problems f o r P a r a b o l i c Equations, J . of Math. & Mech., 8 (1959) 499 [10J B. SHERMAN, General One Phase Stefan Problems f o r the Heat Equation, Siam J . of Appl. Math., 20 (1971) 555 (June) [ I I ] B. SHERMAN, L i m i t i n g Behavior i n Some Stefan Problems as the Latent Heat Goes to Zero, Siam J . of Appl. Math., 20 (1971) 319 (March) [12] H. HASHEMI & C. SLIEPCEVICH, A Numerical Method f o r S o l v i n g Two Dimensional Problems of Conduction With Change of Phase, Chem. Eng. Symp. S e r i e s , 63_ (1967) 34 [13] G.H. MEYER, M u l t i d i m e n s i o n a l Stefan Problems, Siam J . Num. Ana l . 10 (1973) 520 (June) - 124 -[14] Y.K. CHUANG & J. SZEKELEY, Use of Green's Functions for Solving Melting or Solidification Problems, Int. J. Heat & Mass Trans., 14 (1971) 1285 [15] J. DOUGLAS & T. GALLIE, On the Numerical Integration of a Parabolic Differential Equation Subject to a Moving Boundary Condition, Duke Math. J., 22 (1955) 557 [16] G. BLUMAN, Applications of the General Similarity Solution of the Heat Equation to Boundary-Value Problems, Quart.Applied Math. 31^ Q9741 403 [17] M. ABRAMOWITZ & I. STEGUN, Handbook of Mathematical Functions, Nat. Bureau of Standards - Appl. Math. Series, (1964) [18] G. DUFF & R. NAYLOR, Differential Equations of Applied Mathematics, N.Y. (1966) [19] M. BILELLO, Water Temperatures in.a Shallow Lake During Ice Formation, Growth and Decay, U.S. Army Cold Regions Research and Engineering Laboratory, Dec. 1967 [20] L. SIROVICH, Techniques of Asymptotic Analysis, Springer Verlag (1971) [21] S. CHURCHILL, Modern Operational Mathematics in Engineering, N.Y. (1944) ,
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Approximation techniques for the Stefan problem
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Approximation techniques for the Stefan problem Katz, Hart Victor 1974
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Approximation techniques for the Stefan problem |
Creator |
Katz, Hart Victor |
Publisher | University of British Columbia |
Date Issued | 1974 |
Description | The macroscopic description of matter undergoing a phase change (the Stefan Problem) can be formulated as a set of coupled, non-linear partial differential equations. For the case of one space dimension, the thesis develops three approximation methods to solve these equations. (a) Asymptotic Expansions With the Green's functions to suit the given boundary conditions, the system can be transformed into a set of integral equations. For the case where the initial phase grows without limit as т → ∞ ,the large т expansions for the integrals are calculated and the large т behaviour of the interphase boundary found.(b) Perturbation Expansions When the latent heat of fusion of a material is large relative to the heat content of that material, a uniformly valid perturbation expansion in a parameter related to the ratio of these heats is possible. The first few terms of the expansions for the temperature distribution and the position of the interphase boundary are calculated and found to be in good agreement with the few known exact solutions and numerically calculated solutions. (c) Numerical Techniques Rather than use a traditional finite difference scheme, only time is discretized and an analytic expression for an approximate temperature is found for each time step. This method gives good results for the temperatures over the time intervals required by the physics of the problem. Finally, the method is applied to describe the freezing of a shallow lake. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-28 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080116 |
URI | http://hdl.handle.net/2429/19209 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1974_A1 K38.pdf [ 4.84MB ]
- Metadata
- JSON: 831-1.0080116.json
- JSON-LD: 831-1.0080116-ld.json
- RDF/XML (Pretty): 831-1.0080116-rdf.xml
- RDF/JSON: 831-1.0080116-rdf.json
- Turtle: 831-1.0080116-turtle.txt
- N-Triples: 831-1.0080116-rdf-ntriples.txt
- Original Record: 831-1.0080116-source.json
- Full Text
- 831-1.0080116-fulltext.txt
- Citation
- 831-1.0080116.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080116/manifest