Citation |

- Permanent Link:
- http://ufdc.ufl.edu/AA00038052/00001
## Material Information- Title:
- Source-and-sink method of solution of moving boundary problems
- Creator:
- Akbari, Mehdi, 1954-
- Publication Date:
- 1993
- Language:
- English
- Physical Description:
- xii, 229 leaves : ill. ; 29 cm.
## Subjects- Subjects / Keywords:
- Boundary conditions ( jstor )
Heat ( jstor ) Heat flux ( jstor ) Heat transfer ( jstor ) Liquids ( jstor ) Melting ( jstor ) Solids ( jstor ) Subroutines ( jstor ) Surface temperature ( jstor ) Temperature distribution ( jstor ) Dissertations, Academic -- Mechanical Engineering -- UF Mechanical Engineering thesis Ph. D - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt )
## Notes- Thesis:
- Thesis (Ph. D.)--University of Florida, 1993.
- Bibliography:
- Includes bibliographical references (leaves 223-228).
- General Note:
- Typescript.
- General Note:
- Vita.
- Statement of Responsibility:
- by Mehdi Akbari.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for non-profit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
- Resource Identifier:
- 029939746 ( ALEPH )
29713762 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

SOURCE-AND-SINK METHOD MOVING BOUNDARY OF SOLUTION OF PROBLEMS BY MEHDI AKBARI A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1993 TO MY WIFE VAHIDEH, AND DAUGHTERS SARA AND MONA ACKNOWLEDGEMENTS committee, study am most grateful Dr. C. K. Hsieh, as well as many hours to the chairman who provided of guidance of my supervisory the inspiration for this and encouragement. His enthusiasm has been contagious, and I hope to emulate his attitude in the future. I wish to express my sincere thanks to Drs. R. A. Gater, G. Emch, counsel R. Abbaschian, and understanding and H. A. Ingley have been essential whose unselfish in this study. Last but by Vahideh, no means least, for her love and spiritual I would support like to thank during my wife, this period of time of work on the dissertation. iii G ï¿½ good TABLE OF CONTENTS Page ACKNOWLEDGEMENT LIST OF TABLES LIST OF FIGURES NOMENCLATURE ... ABSTRACT CHAPTERS * . . 0 0 0 0 0 * 0 * 0 0 a 0 0 0 0 0 0 0 * 0 0 0 * * * 0 0 * t 0 6 0 0 0 * . 0 *0 0 0 0 a 0 * *60 0 0 0 0 0 0 0 0 0 0 0 0 * 0 * 0 0 9 0 * * 0 * 0 0 * 0 0 0 * 0 0 0 0 0 6 S S S 0 * 0 * 00 0 0 0 0 0 0* S 0 0 0 S 6 * * 0 0S 0 0 0 0* * . 0 0 0 0 0 0 0 0 0 0 0 0 9 0 0 0 0 0 * 0 0 0 0 * * * 0 * *0* 0 0 0 0 0 0& 0 * 0& 000 0 * * 0 0 0 0 0 0 * 0 0 & 0 O* 0 *0 *0*a 00 0 *00a 0 0 0 0* 0*0O* * 0 0 000 0* INTRODUCTION ............. LITERATURE REVIEW 0 0 0 0 0 0 0 9 0 & 0 0 0 0 *0 0 & 0 050* 0 0 0 0 ** *0 0 0 0 0 0 *00 * 0 0 0 18 .............1 SOLUTION OF INVERSE STEFAN PROBLEMS BY SOURCE-AND-SINK METHOD ............... Motivation ........... General Analysis ..... Pre-melt Stage ..... Melting Stage ...... Example Problems ..... Numerical Solution ... Critique of The Method Results and Discussion 0 0 * 00 00 0 0 * 0 0 0 0 0 0* 0 * 05*0 0 0 0 0 0 *00*00 0 00 0*S0 3.7 Extension and Concluding Remarks * 0 0 0 0 9 0 0 0*0 0 0 0 SOLUTION OF ABLATION PROBLEMS WITH ONE MOVING BOUNDARY BY A SOURCE-AND-SINK METHOD . . . . . 0 . a * 0 * * * 0 0 . 0 0 . 0 0 * 0 0 . 0 0 0 . . . . 4.1 Solution Methodology ........ Pre-ablation Stage ........ Ablation Stage ............ Equivalent Problem ........ 4.2 Illustrative Examples ....... 4.3 Numerical Solution of Ablated 4.4 Results and Discussion ...... 0 0 0 0~ a 0 0 0 0 Front S0 0 0 0 SOLUTION OF ABLATION PROBLEMS WITH TWO MOVING BOUNDARIES BY A SOURCE-AND-SINK METHOD . . . ... ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 iv . iii . .vi .vii * .ix * .xi I II III 3.1 3.2 3.3 3.4 3.5 3.6 IV .21 .21 .25 .28 .28 .31 .33 .39 .40 .56 0 9 0O9 0O 0 0 0 0 0 0 0O0 .58 .58 .59 .62 .63 .65 .67 .69 0 ~ . . 0 0 0 ~ . S ~ 0 0 ~ . . S 000 VI APPENDICES A 5.1 General Analysis 0... Pre-melt Stage 0600 Melting Stage ..... Ablation Stage Equivalent Problem 5.2 Examples ............ 5.3 Numerical Solution of Problem ............ 5.4 Numerical Solution of and Energy Storage 5.5 Numerical Examples 00 DISCUSSION AND RESULTS 006 EXTENSION OF THE SSM ..... 0~~~~ 0a 0 0 O 0 0 0 0 000 0 0 0 The Ablation Temperature 0 0 0 6 0 0 0 0 0 0 0 0 000000000000 0 0 0 0 0 0 0 6. 0 ï¿½ 0 0 0 0 0.* 0 0 6ï¿½ï¿½ 0 0 0 0 0 0ï¿½* 0 0 0 0 0 0 00 0 0 0 0 60 0 00 00 0 0 Page S..91 S.ï¿½91 .ï¿½ï¿½94 ï¿½..95 ..98 ..101 ..104 ..109 ..110 ..130 ....60... 0.135 B SSM FORTRAN PROBLEM . C SSM FORTRAN PROBLEM .. D SSM FORTRAN PROBLEM .. E SSM FORTRAN PROBLEM .. F SSM FORTRAN G SSM FORTRAN REFERENCES ............. BIOGRAPHICAL SKETCH .... PROGRAM PROGRAM PROGRAM PROGRAM 000 0ï¿½0 0 PROGRAM PROGRAM PROGRAM FOR FOR FOR o 0 0 6 FOR FOR FOR INVERSE INVERSE" INVERSE 0 00 0 0 0 INVERSE 600*00*0 INVERSE . . . ï¿½ï¿½0ï¿½0 ' "0 STEFAN 0ï¿½ï¿½ï¿½ï¿½FAN STEFAN a 0 0 0 0 0 0 0 STEFAN STEFAN 0 0 . . ï¿½.0 0 ï¿½0 0 ï¿½ 0 0 0 0 0 6 ï¿½ 0 6 0 .136 ï¿½*.000000 000147 0 0 . 0 0 0 0 0 0 * .154 0 0 0 0 0 0 ï¿½ 0 0 0 " ABLATION PROBLEM ... COMBINATION PROBLEM 0 .0 0 0 0 0 00 0 0 0 0 60 0 0 0 0 0 00 0 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 06 0 0 0 0 0 0 60 0 0 0 0 0 0 * 0 0 0 0 0 0 0 0 0 0 a 0 0 0 0 0 0 0 * 0 0 0 .165 .172 .186 .223 .229 LIST OF TABLES Page Table Expressions to account for Properties Conditions Comparison for first effects of aluminum ............... tested in eight examples .. between true and retrieved Stefan-Neumann example solve of boundary 0 0 *0 0 0 0* 0 0 0000 0ï¿½ 0 0 0 * 0 0 ~ ~~ 0 * 0 ~0 0 0 0 00 conditions d for boundary Comparison between true and retrieved conditions for second Stefan-Neumann example solved for boundary flux then temperature ............... . 0 0 ï¿½ 0 ï¿½ ï¿½ . . . . . . . . 46 Comparison between true and for third Stefan-Neumann ex retrieved ample solve conditions d for boundary Comparison for fourth between true and retrieved conditions Stefan-Neumann example solved for boundary flux then temperature ............................ Comparison between true and retrieved for fifth inverse example problem sol conditions ved for Comparison between true and for sixth inverse example pi retrieved conditions roblem solved for ï¿½ . . . 0 0 ï¿½ ï¿½ ï¿½ .51 Comparison between true and retrieved conditions for seventh inverse example problem solved temperaturet....................ween.trueï¿½andï¿½ret Comparison between true and retrieved conditions for eight inverse example problem solved for S 0 0 0 0 & 0 0 0 0 Conditions Conditions tested tested in four examples ..... ....... in three examples . . .. 0 0& 0 0 0 0 0 * * 0 0 0 0 0 0 500* 0 0 0 vi 3.1 3.2 3.3 3.4 3.5 .23 .41 .42 3.6 3.7 3.8 3.9 3.10 3.11 for 4.1 5.1 .53 .70 111 conditions ........................ temperature ....................ï¿½.ï¿½ ..ï¿½.ï¿½ ..ï¿½.ï¿½ ..ï¿½.ï¿½..ï¿½.ï¿½..ï¿½.ï¿½..ï¿½ï¿½.. 43 temperature ...............................................ï¿½ ï¿½ï¿½ ï¿½ï¿½ ï¿½ï¿½ ï¿½48 ï¿½ . ï¿½ ï¿½ ï¿½ 49 ï¿½.........5 temperature .............................. ï¿½ï¿½ ...... heat flux ....................................... ï¿½ . . ï¿½ï¿½ ï¿½ ï¿½ 52 heat flux ...................................... LIST OF FIGURES Page Figure Linearization of solid-liquid interface position curves for the numerical solution ............ * . o o35 System analyzed ............................................61 Trend of ablated surface position for a medium initially at phase-change temperature imposed with a constant heat-flux condition ............73 Accuracy, stability and convergency test of the SSM in the solution of ablation for a medium initially at phase-change temperature imposed with a constant heat-flux condition ...............75 Stabili in the medium conditi ty and convergency test of the SSM solution of ablation for a subcooled imposed with a constant heat-flux o n . . . . . . . . . . . . . . . ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½ ï¿½.. .7 8 Comparison of ablated surface position with two methods in the literature for a subcooled medium imposed with a constant heat-flux condition ................................................. 80 Comparison of ablated surface position with the moment method in the literature for a subcooled medium imposed with a linear heat -flux condition ........................................... 83 Stability and convergency test of the SSM in the solution of ablation for a subcooled medium imposed with a linear heat-flux condition ................................................. 85 Comparison of ablated surface position with the moment method in the literature for a subcooled medium imposed with a quadratic heat-flux condition ....................................... 87 Stability and convergency test of the SSM in the solution of ablation for a subcooled medium imposed with a quadratic heat-flux condition .................................................89 vii 3.1 3.2 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 System analyzed .. . ........... ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½2 Figure 5.1 5.2 5.3 ....*. ** . . . .106 .............114 Trends of solid-liquid interface and ablated surface positions for a combination problem of subcooled medium imposed with a quadratic heat-flux condition ..... 0................. . Temperature distributions in the medium at different times during ablation for a combination problem of subcooled medium imposed with a constant heat-flux condition .............................120 Stability and convergency test of the SSM in the solution of ablation for a combination problem of subcooled medium imposed with a constant heat-flux condition .............................122 Stability and convergency test of the SSM in the solution of ablation for a combination problem of subcooled medium imposed with a linear heat-flux condition ................................124 Stability and convergency test of the SSM in the solution of ablation for a combination problem of subcooled medium imposed with a quadratic heat-flux condition ............................126 5.5 5.6 5.7 5.8 5.9 5.10 viii Linearization of solid-liquid interface and ablated surface position curves for the numerical solution .......................... Trends of solid-liquid interface and ablated surface positions for a combination problem of subcooled medium imposed with a constant heat-flux condition ......................... Trends of solid-liquid interface and ablated surface positions for a combination problem of subcooled medium imposed with a linear heat-flux condition ......................... 5.4 Overall accuracy test of the SSM in the solution of the combination problem of subcooled medium imposed with a constant heat-flux condition . . .............................. ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½129 Page System analyzed ... . . . .............................. ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½ï¿½93 .ï¿½ï¿½.o.ï¿½ ... 116 ï¿½.............00 09 *118 NOMENCLATURE Specific Equation heat Temperature function used for boundary condition Flux function used for boundary condition, Green's function, Heat flux elements in coefficient matrix Thermal conductivity Latent heat of fusion Latent heat of vaporization Normal unit vector Distance from sink or source to sense point Interface Stefan positions number Dummy variable for integration, interface speed Temperature, temperature vector Time velocity of the interface Position, source point Thermal diffusivity Dirac delta function, Kronecker delta ix Lf Ste R1, R2 Density Dummy variable for time Subscripts ai m,v Phase nd Superscripts change Upper time limit index Time index Transition Melt and vapor fronts Initial 1.2 N1, N2 time Dissertation Presented of the Universit Requirements y of Florida for the Degre in Partial me of Docto] Fulfillment of the of Philosophy SOURCE-AND-SINK METHOD OF SOLUTION OF MOVING BOUNDARY PROBLEMS By Mehdi May Chairman: Dr. Chung K. Hsieh Major Department: Mechanical Akbari 1993 Engineering Source-and-sink methods have been developed for the solution of inverse of ablation diffusion and Stefar problems, i problems. ablation Green ' s problems, functions and combination have been used and integrodifferential positions equations are solved equations by using are derived for the interface in the phase-change medium. local linearization These position and the boundary heat flux if it were treated as unknown. The results have shown to be accurate, convergent, and stable. The methods developed for the solution of the inverse diffusion problems have been used to find the boundary conditions for the i approaches: approach. conditions inverse without solution approach They and to be useful are solved by two a time incremental to find the boundary to be supplied at both sides efficient of the phase in that they require change fewer interface. equations to The methods be solved are for the xi and the temperature of the interface Stefan problems. a series Both have shown reliance on the flux information Abstract of to the Graduate School unknowns. solution can also be easily of the problems. In the ablation use of the source-and-sink and combination problems, method heat transfer to analyze is solved fixed imposed domain so that the original on the moving boundary is taken heat-flux condition to be the condition that is imposed on the solid-vapor is then used or liquid-vapor together interface. with the temperature This flux condition at the interface to solve for the interface position as well as the condition imposed on the fixed used boundary. in the temperature Finally, equation these positions to complete and conditions the solution. method has been used successfully in solving examples encompassing one-, two-, and three-phase ablation with the boundary of the medium imposed with constant, linear, and quadratic conditions. Numerical solutions for these examples as well as the trends of the moving boundary positions are discussed in detail. xii the in are The flux developed f or the The algorithms CHAPTER I INTRODUCTION Phase changes occurs when a body is exposed to a large heat flux so that it melts or vaporizes at the surface. With continuous heating, the melting or vaporization front moves inward with time. Since priori, boundary problems treatment the position of the phase-change this type of problem problem is commonl in the heat transfer find applications and processing, th in thermal iermal protection interface v known literature. energy is unknown as the moving Phase storage, of re-entry =change material vehicles in space technology, laser drilling and cutting in manufacturing, freezing and thawing in soils and foodstuff, and photocoagulation opthalmological procedures, among others. It is difficult to solve a phase-change problem because of the presence such of the nonlinear a problem, the number condition of phases at the interface. encountered Usually for in the medium depends appears triple place. solid, occurs on the ambient in the medium. point pressure Thus for of the substance, On the other liquid, for and a medium hand, gas and the temperature a medium only for phase states heated sublimation change may appear of low thermal above range or cooled below and fusion the triple simultaneously, conductivity that the may take point, and this and exposed to In practice, in large heat flux. the changed phase may b e removed as in the re-entry where the substance is removed by a fluid dynamic force. This material removal may also be effected naturally such as in sublimation, where the sublimed In either vapor event, is dissipated the imposed heat by diffusion flux follows to the surroundings. the boundary motion, and such problems the literature. A different are commonly state referred of affairs to as ablation is encountered problems in when the changed phase remains stationary and continues to occupy the space taken its previous state. melt is not removed Thus for example by an external in melting force, a solid, it will adhere if the to the solid from which it changes phase. Then in the absence convection in the melt, the heat transfer in the medium can be analyzed first D by solving resented a Stefan a formal problem in solution recognition of the problem of Stefan in the literature Stefan that, (see problem next chapter is basically in the former, the boundary whereas, in condition the latter, for a survey different he problem is is always the imposed of the literature). from the ablation solved ir imposed o condition a fixed i the fixed follows th( problem domain, The in and boundary e boundary motion. diminishes the recedi have been Thus in the ablation in size with time, an( ng boundary numerous problem, d the impos, due to the shrinkage efforts in the literature the domain ed condition for analysis is acting of the medium. for the solution on There of the Stefan problems. Much less research, however, has been devoted the solution of the ablation problems. Also for those ablation that have been solved and documented by of who open to melted This soon as it is formed. occurs [1 problems in the literature, most moving study found are for the analysis boundary appears for problems of generic in the medium. involving at the solid-liquid two moving and liquid-vapor ablation There in which has been boundaries such interfaces. only a lack of as those They are encountered when the solid, liquid, and gas states appear simultaneously in the medium. Three dissertation. phase-change In the first problems prob1 em, will be an inverse solved Stefan in problem be solved. This problem differs from a conventional Stefan problem in the imposed sense on the that, boundary for the conventional is fully specified. problem, the condition This condition is then used together with others to solve for interface motion as well the temperature problem, however, distribution the interface in the medium. motion is given, For the inverse and this motion is as a part of the input to solve for the boundary condition is necessary found to provide in material for this interface treatment and processing, motion. where This problem the property the material may be a function of the phase-change rate at the interface. Solution of the problem thus enables a better control the property of the material. It should above be noted has recently literature. dissertation However is that the inverse been solved , the method more efficient by Stefan Zabaras developed because problem et al. [I1 for solution it requires described in the in this no information for the heat flux at the interface. As will be shown, two methods are developed for the solution of the problem: one expands and the other uses this will used as that is of of the one condition in polynomial or an inf inite series, the interface a time incremental position information at discrete times to solve for the coefficients in the polynomial or series expansion, while the latter uses the interface position at consecutive times to solve for the conditions incrementally. Both can be used to solve for the conditions accurately as supported by eight examples encompassing the search for constant and boundary. time-variant The methods temperature and are also well adapted to flux conditions a numeral at the solution. The second phase-change problem solved in this work is ablation problem that contains only one moving boundary. problem is solved fixed domain. In by a new method, 1 this method, which the ablated treats phase the problem in is considered remain moving in the space boundary can that is ablated be taken in the to be process. an interior As such, the phase-change interface, and the heat flux that is derived at this interface be used to match the condition imposed at the moving boundary. with the additional at this interface, the position prescription the flux condition of the interface of the temperature there for phase change can be used to determine as well as the hypothetical condition that is imposed conditions solution. at the fixed boundary. can be used in the temperature Finally these equation positions to complete Solution of the ablation problems has rarely been attempted a fixed domain. The method presented in this work thus represents new endeavor, which will be shown to be simpler than most methods described in the literature. To validate the method used in the and solution of the ablation an This a to can Then and the in The former utilIizes method. are provided four examplIes problems, they include with constant the analysis and of time-variant one- and two-phase flux conditions ablation at the imposed moving boundary. Of those four, the ones imposed with constant heat flux can be checked is thus selecte with an exact :d for checking solution. the soluti A one-phase on over the ablation entire p problem eriod of the ablation, testing while the solution a two-phase in ablation a quasi-steady problem state. is selected Both yield for accurate results stabi lity as shoi 4 tests later in this dissertation. have also been made rendering Convergence further assurance success of the method. above The good provide addresses combination results impetus ablation obtained in the second for the solution with of the ablation two moving and Stefan problem of the third boundaries. problems, described problem This an area that has rarely been studied in the literature yet is important in practice when solid, medium. moving liquid, Again boundary and gas the problem is treated states is solved appear in as an interior simultaneously a fixed domain, phase-change in the and the interface. flux condition at this boundary is again used with temperatures at the interfaces to determine the heat flux at the fixed boundary solution. as well as the interface Unlike its predecessor positions for which to complete the two simultaneous integrodifferential time step, equations for the combination are solved problem, simultaneously three equations for each will be solved three simultaneously. examples Again including the results imposition of are good as constant, evidenced linear, flux conditions. the and for that is The the by and wn quadrati c of all will be used for the solution the problems to a through in this work. development This method [2-5], where has recently the method has been subjected shown to be particularly addressed in suited i this wo for the solution ,rk. In this meth of the phase-change od, the problem is problems solved using one set of governing equation, initial, and boundary conditions, literature, equations, region. whereas in the conventional the problem is solved with each set focused The conventional methods methods by using ; developed different on the solution thus cannot of compete in the sets of one phase with the source-and-sink sink method, only it is in the solid assigned conventional methods, for simplicity. one temperature or liquid region equation depends in which in the source-and- will be derived; on the position This is in sharp several contrast equations whether that is to the must be derived, again one for each phase region. The time and effort saved in the source-and-sink method can thus be readily appreciated. With the inclusion of the three separate problems in this dissertation, approach. Thi each problem; coverage of methodology, felt that, the presentation of us one chapter the statement analysis, e only through of materials follows will be self-complete the problem, .xamples, and results such an organizatior a nontraditional of with the solution It is can be general and discussion. 1., the material presented in a coherent manner without severe fragmentation. what follows in the next chapter, a literature review will be given It is then followed by by method Moreover, in the equation. temperature will be devoted and the chapters for the presentation In A source-and-sink method the presentation of the three f irst. problems. conclude to and recommendations this dissertation. As usual, comp ilI ed all computer in the appendix. this research the future. will be useful programs It is hoped developed in this work that the insight in the analysis of the phase gained change are from in will be given conclusions Finally, CHAPTER II LITERATURE REVIEW Heat conduction boundary problems. These change problems constitutes a large class of to as Stefan Clapeyron problems [6], and were first solved the thickness in 1831 by of the solid Lame and crust by freezing of water under a constant temperature condi t ion. able to find the thickness to be proportional to the square root of time but did not find the constant of proportionality. Some thirty years later, Neumannn presented an exact solution to the problem his unpublished the Stefan solution w lecture problem notes, was solved and for the first in its entirety change in time in history, [7-8]. space Neumann's imposed a constant temperature condition. His work is important because, irrespective search for additional solved Neumann exactly today of the exact remain some one hundred numerous solutions those and thirty efforts . Stefan that years made by many problems were originally others that can be analyzed ago. While obtained the exact by Neumann, solution the problem to the Stefan p n has been named )roblem after was first Stefan in recognition of in 1889 [6-9]. who expressed 1 his published Practically the temperature work appearing the same problem of the phase-chan in the open was solved ge medium i literature by Stefan, n terms of moving with phase are commonly referred who determined They were with as for the phase in a semi-infinite in by a similarity variable, to be proportional and the position to the square of root of time. the moving Thus, in boundary essence, the moving boundary transformation. general Stefan problems, problem Since it has that was solved by the transformation now been firmly can be solved means of a similarity can not be used for established exactly, that, it must be for in an unbounded domain, of medium of constant properties, and imposed a constant temperature condition. severe limitation imposed by the similarity transformation, it is not surprising to see a wide variety approximate monographs what solutions developed have been prepared follows, reviewed; some popular interested readers in the literature. for presentation solution are referred of these methods to 10-14 will Books methods. be briefly for details. One of the early methods developed for the solution of the Stefan problems imposed with heat flux and convection conditions the power-series position expansion and the temperature method. In this method, in the medium are expanded the interface in terms of power series series [10] or complimentary are then substituted conditions to determine error functions into the governing the coefficients equation in the series [11-14]. These n and boundary expansion. Power series method works well for the short-time solutions. At that time, the series ai the series re sufficient converge r4 t to yield apidly, ai accurate ad only a results. few terms of The method also works for the solution in the neighborhood of singularities occur as a phase degenerates. For large time, more terms are needed in the series, With the with of and In is that may X14-t, and the method rapidly loses its appeal. 10 Stefan problems integrodifferential can also be solved equations. by reducing Lightfoot [15] solidification as a moving heat source front. With the use of Green' s equations function, he was able to derive one for the temperature two integrodifferential and the other for the interface position. Then by assuming the position a function of the square root of time enabled him to retrieve the Stefan-Nuemann solution. Lightfoot' s dimensional method Stefan has been extended problems involving to the solution phase change in of two an infinite wedge [16]. interface correction To avoid position term a tedious was taken was introduced iteration to be a hype into the solution in two dimensions, lrbolic function a to account the nd for the property variation in different phase regions. Indeed, Lightfoot' s method has been limited to phase change of equal properties. limitation has been lifted by Kolodner [17] by using double source and sink at the phase-change interface. Integrodifferential equations of Volterra or Fredholm also be derived in the solution of the Stefan problems by a Fourier transform inverted is useful [18-19]. in the solution. in providing This occurs when In general, an exact solution the transformed the integral in integral equations equation form. are method However, the integral equations must still be solved numerically solution. One of the unique occurrence continuously of multiple with time. features phase The boundary of the Stefan regions position whose problems domain of the domain is the changes is also of the solution. Boley them took to the This type can for introduced an sought as a part [20] 11 embedding a large method, fixed in which domain the time-variant for the solution domains are embedded of the problem. This introduces an unknown must be solved heat flux at the fixed together with the moving boundary, boundary and this flux position to complete the solution. Boley's embedding method has been extended to the solution of multidimensional problems without internal generation problems however, convenient [21]. solved It should be noted in this dissertation a source-and-sink method to use in the solution that all the moving also work in will be employed of the combination a fixed which boundary domain; is of Stefan more and ablation problems. Stefan [22-29]. I problems can also be solved n this method, by a quasi-steady an asymptotic solution expansion is derived by dropping to find derived the unsteady a long-time by dropping y temperature solution. A the interface terms in the governing quasi-stationary velocity terms solution in equations is also one of the governing short-time equations solution. and the interface Asymptotic expa flux condition mansions to establish of the temperature and the moving boundary position are then constructed by using these solutions for limits Asymptotic expansions problems with nonlinear melt [22]. Basically work best in the solution r conditions a a perturbation ind convective motions technique, it offers of Stefan of the insight into the physics tedious mathemati of the problem. cally; However, even the determination the method is very of the first-order terms in the temperature expansion has proven to be an intractable task [26]. in heat [26]. 12 Higher-order distinct terms drawback Coordinate become of the method. transformation increasingly offers difficult another to obtain, approach to solution embedding of the Stefan method, problems the transformation [30]. n works Somewhat in a fixed similar domain to the without the need of introduction of any new unknowns such as the boundary heat flux in the embedding method. In fact, the time-variant domain is mapped into an invariant domain in the transformation method; moving boundary is thus immobilized in the solution. solve It is noted the Stefan that use of the transformation problem. It only provides does not by the itself convenience of working in a fixed domain. The coordinate transformation has thus been used together accuracy solution with other methods and also in the numerical for the improvement methods for facilitation [31-37]. Like the Stefan problems problems. Stefan fixed also fall However, problem domain, problems into the general the ablation p because whereas briefly category )roblem in the latter in the forme is reviewed above, of the moving more complex the problem c< r the domain ablation boundary than the an be solved in is continuously changing with time. There have been numerous efforts developed the solution however, has bE of the Stefan een devoted to problems; the solution only very 1 i ini ted of the ablation work, problems, which will now be reviewed as follows. Solution of the ablation problems dates back to Landau who first used the coordinate transformation to change the variable a fixed domain. An the the of of for [38] domain encountered in the ablation problem to 13 exact, medium in his medium quasi-steady state solution imposed paper. after with the constc Also studied dropping for the case of heat flux condition is the ablation a semi-infinite has been given for the semi-infinite assumption. Using Laplace transformation, a time-variant the ablation Landau was able to derive heat flux imposed solution was only on the fixed boundary. solution for However, heat flux condition. include Specifically, (i) Stefan number two limiting approaching cases zero were examined (negligible capacity) and (ii) Stefan number approaching infinity (negligible latent un i ty, heat). Landau For the employed case of the Stefan a finite difference number method of the order to solve ablation problem with results method presented graphically. has also been used by Rogerson and Chayt to find the exact melt-through time for ablation of a slab imposed a constant heat flux condition on one side and an insulated condition conduction thermal pr on the equation operties other side. 1 and showed th of the material Rogerson e results engaged integrated to be independent the heat of the change. The integral approach has been used in the solution of the phase-change residuals, problems the method [9,40,41]. was first Essentially introduced a method by of weighted von Karman Pohlhausen in the approximate solution of boundary layer equations. The heat integral reasonably a detailed accurate analysis approach results. is simple However, to use; it also provides handicapped Specifically, the form of the profile the steady-state the pre-melt derived for the constant that heat Landau's of with the [39] in phase and of the temperature it is somewhat field. in the is limited by of the temperature accuracy 14 initially chosen for analysis. Also the approximation can not be systematically improved for accuracy. The moment method proposed by Zien [41] carries promise improvement method over the classical has been used to solve time-dependent heat-flux L heat integral one-dimensional conditions. method. ablation Both pre-melt The moment imposed with and ablation and the integral multiplying profile was Again, of the original the integrand s chosen by the heat balance heat equation powers and substituted integral was carried of temperature. A into the integrated was used, out after temperature version of the heat conduction equation. The heat balance integral based the approximate for the boundary temperature heat flux. profile Although was then used the method as the expression appears to work for the general case of the time-variant heat flux condition, it is expected that the choice of an approximate temperature profile works for the nonmonotonic heat flux may not be applicable to the case of more general time-variant conditions. The ablation problems can also be solved by a variational approach. Lagrangian Biot and Agrawal thermodynamics to applied the solut- the variational ion of ablation analysis problems variable thermal properties [42]. They considered one-dimensional heat transfer in a semi-infinite cylinder imposed with a constant heat flux condition. solved. In the procedure, Both pre-ablation the governing and ablation equations were stages were transformed and the Lagrangian heat-flow equation was derived which provided relationship between the surface velocity and the heat penetration relation between these of solutions were derived. on that and with was also two quantities depth. Another 15 found solved by using the simultaneously energy for equation. These equations the heat penetration depth were then and the velocity of the ablated surface. The variational method has been applied for approximate analysis of a and radiative blation of heating a semi-infinite [43,44]. solid Solutions subject were obtair to convective ied in closed form for both the pre-melt and melt-removal heating regimes. these studies, a cubic temperature profile was taken. The surface temperature, melt removal the thermal were treated penetration depth, as unknown and and the depth of the were determined as functions of time. The ablation Blackwell [45] has solve an ablation a moving coordinate problem employed problem has been solved the exponential in one-dimension. system which numerically. differencing He proposed was attached Recently method to the use of to the ablated surface. T approach, th convection oi all elements hen, by invoking ie element n f the moving and control with the exception The node point the matrices grid, use of could and energy volumes were of the last ablating at the moving interface a finite be defined storage. moving element was fixed control-volume for conduction, In this method, at a uniform and control in space. velocity volume. As such, the last ablating boundary. The method e 1 ement had one moving has been applied to boundary and the solution of one fixed a steady- state ablation problem for which an exact solution was available. Comparisons conduction were also made with central terms and upwind differencing differencing for for the convective the terms; appear to be better In n umer ical1l1y. the exponential schemes 16 Ablation problems have also been solved by a finite element method with deforming spatial grids [46]. In this method, classical finite-element the continuous deformation the ablated surface. Th equations ar of the grid e transformed for lis is done at t a precise I he expense to account for localization of of chores of construction of an additional convective matrix. Recently, as laser 1 ablation beam. has been applied Masters [47] to model used the finite intense heating difference method to solve the ablation of a one dimensional slab imposed intense uniform heat flux and analyzed the effect of melt on the temperature solutions w distribution ere found foi the temperature studied the abl history ation during the heat pulse. r the velocity duri ng due to ablation. of the surface Abak i an e laser a continuous-wav The steady-state recession and Modest and [48] beam irradiating on a moving semi-infinite and semitransparent solid. Using integral method, they were able to derive a set of nonlinear partial differential equations wh ich were solved numerically for the groove depth and shape due to ablation. They also considered the ablation a moving slab caused by irradiation from continuous-wave pulsed laser distribution processing a vaporization beams [49]. and derived The laser Ls investigated occurred at a solution for the temperature has also been used in material by Dabby ar the surface; ad Paek [50]. however, In below their work, the surface, the material was heated by absorption of the laser radiation, which might reach Explosion a temperature may thus take place, higher which than that provides a for means vaporization. for material in the drilling the such with an of an and removal process. 17 The review heat conduction given above problems. is strictly In these problems, of regular geometry, the system governing specified, equation, initial and the problems and boundary are solved mainly conditions are fully for the determination of the temperature field. In the inverse problems, however, roles of known and unknown quantities are exchanged. There been abundant of the regular the solution studies doc problems; of inverse cumented in the literature not much work, problems, an however, has d for those for the solution been devoted to that have been solved and documented in the literature, nearly all of them are for heat conduction study for inverse without problems phase change [51-60]. change. There is a lack of Most inverse heat conduction problems deal with a situation an extra and this temperature temperature is available is used together at one point with others in the domain, to find the condition imposed on the boundary. Stolz [51] first solved problem formulated conduction principle. condition found a small numerically. as though problem An integral and it time step [52] was able to improve In his work, were a direct he equation by the problem. was able to was derived numerical problem a linear use the superposition The method inversion. used were too small. an accurate by using solution. a procedure involved minimization of the sum of the squared difference between the actual and the calculated temperatures at the location [53] developed where the temperature data for the solution the have with phase where such was solved, inverse Since was heat was solved to be inefficient for the unknown surface if the time steps must still be used for this method was Yet, Beck that were given. Burggraf an exact 18 series solution capacitance temperature point. data to a one-dimensional approximation serving and heat flux histories Approximate were used for input inverse as results were found in solution. the problem leading were provided if discrete Beck [54] with the lump term. The at an interior or experimental moved one step further by solving an inverse problem in which the material properties were treated as functions of temperature; problem element solved solver becomes nonlinear. has also been reported A two-dimensional inverse in the literature finite [55]. can be used for the solution of the heat flux imposed on the surface of nuclear measurements fuel rod with the use of the interior temperature for input. a handful of studies are found for the solution of inverse heat transfer with phase change. Macqueene et al. [61] proposed inverse welding finite element process. method In their to determine method, the efficiency the latent of an arc heat of fusion is taken into consideration by the variation of the elemental specific heat when the point. average Conduction temperature of the element in both the solid and reaches liquid the melting regions was accounted for. Katz and Rubinsky [62] proposed a front-tracking finite-element method for the solution of one-dimensional inverse Stefan problems. His method was applied for the determination the position of the solid-liquid interface and the transient temperature distribution in the solid region during stationary welding. An inverse method was also used by Landram [63] for the analysis of the thus the It Only an of arc interface pos i tion and the energy transport 19 mechanisms during to be important welding. during The vaporization energy the motion loss of the solid-liquid was found interface. The interface shape, being close to hemispherical, gives clear indication transfer. of a nearly An analytic one-dimensional solution to inverse radial Stefan symmetry problems for the heat in Cartesian and spherical geometries was provided by Rubinsky and Shitzer In their analysis, the inverse Stefan problem was characterized two boundary conditions at the moving front. This is in sharp contrast to the technique developed in this dissertation Chapter 3) which needs only one condition at the interface to solve for the boundary at the phase cl derived solved guessed problem. condition. hange by integrating by the method solution Series In their temperature. the governing of analytic was taken solution w; work, the medium An integral equation, iteration to be the long-ti as then developed by I was initially equation which in whi me sol induct was then was, in turn, ch the first lution to the tion following a number of iterations. A boundary element analysis with constant elements has been developed inverse analysis by Zabaras solidification developed by et al. problem for the solution [1. Beck [52-54] They and Burggraf of a one-dimensional used the sensitivity [53] for inverse heat transfer able to solve region given solution. two separate and the other by Rubinsky Using inverse in the liquid and Shi tzer, an integral Stefan region. two conditions formulation, problems, they were one in the solid Thus similar to the ones must be provided at front, also an inefficient method. E64]. by (see the freezing 20 It should be mentioned that a source-and-sink method recently been developed in the heat transfer literature, which particularly suited for the solution of regular and inverse Stefan problems. to apply In a series this method of papers, Hsieh to the solution and his associates of regular were able and inverse one- and two-phase without melting subcooling and solidification and superheating problems and imposed for medium with and with constant and monotonic temperature has been applied and heat flux conditions to the solution of phase change [2,3]. imposed The method with cyclic conditions developed [4,5,65]. Two inverse with the problem solution formulated with techniques have also been a source-and-sink method as shown in reference to the boundary 3 element The method method has shown as reported by to be closely Hsieh et al. related [66,67]. a further extension of these has is study is works. The present CHAPTER III SOLUTION OF INVERSE STEFAN PROBLEMS BY A SOURCE-AND-SINK METHOD This chapter presents the development of a source-and-sink method to solve inverse Stefan problems. In these problems, conditions boundary. the phases are specified Typically, is given at the moving the temporal location and this location rather than of the interface is used together the fixed between with others to determine the boundary temperature or heat flux that is required to provide for this interface motion. In what will be given designed follows first. for the solution in this chapter, It is followed the motivation by of the inverse a general MStefan for this study analysis problems that i s in three dimensions. dimensional examples extensions This analysis inverse are problem provided of the method and is then used in the solution examples. discussed are included Finally, in results details to conclude and of one- for these possible this chapter. 3.1 Motivation Heat diffusion in a medium with constant properties governed by the partial differential equation V 2T(ft) + k - aT(ft) Ot the is YER t>0 21 (3.1) 22 where all notations have their usual meaning. For this medium, conditions imposed on the boundary are usually one of the following types T(f,,t) = Fi(Yi,t), ri E Bi T(fi,t) = Gi(iBt) Oni ki OT(f,t) o9ni T(F,t) - 1 Hi(fit), (3.4) which represent conditions, the respectively. familiar I In (3.3) Dirichlet, and (3.4), Neumann, denotes and Robin an outward drawn normal. Then, with the additional initial condition given T(FO) = Ti(i) the temperature solution can be expressed by means of Green's function as [68] +0 f G(fit I r',O)Ti(r')dV' R/ J f G(ft I r',r)u"'(r',r)dV'dr + OR' Here, the braced conditions given term is used to account earlier. Their expressions for the three are 1 i sted boundary in Table 3.1. A distinct feature is found in the Green's function method above--the effects of the initial condition, heat generation conditions are embodied the (3.2) (3.3) as (3-5) T(ft) (3.6) 1 (or Yi E Bi destruct ion), respectively in and boundary 23 Table 3.1 Expressions conditions to account for effects of boundary U BOUNDARY CONDITION 4 t)F1 j t) t) -G1 ï¿½l, t) i t)-=Hi (T~i t) I { EXPRESSION 0 a G (r, osl C tIT'~, V 0/ t fG(T tI , os/ ,7) dS;dtc dn. ï¿½T (T/i, t) h adnrI) + dn1 dS1 dt ,k i w (7-4j I Fj (74i )n f D (-I jf ' ki '... 24 the first, second, and third terms on the right of equation (3.6). Then, in the case of regular problems, once these conditions fully specified, t effort, the Green's ;he temperature function can be easily can be obtained found. by using In this the concept of point charges [68] or by solving auxiliary problems [69]. The format the solution of (3.6) of the inverse turns out to be particularly problems. For the inverse suited problems, left hand side of this equation can be used to represent the extra temperature points. Th: information, that is provided en this temperature such as the init either at the boundary can be used to determine ial condition, heat ger or at interior the missing =ration, or boundary conditions. In these efforts, the missing quantities be expanded substituted either in into their a polynomial respective or an infinite integrals in (3.6). series, which is The resulting equation series initial is then solved to complete condition numerically the solution. for the coefficients This method and the heat generation. clearly works in the for the for the boundary condition, assume a particular since it is unknown 'type' of condition a priori, one must first on the boundary. For convenience, one could use equation (1.2) or (1.3) for this condition. Green'I s function is then found with this assumed exchanged condition. Notice in by Hsieh condition [70]. can be That i s a Dirichlet condition can be accomplished by means Neumann condition and vice versa. The search for the condition thus not restricted yet, an incremental by the type of the conditions assumed. Better approach can also are for the can However, that is imposed to if necessary as shown say, that the boundary and Shang of is -- - - - - - - -- 7 lie be developed to track solution 25 the boundary conditions accurately as will be shown later. concept problems above will now be applied in which the interface for the solution positions of inverse are given, Stefan and these positions are used to find the boundary conditions that must be imposed to cause the interface motions. 3.2 General Analysis For the sake of illustration in what follows the Stefan problems consist of two stages: a pre-melt stage, when heat is added to the surface of a subcooled med i um to raise its temperature to the phase-change temperature; and a melting stage, when the medium changes phase and the melting starts at the surface interface moves inward with time. It is assumed that the properties for different has a distinct medium. will be phases are constant melting temperature Convection developed is negligible. for the and of equal value. ; that is, no mushy For generality, solution of both The medium zone in the the analysis melting and solidification problems. The analysis can also be extended for the solution unequal moment, medium of Stefan phase problems properties the simplified shown in Figure in multiple phases and in as will be discussed problems 3.1. will be solved The formulation a medium later. with For the by considering of these the problems follows below. The and the Figure 3.1 System analyzed Pre-melt Stage Governing equation-- FE (L US) V 2To(f,t) = 1T0(rt) a t , (3.7) to> t> Initial To(F,0) = Ti(f) Melting condition-- (3.8) Stage Liquid Region: Governing V 2TL(f,t) = equation-- 1 aTL(ft) at (3.9) t>to Solid Region: Governing V 2Ts(f,t) equation-- -1 OTs(ft) at 28 t>to (3.10) Initial condition-- Ts(fto) = To(ft0) Interface TL(Fft) =TV Conditions: Ts(rf t) OTs(rf,t) On OTL(rf,t) On Svn(t) vn(t) =V.n Here F 1 interface denotes motion. the interface position For the inverse and Stefan v,, the history problems of the of interest in this study, this history is used for the determination of the missing boundary conditions. The problems as posed can be solved by use of the Green's function method described in the preceding section. However, direct use of this method would require (3.6) to be applied to two separate regions, liquid and solid, and the solution so obtained not be as efficient thus used [2,4,15]. as one desires. In this method, A source-and-sink the melting interface method is is taken a moving be a moving conventional represent heat-sink heat-source methods in the temperatures front and front. which a freezing Then, different in different interface in sharp equations regions, only i s taken contrast to to are used to one equation it is in the solid 29 (3.11) (3.12) (3.13) (3.14) to be may will be derived. Whether or liquid region is 30 determined equation. simplified by the position The solution with this method that is assigned of the inver 1. Following in the temperature problems this approach can then be the melting stage is solved by considering an equivalent problem as follows: Governing equation for the equivalent problem: FE (LUS) V 2T(ft) PL vn(t)b(f - f) 1 0T(ft) (9t Initial condition for the equivalent problem: T(fit0) = T0(ft0) Interface conditions for the equivalent problem: T(?1,t) = Tm, vn(t) = Vii (3.17a, b) where denotes a Dirac delta function. The signs preceding this function are used for freezing (+) and melting (-). can be shown readily that (3.15) reduces to (3.9) (3.10). Furthermore, by integrating (3.15) across the interface rf-e to Ff+e and forcing e to be zero in a limiting process, equation pill (3.15) reduces box at the interface to (3.13). This in Figure 3.1. can be proved by using Other equivalences are apparent. t>to (3.15) (3.16) It and from the for b(lz- Ff) this stage 31 The equivalent which problem the heat generation m can be solved term is changed by referring to (3.6) to the interface in motion term as t L SG(ft t0LuS LI LuS - ff)dV dr + (3.18) where the plus sign is used for freezing and minus sign for melting. Finally, equation the boundary conditions to the interface position, can be found rf, and T(FI by setting t) F in this to the melting temperature, Tm 9 as L G(Ff,t I f',t0)Ti(f')dV' t4' toLuS G(?ft I F',r)v - r1)dV'dr + .j{ The missing implicitly. boundary conditions In this effort, can then be found the time when melting by solving starts (to) them can be determined be taken generation by solving the pre-melt to be the special term is zero. Solution problem, whose case of (3.6) in this stage solution in which can again the heat is thus elementary. 3.3 Example Problems The analysis above is now used to solve example problems which are in semi-infinite domain in which the interface motion is given. For the sake of generality, the analysis will be developed determining either the Dirichlet condition or the Neumann condition the interface motion may Tm =I LuS (3.19) for G(Fjt I f',t0)Ti(f')dV' T(fit)= 1 71(7)6(f/ that is imposed on the boundary at x=O, and 32 be the result of either melting or solidification in a one- and two- phase medium. Also for illustration purposes, the heat flow is dimensional, the initial temperature being un i form. Extensions more dimensions are provided in Appendix For the problem given, Green'Is function can be found to be G(x,t I x',r)= 1 [exp( 2 i-ra(t - -r) (X - x)2 4a(t -- )) where the plus and minus signs are to be used when the flux and temperature respectively. condition is assumed to appear at can then be obtained the by using boundary, (3.18), which is recast in a general format To(x,t) Tm ^* t-to H(t - to) Ste f dR~r+t0) d- G(xt I R( r+to),r)dr where all temperatures, including Tm , are measured in excess of the initial temperature, To(x,t) t E(s) = Here F(s) x F(s) 2a t-s G(s) and G(s) denote the assumed temperature and heat flux c respectively. one A to (3.20) The temperature as T(x,t) Tm (3.21) and x2 exp[ - 4( ) ]ds in which (3.22) (3.23a,b) condition Also H(t-tï¿½) - 0 for t>to t CTM Ste -L where Ste is known as the Stefan number. With the use of the circumflexed Heaviside function given by (3.24), equation (3.21) holds for all time and for both one- and two-phase problems. Equation (3.21) can be used to determine the unknown boundary condition by invoking T0(R(t) a=1- t) use of the condition t H(t - to) - Ste -t 0 dR(r+t) G(R(t),t d7- t the interface to IR(r+to),r)dr where the plus and minus signs on the left hand side are to be used for freezing and melting, respectively. 3.4 Numerical Solution A local In this effort, can be used time to solve range (3.26) is divided numerically. into small increments see Figure in which 3.2. each increment, the interface position can be taken integral is taken to be linear; out of the integral for written as a summation as 33 (3.24) (3.25) as (3.26) linearization the entire Then the dR/dt and the convolution Figure 3.2 Linearization position of solid-liquid curves for interface the numerical solution 35 R(t) R(tN) R(t N-1) R(tn) R(tn-I) R(t2> R(t 1) R (t) R (tn -1) -I+n(A-1t (A 0 (AR)n (A t)I to* tI t2 tn 1 tn tN1 tN t I I I I I I I 0 -to T -to 36 TO(R(tN)'tN) =l H(tN- to, Tm Ste m n=l dR(t-) dt t to G(R(tN),tN - to I R(-+to),r)dr] tn_1 - to Here the signs are to be selected following the statement (3.27) below (3.26). evaluated For the invers with the input problems, the right of the interface motion hand side data. As can be for the left hand side, if the boundary conditions are represented power series, the number then the number of the terms that of terms are taken on this side must be equal in the series. to In practice, equation (3.27) can be written by means of matrix elements melting problem for the series method N Sbmnan = Cm n=l N + dmnen n=l where m=N; 10 tN > to; for N=I ,2, .. m bmn = R(tm) 2a sn-1 (tin -S3/ exp[ R(t7n)2 4a(tin - s)"ds for in> n; N F(s) = L ans"-1 n-1 tm 1f sn-i k J (tm- S)1/2 0 R(tm)2 exp[ - 4a(:tm-s )]ds for m n; N G(s) = Lann-1 n=l (3.29a,b,c) (3.30) cm - iaTm by as for (3.28) and m tn - to Lif G(R(tm),tm - toR(r + t0),r)dr for m> n tn_1 - to dR(t-) dt Notice closed series that the integrals form if the missing as shown coefficient in (3.29) matrices Band D in these conditions or in terms are of of a lower can be performed in the series. a Fourier triangular in power The structure, characteristic permitting the solution of (3.28) to be carried simply by using forward substitution, a simple numerical procedure. The boundary conditions can also be determined by use of incremental tN tN-1 approach. In this E(s) exp (tN- S)1/2 method , (3.27) is recast ]ds = 4 !Tm{ï¿½ 1 H(tN - to)N Ste n-1 ti - to [dR(t-) I G(R(tN),tN - t0 I R(r+to),r)dr]} - tn_1 - to N-1 {ï¿½zF n=1 where to 1 tn f tn__ E(s) (tN _ s)1/2 exp[- ]ds} (3.33) to zero all the time in the integral dmn =1 (3.3 1a,b) (3.32) equations are expressed out an as 0 for on the left and s set 38 the last integral on the right. In practice, equation (3.33) can be recast F(tN) G(tN) PNG N Tm + EdNnen n=l N-1 -E n=l F(tn)Pn,F G(tn)Pn,G (3.34 ab) where tN> to ;N=I ,2, ï¿½ , dNn can be obtained by referring (3.31); tn R(tN)f 2a f tn - tn I tn-1 1 (tN s 1 (tN-S)/ exp[ - R(tN)2 4a(tN -s) R(tN)2 exp[ 4a(tN-s) (3.35 a,b) In (3.34), the summation vanishes when N-1=0. Then starting so on, the conditions can be evaluated incrementally. this effort, are used as the F and G values the input found in the computati from the previous time steps on of the right hand side of (3.34), which immediately gives the F and G values at the succeeding new time step. The algorithms This continues until can be developed the desired for a rapid time solution i s reached. of the conditions. Computer programs developed for the solution of inverse Stefan problems given in this chapter are provided in Appendices through E. as and to PnF PnG ]ds N=I, and from In 39 3.5 Critique of the Method Using Green's functions, the present method is limited to the solution of problems whose Green'Is function can be obtained analytically. This excludes those problems whose boundary conditions nonlinearity Kirchhoff tr and governing equations can not be resolved 'ansformation). are by using nonlinear and transformation use of Green's function, projects an impression of being related to the boundary element method that [64]. Yet method, th the liquid t has been undergoing t, they are e boundary and solid . functionally element regions, through development different. equations whereas are written s in the present in recen" ;eparatf metho( t years element ely for d, only one equation solved (3.27) in the prese is derived. nt method The number of is thus reduced the equations by half, which particularly true in the solution by the incremental approach described earlier. More important, in the boundary element method, the problems cannot be solved without the information for the heat fluxes appear on both sides of the interface [64,1]. Such fluxes, however, are unnecessary in the present work. Instead, equation (3.26) embodies conditions (3.12) and (3.13) in the form of a single integrodifferential equation. There is no need for the satisfaction of the flux conditions; as a result, the present method effective because it requires less information in the next section, for input. the present method also accurate. In fact, is exact. such accuracy is not unexpected because The only approximation Yet with the such (e.g., it In boundary to be is that As will be shown is more is (3.26) in the analysis is equation 40 the local linearization which has been applied to the interface motion and the boundary in the present C constant elements. analysis, condition; they The analysis see (3.33) and (3.35). have been approximated can thus be improved for In fact, by using accuracy by using elements higher order as in references elements [59,71]. such as linear and quadratic 3.6 Results and Discussion For the numerical experiments performed in this study, interface motion data are taken from reference 5. Aluminum will be used for tests; its properties are given in Table 3.2. The inverse solution techniques developed in this chapter are used to solve eight example problems of which four having exact solutions. retrieved the exact interfaces interface conditions conditions solutions for these for error. move as a function motion that Stefan-Neumann data examples In square of are used to retrieve appear on the surface problem, and Table can thus be compared these four examples, root of time, the constant t( This problem 3.2 provides with the and the emperature i s known a summary as the of the conditions tested The first in the example examples. deals with a general two-phase Stefan- Neumann problem; the medium is initially subcoo led to 300 K, which is lower than the melting temperature (see Table 3.3). For this example, the temperature temperature is thus assun at the boundary ned to appear a is unknown; t the boundary an unknown S, and the interface motion data are used to find this temperature. two sets of results are the The The results are listed in Table 3.4, where 41 Table 3.2 Properties of aluminum Melting temperature, T. (K) 932 Vaporization temperature, T, (K) 2,543 Heat of fusion, Lf (J/kg) 389,600 Heat of vaporization, L, (J/kg) 9,462,000 Thermal conductivity, k (W/m K) 200 Density, p (kg/rn3) 2,710 Specific heat, c (J/kg K) 1200 42 Table 3.3 Conditions tested in eight examples PROBLEM INPUT TRUE CONDITION TYPE OF BOUNDARY DESCRIPTION DATA F(t):K; G(t):W/m2 CONDITION ASSUMED SOLVED 1 two-phase R(t)-4t F(t)=1000 tempeatur temperature Tj=300K TpS300K 4 one-phase R(t)-4t F(t)= 1000 heat flux temp ure Ti=932K=T 5 one-phase reference F(t)= 1000+5t temperature temperature Ti=932K=T, [5] 6 two-phase reference G(t)=6.388034x 106 heat flux heat flux T300K 7 two-phase reference F(t)f932+40t temperature temperature Tr=300K 8 one-phase reference G(t)=3xl06+5xl0'ta heat flux heat flux T1932K=Tm [5] LIII I I 43 Table 3.4 Comparison between true and retrieved conditions for first Stefan-Neumann example solved for boundary temperature ________Boundary Conditions, T (K) Imposed True Condition Retrieved_____ Series Method N=3 ii' i 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1003.22 1002.97 1002.72 1002A8 1002.25 1002.03 1001.82 1001.61 1001.42 100123 1001.06 1000.89 1000.73 1000.58 1000A3 1000.30 1000.17 1000.06 999.95 999.85 II I aI N=4 1003.17 1002.87 1002.58 1002.30 1002.04 1001.79 100136 100134 1001.14 1000.95 1000.78 1000.63 1000.50 1000.38 100028 100020 1000.14 1000.09 1000.07 1000.07 Incremental Method 1003.49 1000.41 100020 1000.11 1000.03 1000.05 1000.04 1000.02 1000.68 999.88 999.97 999.99 999.99 999.99 999.99 999.98 999.98 999.98 999.90 999.94 F (S) =ï¿½ ansn-1 nul For N=3, a1=1003.490405 a = -0.268732 a%= 4.350448xl03 For N-4, a1=1003.490405 -0.322478 6.264645x10"3 6.552981x10"s Time (sec) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 44 given--one incremental for the series solution In the series method. solution solution method, method a power and the other series of degree for the N-1 is used to represent the temperature (see (3.29b)), and two values are tested. I time intervals In this method, the interface are used for input. Thus, position da for example, ta at equal for N=3, R data at times determine equal the series, to 0, which 6, 12, and 18 seconds is, in turn, are used to used to generate the temperature v the retrieved alues listed 1 temperatures for 20 time steps at the boundary in the table. by both series Comparing with the true condition (exact temperature) listed to the left indicates are in good agreement. these series are listed In this case, at the bottom the coefficients of the table. found for From the temperature values tabulated, the results for the low power series appear trend to be persists 10 (results temperatures as good as those of the high even with the test of not shown). a higher This gives have been converged. power power series, series the indication As for the incremental and such of degree that the solution method, the results are also good; errors are of the order of 10-3 at large time. For the incremental method, the boundary conditions are foun given. d at exactly The time step the same times for the solution en the interface positions are is thus identical to that for the interface data input. Also notice that the convergence stability difference solution o that are normally encountered methods )f integral are nonexistent equations. Also in in the conventional the present as in the finite incremental case of the series the incremental solution of N they and are generally results solution method, 45 better at large time than small time, which will be further discussed later. In the example above, a temperature condition is imposed, same 'type' of condition is assumed in the process of the inverse solution. Since the type of the condition that is imposed the boundary is unknown a priori, it may well be the heat flux condition that one assumes, and the question to be addressed now is whether it is still possible to retrieve the temperature condition via the assumed flux condition. Use will now be made of equation (3.21) to determine flux condition is fou this temperature co ind, and the results ndition once the boundary are listed in Table As shown good. but in the table, the series the incremental solution results solution results are not as accurate are still as those i sted in Table 3.4. This is certainly a result of the errors being accumulated evaluation flux. first in the evaluation of the temperature While this example using serves of the heat flux next in the the previously well determined to illustrate heat that the conditions condition are still may exchangeable, lead to large such errors, a two-step particularly solution of the in the series solution according approach compared method, and should to experience, of solving fl thus be avoided the computer Lux then with that in the separate, in practice. time saved temperature one-step is In fact, in this two-step insignificant approach as of direct evaluation of the flux and temperature. A slight modification is made in the next two examples: time the medium is not subcooled, the initial temperature temperature of the medium being and 4 the and on 3.5. this (see examples 3 equal to the melting 46 Table 3.5 Comparison between true and retrieved conditions for second Stefan-Neumann example solved for boundary flux then temperature Retrieved Series Method N-4 N=5 i i i i 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 585.89 715.58 799.99 862.06 909.11 945.07 972.34 992.57 1006.96 1016.46 1021.85 1023.76 1022.76 1019.34 1013.95 1007.01 998.89 989.05 980.54 971.00 620.85 758.46 844.90 905.50 948.70 979.11 999.67 1012.44 1019.03 1020.76 1018.76 1014.02 1007.48 999.99 992.40 985.51 980.13 977.04 977.04 980.93 __ _ __II Na Incremental Method 1023.06 1007.71 1003.60 1000.88 1000.04 999.68 999.52 999.44 998.11 998.30 998.37 998.44 998.51 998.59 998.66 999.45 999.16 999.11 999.07 997.15 G(S) =ï¿½ ansn-1 nsl For N--4, For N=5, a,=8001187.956002 a2= -387201.858263 a3= 1581.702473 a,= 148.390536 al=8945600.085189 a2= -541131.047558 2763.120508 324.034498 a,= 3.817008 [________Boundary Conditions, T (K) Time (sec) Imposed True Condition 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 47 in Table one-phase 3.3). problem. The Stefan-Neumann It will be prob1 ems tested solved that thus become the results unaffected series by the change of tests to the one-phase are made and the results problem. A are listed gain, the in Tables and 3.7. Again, The two- step reinforcing the one-step solution solution results the recommendation results are poorer made earlier are good (Table (Table 3.7), in the testing 3.6). further of the two-phase problems. Having problems, at satisfactorily ttention is completed now directed testing of the Stefan-Neumann to the solution of inverse Stefan problems whose interface motions must be met by imposing time-variant temperature and flux conditions. There are no exact solutions for these problems, and the interface motion taken from Choi [5] who have solved the regular (forward) version the problems with great accuracy. The interface position then used to retrieve listed results in Tables the boundary 3.8 through for the direct retrieval 3.11. conditions Tables of the linear and the results 3.8 and 3.9 give temperature and flux conditions F(t) = 1000 + 5t G(t) - 6.388034 x 106 + 2.82752 x 105t while Tables 3.10 and 3.11 give the results for the direct retrieval of the linear temperature are same 3.6 the data are of data are are the heat (3.36) (3.37) heat flux condi ti ons and quadratic 48 Table 3.6 Comparison between true and retr for third Stefan-Neumann example temperature ieved conditions solved for boun [I ___ Boundary Conditions, T (K) unposed True Condition Retrieved .. ... I ' Series Method N=3 N=4 Incremental Method a a a a 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1033.83 1030.31 1026.931023.69 1020.59 1017.63 1014.81 1012.13 1009.59 1007.19 1004.93 1002.81 1000.82 998.98 99728 995.71 994.29 993.00 991.86 990.85 1033.12 1028.96 1025.01 1021.29 1017.80 1014.55 1011.55 1008.81 1006.33 1004.13 1002.21 1000.57 999.23 998.20 997.47 997.07 997.00 997.26 997.87 998.82 1037.49 99822 998.58 998.85 999.05 999.20 999.24 999.38 999.44 999.49 999.53 999.59 999.60 999.65 999.65 999.69 999.69 999.73 999.71 999.78 a a a a N F(s) = asn-1 nul For N=3, a1=1037A96339 a27= -3.728416 a*= 6.982696xl02 For N-4, aj=1037A96324 ai -4.474102 a3= 1.005513x10' ak,= 1.324365x103 dary Time (sec) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 l lNO 49 Table 3.7 Compar i son for fourth flux then t between true and retrieved conditions Stefan-Neumann example solved for bou ;emperature ndary I______BoundaryConditions,_T_(K) Imposed True Condition Retrieved Series Method N-4 ii i i 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 906.32 953.87 980.48 997.67 1008.95 1016.02 1019.89 1021.28 1020.73 1018.66 1015.44 101139 1006.79 1001.89 996.94 992.16 987.77 983.98 980.99 979.00 N=5 917.43 96625 991.90 1006.84 1015.10 1018.67 1018.86 1016.61 1012.66 1007.67 1002.18 996.73 991.80 987.87 985.38 984.80 986.56 991.12 998.92 1010.44 J I I a Incremental Method 1056.39 996.60 1002.26 999.33 999.70 999.45 999.16 999.22 999.09 999.23 999.05 999.32 998.89 999.72 998.27 1001.19 994.04 1014.19 951.97 1162.46 N na sn-1 nul For N=4, For N=5, a,=2805963.114624 a2= -204888.890714 a3. 2268.882336 4a-= 99.892255 a,=3137161.977178 a2= -286340.839757 a3= 3963.574559 a4= 218.130426 a,-- 1.467676 Time (sec) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 G(s) mummi 50 Table 3.8 Comparison between true and retrieved conditions for fifth inverse example problem solved for temperature I____ Boundary Conditions, T (K) Imposed True Condition iRetrieved Series Method N=3 Error N-4 Err Incremental Method I I 'U I I 1005.00 1010.00 1015.00 1020.00 1025.00 1030.00 1035.00 1040.00 1045.00 1050.00 1055.00 1060.00 1065.00 1070.00 1075.00 1080.00 1085.00 1090.00 1095.00 1100.00 1073.05 1064.27 1056.80 1050.64 1045.81 104229 1040.09 1039.21 1039.64 1041.40 1044.47 1048.86 1054.56 1061.59 1069.93 1079.58 1090.56 1102.85 1116.47 1131.39 6.77 5.37 4.11 3.00 2.03 1.19 0.49 -0.07 -0.51 -0.81 -0.99 -1.05 -0.97 -0.78 -0.47 -0.03 0.51 1.17 1.96 2.85 1078.13 1061.16 1048.37 1039.36 1033.71 1031.04 1030.94 1033.01 1036.86 1042.07 1048.26 1055.02 1061.95 1068.65 1074.72 1079.77 1083.38 1085.17 1084.72 1081.65 7.28 5.07 329 1.90 0.85 0.10 -0.39 -0.67 -0.78 -0.75 -0.64 -0.47 -0.29 -0.13 -0.02 -0.02 -0.15 -0.44 -0.94 -1.67 1003.61 1007.83 1012.50 1017.33 1022.27 1027.17 1032.24 1037.26 1042.29 1047.31 1052.35 1057.36 1062.38 1067.38 1072.39 1077.40 1082.37 1087.50 1092.06 1098.80 II Ii ,l II I a a = ansn-I nul For N=3, For N-4, a1=1083.1636095 a2= -10.7640380 a3= 0.6587919 a,=1099.6680063 a2= -23.9563496 2.4866623 a4= -0.0666941 Time (sec) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 Error -0.13 -0.21 -0.24 -0.26 -0.26 -0.27 -0.26 -0.26 -0.25 -0.25 -0.25 -0.24 -0.24 -0.24 -0.24 -0.23 -0.24 -0.22 -0.26 -0.10 F(s) 51 Table 3.9 Comparison between true and retrieved conditions for sixth inverse example problem solved for heat flux Boundary Conditions, q (Wlm) Imposed True Condition Retrieved Series Method N=4 Error N=5 Error Incremental Method 1 i I I I i 7589730.94 7660418.94 7731106.94 7801794.94 7872482.94 7943170.94 8013858.94 8084546.94 8155234.94 8225922.94 8296610.94 8367298.94 8437986.94 8508674.94 8579362.94 8650050.94 8720738.94 8791426.94 8862114.94 8932802.94 8076200.72 8147001.48 8206716.05 8257180.74 8300231.83 8337705.60 8371438.35 8403266.35 8435025.90 8468553.29 8505684.79 8548256.71 8598105.32 8657066.91 8726977.77 8809674.19 8906992.45 9020768.84 9152839.65 9305041.17 6.41 6.35 6.15 5.83 5.43 4.96 4.46 3.94 3.43 2.94 2.52 2.16 1.89 1.74 1.72 1.84 2.13 2.60 3.28 4.16 7558675.50 7694176.58 7815212.59 7920261.67 8008672.53 8080664.49 8137327.48 8180622.00 8213379.15 8239300.64 8262958.75 8289796.37 8326126.99 8379134.69 8456874.13 8568270.59 8723119.93 8932088.60 9206713.65 9559402.73 -0.40 0.44 1.08 1.52 1.73 1.73 1.54 1.19 0.71 0.16 -0.40 -0.93 -1.33 -0.61 -1.42 -0.94 0.02 1.60 3.88 7.01 7184872.69 7371675.43 7493717.81 7596454.24 7689904.72 7775533.53 7859641.15 7940899.93 8020420.49 8098888.64 8176244.61 8252784.49 8329871.72 8402584.80 8480061.68 8553180.08 8628203.50 8701314.04 8775467.98 8847997.32 a a a a a i m For N-4, G (S) =f anbs-' n s For N=5, a,=3397040.84 a2=2247925.15 %=-353114.747 a4= 19587.0623 a,=8734725.30 a2=-3092029.40 a3= 1293821.95 a4=-188030.644 as= 9286.34410 Time (sec) 4.25 4.50 4.75 5.00 5.25 5.50 5.75 6.00 6.25 6.50 6.75 7.00 7.25 7.50 7.75 8.00 8.25 8.50 8.75 9.00 Error -5.33 -3.77 -3.07 -2.63 -2.31 -2.11 -1.92 -1.77 -1.65 -1.54 -1.45 -1.36 -1.28 -1.25 -1.15 -1.12 -1.06 -1.02 -0.98 -0.95 52 Table 3.10 Comparison between true and retrieved conditions for seventh inverse example problem solved for temperature [IBoundary Conditions, T (K) Imposed True Condition Retrieved Series Method N=3 Error N-4 Error Incremental Method a i i i a i U 972.00 1012.00 1052.00 1092.00 1132.00 1172.00 1212.00 1252.00 1292.00 1332.00 1372.00 1412.00 1452.00 1492.00 1532.00 1572.00 1612.00 1652.00 1692.00 1732.00 976.35 1014.88 1053.63 1092.58 1131.74 1171.11 1210.69 1250.48 1290.48 1330.70 1371.12 1411.75 1452.59 1493.64 1534.91 1576.38 1618.06 1659.95 1702.06 1744.37 0.44 028 0.15 0.05 -0.02 -0.07 -0.10 -0.12 -0.11 -0.09 -0.06 -0.01 0.04 0.11 0.19 0.27 0.37 0.48 0.59 0.71 1020.27 1044.41 1071.70 1101.84 1134.51 1169.40 120612 1244.64 1284.37 1325.09 1366.49 1408.27 1450.12 1491.73 1532.80 1573.01 1612.05 1649.63 1685.42 1719.12 4.97 3.20 1.87 0.90 0.22 -022 -0A7 -0.58 -0.59 -0.51 -0A0 -0.26 -0.12 -0.01 0.05 0.06 0.003 -0.14 -0.39 -0.74 969.93 1006.51 1042.85 1079.43 1116.09 1153.87 1192.15 1230.89 1270.11 1310.73 1348.84 1389.93 1429.95 1470.27 1510.54 1550.77 1590.93 1631.04 1671.39 1710.13 a a a I I a * F(s) =n ansnnal For N=3, a,=938.0379452 a2= 38.2157864 a3= 0.1050509 For N--4, a,=999.5873340 a27 18.8517864 a3= 1.8834608 a4- -0.0513597 Error Time (sec) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 -0.21 -0.54 -0.87 -1.15 -1.40 -1.55 -1.64 -1.69 -1.69 -1.60 -1.69 -1.56 -1.52 -1.46 -1.40 -135 -1.31 -1.27 -1.22 -1.26 WM.j 53 Table 3.11 Comparison between true and retrieved conditions for eight inverse example problem solved for heat flux [______ Boundary Conditions, q (W/m2) Imposed True Condition Retrieved Series Method N-4 Error N=5 Error Incremental Method ï¿½ I I I I I 3003125.00 3012500.00 3028125.00 3050000.00 3078125.00 3112500.00 3153125.00 3200000.00 3253125.00 3312500.00 3378125.00 3450000.00 3528125.00 3612500.00 3703125.00 3800000.00 3903125.00 4012500.00 4128125.00 4250000.00 3007778.76 3021875.02 3040169.58 3062934.56 3090442.06 3122964.20 3160773.08 3204140.84 3253339.57 3308641.38 3370318.40 3438642.74 3513886.50 3596321.81 3686220.76 3783855.48 3889498.09 4003420.68 4125895.38 4257194.29 0.15 0.31 0.40 0.42 0.40 0.34 0.24 0.13 0.01 -0.12 -0.23 -0.33 -0.40 -0.45 -0.46 -0.42 -0.35 -0.22 -0.05 0.16 3033589.06 3035992.94 3043865.79 3057990.24 3078999.09 3107375.31 3143452.06 3187412.63 3239290.49 3298969.30 3366182.86 3440515.14 3521400.30 3608122.64 3699816.64 3795466.96 3893908.40 3993825.95 4093754.75 4192080.13 1.01 0.78 052 0.26 0.02 -0.16 -0.31 -0.39 -0.43 -0.41 -0.35 -0.27 -0.19 -0.12 -0.08 -0.12 -0.23 -0.46 -0.83 -1.36 II 1-. L - I II 2995616.25 3004320.50 3019590.53 3036905.38 3061522.20 3092582.42 3130006.55 3173845.90 3223966.45 3280396.20 3343088.04 3411906.45 3487162.49 3568166.70 3656322.04 3748732.19 3851068.35 3951711.24 4084808.11 4096803.88 For N-4, a,=2997608.684 a2=33190.75084 a3=29232.59062 a4=2902.536744 For N=5, a,=3035721.67 a2=- 162102870 a3= 27733.4034 a4=12343.6674 a,=-1598.21379 Time (sec) 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 425 4.50 4.75 5.00 Error -0.25 -0.27 -0.28 -0.43 -0.54 -0.64 -0.73 -0.82 -0.90 -0.97 -1.04 -1.10 -1.16 -1.23 -126 -1.35 -1.33 -1.51 -1.05 -3.60 nl antsn-1 F(t) = 932 + 40t G(t) = 3 x 106 +5 x 104t2 In these following tables, the computational errors are calculated using definition: Error=1 q where pand q represent retrieved and imposed true conditions, respectively. As described in Table 3.3, those in Tables 3.8 and 3.11 for the medium initially at phase change temperature, while those Tables 3.9 and 3.10 are for the medium initially subcooled at 300 K. The former are thus one phase problems while the latter are two- phase problems. Again both series and incremental solution methods are used for solution Table 3.8, N that is the series as low and their solution as 3, whereas results results are good. converge in seeking For example even with a value the flux condition in of in Table degree 3.9, the series results not shown). converges rapidly In Tables from N=4 to N=5 (higher 3.8 and 3.11 the medium melts as soon as the boundary conditions are imposed, whereas in Tables 3.9 the medium all cases, the starts accuracy to melt at time greater of the results appears than 4 seconds. to be unaffected the time when melting 3.11 follow the takes same trends place. as the The results ones in Table in Tables 3.10 and 3.8 and 3.9. Tests conditions are thus successful. 54 (3.38) (3.39) the (3.40) are in In by for the time variant 55 The inverse solution techniques developed in this study expected to be accurate as mentioned earlier that is close to the end of the previous this work, than small the time. section. accuracy This Yet for the Stefan of the techniques can be attributed problems is better solve at large to the curvature d in time of the interface position curve, which is always large at small time [2]; see Figure will be position 3.2. a slight at small Then, error time. in the numerical associated At large solution of (3.27), with the linearization time, however, the position there of the curve tends fact, as (3.27) r boundary present to be linear; time apidly the linearization error will be progresses, outnumber conditions method can at large the accurate terms under the inaccurate terms to always be evaluated time; see the results diminished. the summation the effect accurately in Tables In in that the with the 3.4 through 3.11. marching This is schemes a distinct reported departure from the trends in the literature in which of other the errors with time. It should also be pointed out that, for the Stefan-Neumann problem chosen for compar i so'n in the present study, there is a singularity of the temperature at zero time. This also contributes to the large discrepancy of the results at small time, which must not be overlooked. It has been firmly established that, in the solution of the Stefan problems, only the Stefan-Neumann problems can be solved exactly. Yet, it is also possible to develop an exact solution an exponential condition imposed on the boundary; such condition, however, has been considered as physically untenable yet, such a condition [7]. are to grow time tend for in the Worse literature - - - - %F gives rise to a 56 constant velocity linearization scheme of the interface, exact a situation in the solution making of the inver the present se problems Perfect results will be obtained, rendering the test of the exponential conditions meaningless. 3.7 Extension and Concluding Remarks The analysis developed for the solution of inverse problems, times for re-melt in this chapter problems can be with multiple and re-freeze readily phases. of the medium extended For such must be closely adapted accounted for, for the development and the analyses [4,65] of the inverse can be readily solution techniques presented in this study. for the solution unequal double of inver for different source of the regular The present -se problems phases and sink fronts problems analysis in which of the medium. must be used in reference 17. can also be extended the properties are For such problems, as given Finally, in the solution it is noted that although problems in one-dimensional, semi-infinite doma in have been solved plane for examples wall) in this work, can also be solved prob 1 ems in finite with the present domains methods. (e.g., For these problems, imposed. requires there are two boundaries Two unknowns the input of and two boundary are thus sought one additional conditions simultaneously, condi ti on are and this in the form of either temperature or heat flux at any interior point close to the boundary where no phase with two phase simultaneously change change takes place. interfaces from both sides On the caused other hand, by separate of the boundaries, motion data for the second for problems heat input the interface this additional [5]. interface will serve as 57 condition. sides In any event, of the interfaces. be applied to the soluti4 no flux information Furthermore, th on of problems i s needed e present method in multiple di at both can also dimensions. Again the inverse solution for these problems can be developed on the basis of the solution problems. of regular versions of these CHAPTER IV SOLUTION OF ABLATION PROBLEMS WITH ONE MOVING BOUNDARY BY A SOURCE-AND-SINK METHOD It is the method boundary purpose for the soluti4 As will be sh of this chapter to on of ablation p own, the essential ,present )roblems feature a source-and-sink with one moving of this method is that the solution will be sought in a fixed domain which does not change with time. In fact, the ablated region is to be treated fictitious boundary boundary. domain will match Then, in which the flux condition that of the imposed with the additional conditior vaporization at the ablated i at the moving temperature given for the moving boundary, the conditions at this boundary overspecified for the fictitious domain. The problem can thus be taken fixed as an inverse boundary problem is sought in the during sense the ablation that the condition period. on the This condition will provide the solution for the necessary conditions at the moving boundary of the problem. 4.1 Solution Methodology a subcooled medium, the problem can be divided into two stages; namely, pre-ablation stage and ablation stage. During pre-ablation order stage, to raise heat is added its temperature to the surface of to the phase-change the medium in temperature. Then with continuous heating, ablation takes place. The ablated 58 as a are for For the 59 surface moves this boundary As shown homogeneous an inward motion. with time and the imposed in Figure 4.1, nd isotropic. it is assumed The thermophysic heat flux follows that the medium cal properties constant. For the ablation problem, the ablated region immediately surface removed phenomenon. upon formation. Radiation It is also assumed is taken that the medium to be a changes phase medium. at a distinct Moreover, I temperature; there is that is, no volumetric no mushy zone in the heat generation in the medium. Under these assumptions, the ablation problem can be formulated Pre-ablatic as follows: )n Stage Governing equation-- FE (AUS) V 2T0(Ft) -1 dT0(ft) 0t t0>t>0 Initial T0(?,O) =Ti(f) Boundary G(t)= -k condition-- (4.2) condition-0T0(f0 ,t) Oni' (4.3) is are is (4.1) Figure 4.1 System analyzed 61 Vn m r ni r oa SI Si SI to S Ablation Stage Solid Region: Governing V 2Ts(ft) - equation-- 1ES ~1 OTs(ft) at ' (4.4) t>to Initial condition-Ts(fto) =T0(ft0) Boundary conditions-- Ts(fs,t) =TV G(t) -- kTs(is,t) G-t)=n +pL V. where all notations have their usual meaning; see nomenclature. rS represents the position for the ablated surface, which is unknown for the ablation TS(F,to) is of ablation. problem; to is the time when ablation the temperature distribution Both to and To(Fto) in the medium are to be found starts; and at the onset from the solution of pre-ablation equation temperature (4.7), stage. which gradient, Notice relates is the only that the boundary the ablation source flux condition, velocity of nonlinearity to in ablation 62 (4.5) (4.6) (4.7) Here the the problem. 63 The ablation problem above is solved by using a source-and- method. In this method, the domain for investigation extended to cover AuS; the ablation front becomes an interior phase- change which interface. the ablated Solution region is thus rendered is fictitious in in the a fixed sense domain, in that it is physically nonexistent; however, the unablated (solid) region real. The conditions imposed at the moving boundary in the original ablation problem phase-change are taken interface. to be those Thus, for the imposed on new problem * the interior in the fixed domain, the interface position and the fictitious condition on the fixed found, boundary they are the unknowns can be used for input to be determined. Once they in the determination are of the temperature totally. A in the solid ,n equivalent region. problem The problem is formulated with is thus solved a source-and-sink method as follows: Equivalent Problem Governing equation-- FE (A US) V 2T(ft) - Initial OT(ft) at pL Fk~,"~ condition-- (4.9) T(fr t0) = T0(ft0) Interface conditions-- sink is is t>to (4.8) -- S) 1 T(?s,t) =TV G(t) k - OT(rs't) On +pLV.i where denotes a Dirac delta function. As has been shown Chapter across 3, equation (4.8) r the interface at FS 'educes t enables o (4.4). the G(t) Also integrating to represent (4.8) the heat flux imposed Equation on the interface (4.8) at the fictitious can be solved by using side. (3.18) in which domains for integration for the first and second terms on the right are changed to AUS. Also minus sign is taken for the Dirac delta function term in which rf is changed to FS G(ft I ?',r)Vii6(?' - Fs)dV'dr + This temperature temperature equation equation (4.10) is forced by setting to satisfy the F to the ablated interface surface position, rs, and T(Fs,t) to the phase-change temperature, G(Ys,t I f',t0)Ti(?')dV' AuS t Lff G(fs,t I ',r)V. 6(?' -Ys)dV'dr + { } toAuS is also differentiated to satisfy 64 (4.10) (4.11) in the T(ft) fS AuS as toAuS (4.12) as (4.13) 6 (-F - -FS) G(fLt I f',t0)Ti(f')dV' TV-= (4.12) Equation as t Lff 65 - ______I __ - J aG (F?t 1 F117) G(t)= - k ArfStan ï¿½Ti(')dV'-TL S an V.n( ai4 - Ys)dV'dr } +pLvii r=rS where differentiation is to be effected at the solid side of the interface. summation Notice term, that the boundary which, according condition to Table is embodied 3.1, should in the be expressed in the form a G( ,t i /T)OT(rt) f On 0o where the following condition i s taken at the fixed boundary: aT(fot) Oni - G(?ft) k for t and - (f ,t) k for t > to where there g(Fot) represents are two unknowns the imposed in (4.13) fictitious and (4.14): heat flux. g(o ,t) and V; Clearly, there are two equations be found to solve them. from the solution of In this effort, a pre-ablation and TO(-Ft0) problem are to as described previously. 4.2 Illustrative Examples Use is now made of the solution methodology given in the preceding examples section to solve will be provided examples and they in semi-infinite include medium med i um. with Four or without subcooling with the boundary imposed with constant or variable flux conditions. (4.14) dfodr (4.15) (4.16) heat 66 For been given a subcooled in Chapter medi un, 3ï¿½ the pre-ablation The pre-ablation stage temperature solution has can be taken from (3.22) as T0(x,t) = T-if 0 G(s) exp[ - s)1/2xp x2 -4-xs]ds 4a(t-S) Thus the time when ablation starts (to) can be found by solving following equation to 0 implicitly ds Also at the moment when ablation starts, the temperature in the solid region is given by the relation T0(x,t0) =Ti (4.19) exp[- x ]ds 4a(to - S)]d which serves as the initial condition for the ablation stage. For the problem at hand, (4.12) can be used to derive temperature T(x,t) = Ti + to fa t kr t0 exp[ - x ]dr 4a(t-r) exp[- 4a(t - ]dr L Also (4.13) t "--to dRl(7) 1 {exp( - (x - R(0))2)+ exp( (x+R(r))2 4ra(t- r) p- 4oat(t -t)- r) (4.20) is used to derive (4.17) the (4.18) as TV = Ti k 4 2'r R 2(t) exp[ 1) 4a(t - 7 67 t 44t0 exp- R-(t)]d4a(t- r)]d dRl(r) dr 1 4ira(t - r){exp( - (Rl(t) - R,(r))2 4a(t --r) ) + exp( - (Rl(t)+RI(r) 2 4a(t- ,r) Equation (4.14) follows _ Rl(t) G(t)- 2--- R1(t) exp[- 4( ) Rl(t) ]dr + RFf( j r=ta g(r) exp[ (t-r)3/2 ) R2(t) 4a(t-- r) dR(-) d7r (t-)312(R R1(r))exp( - (Rl(t) - R1(.r))2 + (Rl(t)+Rl(r))exp(- (Rl(t)+l(r))2 1 4a(t -r) dRl(t) + pL dt Solution of (4.21) and (4.22) subject to the initial condition, to) = 0, yields the results for two unknown functions Rl(t) g(t). 4.3 Numerical Solution of Ablated Front Equations ( integrodi fferential 4.21) and equations, (4.22) which are coupled nonlinear this effort, a local linearization scheme is employed as described in Section 3.4. For the present problem, the entire time range divided treated into small increments as constants. They in which both dRl(t) dt can thus be taken and g(t) are out of their respective integrals, and the convolution integrals in these are changed to summations G(r) 27r t -=t I as (4.21) t :ft ]d'r 4a(t - r) (4.22) and will be solved numerically. In is TV = Ti+ Lp 4-fr equations as 68 G(7) tn g(tn) J r=tn_1 exp[- R2(tN) ]dr 4a(tN - 7) exp[ - dR1 (tn) dt 1 NJ4wa(tN - 'r) exp( - (RI(tN) - Rl(r))2 4a(tN -7-) ) +exp( _((tNR-Rr ))2)' ' 4a(tg 7 (4.23) and to Rl(tN) G(r) 2" ra JIo( tN -T)3/2 4a(tN-) ]dr N +E g(tn) 11=1 Lp N1 tn dR(tn) dt T=tn-1 1 (tN - 7)3/2 {(RI (tN) - R1(r))exp( - (Rl(tN) -2 4a(tN--) + (Rl(tN)+Rl(r))exp(- (Rl(tN)+Rlr))2L dRl(tN) 4a(tN - 7) )jr + pL dt Notice that the removal of dRl(t) dt and g(t) will cause a slight error in the numerical is a gradual L solution function o of the ablated )f t, and front posi tion. so is g(t) during However, the early of ablation. expected The to be small. error associated In practice, with the linearization such an error can always reduced by taking small Tv = Ti + a N n=l L N n=l ]dr G(tN) exp[- ]dr} (4.24) Ri stage is be +1 k4ir time increments. 69 It is also noted that the. ablated position Rl(t) in the Green'Is f velocity, Function can dRl(tn-1) dt time t should and g(t) be relate as shown thus be discounted in these equations to R(tn-1) by in Figure 3.2. as unknown. will be solved means of the The position Numerically, incrementally interface R1 at any the Rl(t) starting from to until the final time tN is reached. This corresponds time marching developed for scheme in the numerical the solution solution. of the ablation Computer problems given programs in this chapter are compiled in Appendix 4.4 Results and Discussion Four examples are tested and they include semi-infinite med i um with and without constant subcooling and time-variant with the moving conditions (see boundary Table imposed with 4.1 for summary). Notice that, in this table, under the column of problem description, the ablated medium in heat-flux region itially condition, is always at phase-change such counted as one phase. temperature as Example 1, there imposed Thus, with for a constant is one ablated phase. It is thus called one-phase prob 1 em in the table. Example 1 deals with an ablation prob 1 em that can be solved exactly medium As shown ablated with in Carslaw a constant and Jaeger velocity [71, for a semi-infinite , the heat flux imposed on the boundary is given by the relation G-[L + c(Tv -Ti) ]pU where is constant and the velocity is related to the ablation position to a f (4.25) as 70 Table 4.1 Conditions tested in four examples Problem Material Heat Flux Remarks Description Condition Imposed G(t)_(W/m2);t_(S) 1. One-phase Aluminum G(t)=5xl0 Exact solution is available for this Ti=932 K=Tn-=Tv case. 2. Two-phase Aluminum G(t)=2x106 No exact solution for this case; Ti=882 K 3. Two-phase Aluminum G(t)=8xl0O t No exact solution for this case; Tj=882 K 4. Two-phase Aluminum G(t)=2xl0 t2 No exact solution for this case; Ti=882 K 71 Rl(t) t It is expected that, for a medium initially at the phase-change temperature, C is equal to pLU product. It is also noted equation state. given as (4.25) is strictly At that time, valid the temperature for the medium in in the solid a quasi-steady is stabilized u W(X - Ut) (4.27) which time. occurs The first for the medium two examples having are ablated thus chosen for a long period to substantiate of points. For the first example, the ablation is due to melting. med i um ablates surface is initially at the melting as soon as the heat is applied. is shown in Figure 4.2, where temperature so that the surface The position the curve of the ablated appears linear error Using curves the exact shown solution 4.3. (4.25) for comparison Clearly, the results yields the are stable converge satisfactorily, errors being less than one tenth percent slightly for all the time increments larger at smallI time, tested. a result Also the of the errors are singularity associated with the heat flux that is abruptly applied on the surface at time zero. In what follows in the numerical computation of all the examples to be presented in this chapter, the time increment i s taken as 0.5 sec, unless otherwise noted. for the algorithms (4.26) that T(x,t) = Tve these The and in Figure to be of one an excellent test Example 1 provides Figure 4.2 Trend o med i um imposed f i ablated surface position nitially at phase-change with a constant heat-flux for a temperature condition (s) I sewi OL GL, 0 000 CA rm 0 c C,) 0 "3 0 rme ~0 0 3, two 9*0 9"0 O'L oz Figure 4.3 Accuracy, stability and convergency test of the SSM in the solution of ablation for a medium initially at phase-change temperature imposed with a constant heat-flux condition (s) ' *wl.L 0 LO'O OVO 003L 9L 01 9 0 Oz 76 developed for the solution of the ablation problems dissertation. In fact, the good results obtained in the first examplI e are somewhat expected because of the linearization used in the solution of the integrodifferential equations for which ablation rate itself is linear. ExamplI e initially solution 2 is different subcooled for checking (see Table because 4.1). the results, it deals with This time, which a medium there is which is no exact have also been tested for convergence in Figure 4.4. Of interest in this figure is the slightly smaller ablation rate at small time, which attributed to the large side of the moving slope boundary of the temperature curve at (not shown in dissertation). the solid At large time, this slope diminishes, and finally it is stabilized evidenced by the linearity of the position curve at large time. Using yields the numerical an ablation data over the last time step in Figure rate of 0.00161 m/s, which 4.4 is in good agreement with that computed This gives using a good equation indication (4.25), error being of the ablation less than 2.5 approaching the quasi-steady state. An effort is also made to check the results of Example 2 with the analytical solutions reported in the literature. As shown Figure 4.5, both the heat integral method (HIM) of Vallerani and the moment method (MOM) of Zien [41] yield results that are in close agr results a the plot. reement with the present Llso fall right Example on the 2 has thus been source-and-sink curve method. and have thus been omitted tested successfully. with 3 and 4 deal with a in this the can be as in [40] in subcool1ed med ium imposed Examples Figure 4.4 Stability and convergency test of the SSM in the solution of ablation for a subcooled medium imposed with a constant heat-flux condition Ablated Surface Position, R1 (cm) 0 0 CA o - I0cn b 001 0o oo Fob0 0 Figure 4.5 Comparison of two methods in medium imposed condition ablated surface the literature with a constant position with for a subcooled heat-flux Ablated Surface Position, R, (cm) 0. 0- - I 3 Ii 0 CA 0 Ci 0 CA0 00 > X i . 0 "/) ic 0 I~CA UC 08 81 time-variant heat-flux Example 3, Like Example phases while 2, in the medium conditions. a quadratic the medium A linear heat-flux is subcooled for both examples. heat-flux is imposed initially; Comparison is imposed in Example there in 4. are two of the present solution is given with the moment in Figure 4.6; method for the linear convergence test based heat-flux on the change condition of time increments for this condition in given in Figure 4.7. As shown both figures, the results for the source-and-sink method are in close agreement with the moment method, and the results converged with a seen in Figures time step as large as 0.5 sec. 4.8 and 4.9 for the case of Similar results are quadratic heat-flux. The tests for ablation problems with one moving boundary have thus satisfactorily. in have been completed Figure 4.6 Comparison of ablated surface position with the moment method in the literature for a subcooled medium imposed with a linear heat -flux condition Ablated Surface Position, R1 (cm) o 0 0 0 _b 0t> xrcn 0 C CA, Ag 3 00 i' 040 0" cc %MOO - Figure 4.7 Stability and convergency test of the SSM in the solution of ablation for a subcooled medium imposed with a linear heat-flux condition (s) 1 lBWI.L OL 0"0 g9 CA 0 0 0 -9 0 9'0 6"0 ZI g'l, 9"3 V': C) 9z oz Figure 4.8 Comparison of ablated surface position with the moment method in the literature for a subcooled medium imposed with a quadratic heat-flux condition Ablated Surface Position, R, (cm) o 0A 0C O h3 CA I0i 0 " 301 , 0 CA' Figure 4.9 Stability and convergency test of the SSM in the solution of ablation for a subcooled medium imposed with a quadratic heat-flux condition |

Full Text |

PAGE 1 SOURCE-AND-SINK METHOD OF SOLUTION OF MOVING BOUNDARY PROBLEMS BY MEHDI AKBARI A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1993 PAGE 2 TO MY WIFE VAHIDEH, AND DAUGHTERS SARA AND MONA PAGE 3 ACKNOWLEDGEMENTS I am most grateful to the chairman of my supervisory committee, Dr. C. K. Hsieh, who provided the inspiration for this study as well as many hours of guidance and encouragement. His enthusiasm has been contagious, and I hope to emulate his attitude in the future. I wish to express my sincere thanks to Drs . R. A. Gater, G. G. Emch , R. Abbaschian, and H. A. Ingley whose unselfish good counsel and understanding have been essential in this study. Last but by no means least, I would like to thank my wife, Vahideh, for her love and spiritual support during this period of time of work on the dissertation. m PAGE 4 TABLE OF CONTENTS Page ACKNOWLEDGEMENT iii LIST OF TABLES vi LIST OF FIGURES Â• Â• VI 1 NOMENCLATURE ix ABSTRACT xi CHAPTERS I INTRODUCTION 1 II LITERATURE REVIEW 8 III SOLUTION OF INVERSE STEFAN PROBLEMS BY A SOURCE-AND-SINK METHOD 21 3.1 Motivation 21 3.2 General Analysis 25 Pre-melt Stage 28 Melting Stage 28 3.3 Example Problems 31 3.4 Numerical Solution 33 3.5 Critique of The Method 39 3.6 Results and Discussion 40 3.7 Extension and Concluding Remarks 56 IV SOLUTION OF ABLATION PROBLEMS WITH ONE MOVING BOUNDARY BY A SOURCE-AND-SINK METHOD 58 4.1 Solution Methodology 58 Pre-ablation Stage 59 Ablation Stage 62 Equivalent Problem 63 4.2 Illustrative Examples 65 4.3 Numerical Solution of Ablated Front 67 4.4 Results and Discussion 69 V SOLUTION OF ABLATION PROBLEMS WITH TWO MOVING BOUNDARIES BY A SOURCE-AND-SINK METHOD 90 IV PAGE 5 Page 5.1 General Analysis 91 Pre-melt Stage 91 Melting Stage 94 Ablation Stage 95 Equivalent Problem 98 5.2 Examples 101 5.3 Numerical Solution of The Ablation Problem 104 5.4 Numerical Solution of Temperature and Energy Storage 109 5.5 Numerical Examples 110 VI DISCUSSION AND RESULTS 130 APPENDICES A EXTENSION OF THE SSM 135 B SSM FORTRAN PROGRAM FOR INVERSE STEFAN PROBLEM 136 C SSM FORTRAN PROGRAM FOR INVERSE STEFAN PROBLEM 147 D SSM FORTRAN PROGRAM FOR INVERSE STEFAN PROBLEM 154 E SSM FORTRAN PROGRAM FOR INVERSE STEFAN PROBLEM 165 F SSM FORTRAN PROGRAM FOR ABLATION PROBLEM 172 G SSM FORTRAN PROGRAM FOR COMBINATION PROBLEM 186 REFERENCES 223 BIOGRAPHICAL SKETCH 229 V PAGE 6 LIST OF TABLES Table Page 3.1 Expressions to account for effects of boundary conditions 23 3.2 Properties of aluminum 41 3.3 Conditions tested in eight examples 42 3.4 Comparison between true and retrieved conditions for first Stefan-Neumann example solved for boundary temperature 43 3.5 Comparison between true and retrieved conditions for second Stefan-Neumann example solved for boundary flux then temperature 46 3.6 Comparison between true and retrieved conditions for third Stefan-Neumann example solved for boundary temperature 48 3.7 Comparison between true and retrieved conditions for fourth Stefan-Neumann example solved for boundary flux then temperature 49 3.8 Comparison between true and retrieved conditions for fifth inverse example problem solved for temperature 50 3.9 Comparison between true and retrieved conditions for sixth inverse example problem solved for heat flux 51 3.10 Comparison between true and retrieved conditions for seventh inverse example problem solved for temperature 52 3.11 Comparison between true and retrieved conditions for eight inverse example problem solved for heat flux 53 4.1 Conditions tested in four examples 70 5.1 Conditions tested in three examples Ill vi PAGE 7 LIST OF FIGURES Figure P&ge 3.1 System analyzed 27 3.2 Linearization of solid-liquid interface position curves for the numerical solution 35 4.1 System analyzed 61 4.2 Trend of ablated surface position for a medium initially at phase-change temperature imposed with a constant heat-flux condition 73 4.3 Accuracy, stability and convergency test of the SSM in the solution of ablation for a medium initially at phase-change temperature imposed with a constant heat-flux condition 75 4.4 Stability and convergency test of the SSM in the solution of ablation for a subcooled medium imposed with a constant heat-flux condition 78 4.5 Comparison of ablated surface position with two methods in the literature for a subcooled medium imposed with a constant heat-flux condition 80 4.6 Comparison of ablated surface position with the moment method in the literature for a subcooled medium imposed with a linear heat -flux condition 83 4.7 Stability and convergency test of the SSM in the solution of ablation for a subcooled medium imposed with a linear heat-flux condition 85 4.8 Comparison of ablated surface position with the moment method in the literature for a subcooled medium imposed with a quadratic heat-flux condition 87 4.9 Stability and convergency test of the SSM in the solution of ablation for a subcooled medium imposed with a quadratic heat-flux condition 89 VI 1 PAGE 8 Figure Page 5.1 System analyzed 93 5.2 Linearization of solid-liquid interface and ablated surface position curves for the numerical solution 106 5.3 Trends of solid-liquid interface and ablated surface positions for a combination problem of subcooled medium imposed with a constant heat-flux condition 114 5.4 Trends of solid-liquid interface and ablated surface positions for a combination problem of subcooled medium imposed with a linear heat-flux condition 116 5.5 Trends of solid-liquid interface and ablated surface positions for a combination problem of subcooled medium imposed with a quadratic heat-flux condition 118 5.6 Temperature distributions in the medium at different times during ablation for a combination problem of subcooled medium imposed with a constant heat-flux condition 120 5.7 Stability and convergency test of the SSM in the solution of ablation for a combination problem of subcooled medium imposed with a constant heat-flux condition 122 5.8 Stability and convergency test of the SSM in the solution of ablation for a combination problem of subcooled medium imposed with a linear heat-flux condition 124 5.9 Stability and convergency test of the SSM in the solution of ablation for a combination problem of subcooled medium imposed with a quadratic heat-flux condition 126 5.10 Overall accuracy test of the SSM in the solution of the combination problem of subcooled medium imposed with a constant heat-flux condition 129 VI 1 1 PAGE 9 NOMENCLATURE c E F G r Ste s T t v X a 6 Specific heat Equation Temperature function used for boundary condition Flux function used for boundary condition, GreenÂ’s function, elements in coefficient matrix Heat flux Thermal conductivity Latent heat of fusion Latent heat of vaporization Normal unit vector Distance from sink or source to sense point Interface positions Stefan number Dummy variable for integration, interface speed Temperature, temperature vector Time velocity of the interface Position, source point Thermal diffusivity Dirac delta function, Kronecker delta IX PAGE 10 p Density r Dummy variable for time Subscripts and Superscripts m,v Phase change N 19 N 2 Upper time limit index n Time index t Transition 1,2 Melt and vapor fronts 0 Initial time x PAGE 11 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy SOURCE-AND-SINK METHOD OF SOLUTION OF MOVING BOUNDARY PROBLEMS By Mehdi Akbari May 1993 Chairman: Dr. Chung K. Hsieh Major Department: Mechanical Engineering Source-and-s ink methods have been developed for the solution of inverse diffusion problems, ablation problems, and combination of ablation and Stefan problems. GreenÂ’s functions have been used and integrodif f erential equations are derived for the interface positions and the temperature in the phase-change medium. These equations are solved by using local linearization of the interface position and the boundary heat flux if it were treated as unknown. The results have shown to be accurate, convergent, and stable. The methods developed for the solution of the inverse diffusion problems have been used to find the boundary conditions for the inverse Stefan problems. They are solved by two approaches: a series solution approach and a time incremental approach. Both have shown to be useful to find the boundary conditions without reliance on the flux information to be supplied at both sides of the phase change interface. The methods are efficient in that they require fewer equations to be solved for the xi PAGE 12 unknowns. The algorithms can also be easily developed for the solution of the problems. In the use of the source-and-sink method to analyze the ablation and combination problems, heat transfer is solved in a fixed domain so that the original heat-flux condition that is imposed on the moving boundary is taken to be the condition imposed on the solid-vapor or liquid-vapor interface. This flux condition is then used together with the temperature at the interface to solve for the interface position as well as the condition imposed on the fixed boundary. Finally, these positions and conditions are used in the temperature equation to complete the solution. The method has been used successfully in solving examples encompassing one-, two-, and three-phase ablation with the boundary of the medium imposed with constant, linear, and quadratic flux conditions. Numerical solutions for these examples as well as the trends of the moving boundary positions are discussed in detail. xi 1 PAGE 13 CHAPTER I INTRODUCTION Phase changes occurs when a body is exposed to a large heat flux so that it melts or vaporizes at the surface. With continuous heating, the melting or vaporization front moves inward with time. Since the position of the phase-change interface is unknown a priori, this type of problem is commonly known as the moving boundary problem in the heat transfer literature. Phase change problems find applications in thermal energy storage, material treatment and processing, thermal protection of re-entry vehicles in space technology, laser drilling and cutting in manufacturing, freezing and thawing in soils and foodstuff, and photocoagulation in opthalmological procedures, among others. It is difficult to solve a phase-change problem because of the presence of the nonlinear condition at the interface. Usually for such a problem, the number of phases encountered in the medium depends on the ambient pressure and the temperature range that appears in the medium. Thus for a medium heated or cooled below the triple point of the substance, only sublimation and fusion may take place. On the other hand, for phase change above the triple point, solid, liquid, and gas states may appear simultaneously, and this occurs for a medium of low thermal conductivity and exposed to a large heat flux. In practice, the changed phase may be removed as 1 PAGE 14 2 soon as it is formed. This occurs in the re-entry where the melted substance is removed by a fluid dynamic force. This material removal may also be effected naturally such as in sublimation, where the sublimed vapor is dissipated by diffusion to the surroundings. In either event, the imposed heat flux follows the boundary motion, and such problems are commonly referred to as ablation problems in the literature. A different state of affairs is encountered when the changed phase remains stationary and continues to occupy the space taken by its previous state. Thus for example in melting a solid, if the melt is not removed by an external force, it will adhere to the solid from which it changes phase. Then in the absence of convection in the melt, the heat transfer in the medium can be analyzed by solving a Stefan problem in recognition of Stefan who first presented a formal solution of the problem in the open literature (see next chapter for a survey of the literature). The Stefan problem is basically different from the ablation problem in that, in the former, the problem is solved in a fixed domain, and the boundary condition is always imposed on the fixed boundary whereas, in the latter, the imposed condition follows the boundary motion. Thus in the ablation problem, the domain for analysis diminishes in size with time, and the imposed condition is acting on the receding boundary due to the shrinkage of the medium. There have been numerous efforts in the literature for the solution of the Stefan problems. Much less research, however, has been devoted to the solution of the ablation problems. Also for those ablation problems that have been solved and documented in the literature, PAGE 15 3 most are for the analysis of generic ablation in which only one moving boundary appears in the medium. There has been a lack of study for problems involving two moving boundaries such as those found at the solid-liquid and liquid-vapor interfaces. They are encountered when the solid, liquid, and gas states appear simultaneously in the medium. Three phase-change problems will be solved in this dissertation. In the first problem, an inverse Stefan problem will be solved. This problem differs from a conventional Stefan problem in the sense that, for the conventional problem, the condition imposed on the boundary is fully specified. This condition is then used together with others to solve for interface motion as well as the temperature distribution in the medium. For the inverse problem, however, the interface motion is given, and this motion is used as a part of the input to solve for the boundary condition that is necessary to provide for this interface motion. This problem is found in material treatment and processing, where the property of the material may be a function of the phase-change rate at the interface. Solution of the problem thus enables a better control of the property of the material . It should be noted that the inverse Stefan problem described above has recently been solved by Zabaras et al . [1] in the literature. However, the method developed for solution in this dissertation is more efficient because it requires no information for the heat flux at the interface. As will be shown, two methods are developed for the solution of the problem: one expands the condition in a polynomial or an infinite series, and the other uses PAGE 16 4 a time incremental method. The former utilizes the interface position information at discrete times to solve for the coefficients in the polynomial or series expansion, while the latter uses the interface position at consecutive times to solve for the conditions incrementally. Both can be used to solve for the conditions accurately as supported by eight examples encompassing the search for constant and time-variant temperature and flux conditions at the boundary. The methods are also well adapted to a numeral solution. The second phase-change problem solved in this work is an ablation problem that contains only one moving boundary. This problem is solved by a new method, which treats the problem in a fixed domain. In this method, the ablated phase is considered to remain in the space that is ablated in the process. As such, the moving boundary can be taken to be an interior phase-change interface, and the heat flux that is derived at this interface can be used to match the condition imposed at the moving boundary. Then with the additional prescription of the temperature for phase change at this interface, the flux condition there can be used to determine the position of the interface as well as the hypothetical condition that is imposed at the fixed boundary. Finally these positions and conditions can be used in the temperature equation to complete the solution . Solution of the ablation problems has rarely been attempted in a fixed domain. The method presented in this work thus represents a new endeavor, which will be shown to be simpler than most methods described in the literature. To validate the method used in the solution of the ablation problems, four examples are provided and PAGE 17 5 they include the analysis of oneand two-phase ablation imposed with constant and time-variant flux conditions at the moving boundary. Of those four, the ones imposed with constant heat flux can be checked with an exact solution. A one-phase ablation problem is thus selected for checking the solution over the entire period of the ablation, while a two-phase ablation problem is selected for testing the solution in a quasi-steady state. Both yield accurate results as shown later in this dissertation. Convergence and stability tests have also been made rendering further assurance for the success of the method. The good results obtained in the second problem described above provide impetus for the solution of the third problem that addresses ablation with two moving boundaries. This is a combination of the ablation and Stefan problems, an area that has rarely been studied in the literature yet is important in practice when solid, liquid, and gas states appear simultaneously in the medium. Again the problem is solved in a fixed domain, and the moving boundary is treated as an interior phase-change interface. The flux condition at this boundary is again used with the temperatures at the interfaces to determine the heat flux at the fixed boundary as well as the interface positions to complete the solution. Unlike its predecessor for which two simultaneous integrodif f erent ial equations are solved simultaneously for each time step, for the combination problem, three equations will be solved simultaneously. Again the results are good as evidenced by three examples including imposition of constant, linear, and quadratic flux conditions. PAGE 18 6 A source-and-sink method will be used for the solution of all the problems in this work. This method has recently been subjected to a through development [2-5], where the method has shown to be particularly suited for the solution of the phase-change problems addressed in this work. In this method, the problem is solved by using one set of governing equation, initial, and boundary conditions, whereas in the conventional methods developed in the literature, the problem is solved by using different sets of equations, with each set focused on the solution of one phase region. The conventional methods thus cannot compete with the source-and-sink method for simplicity. Moreover, in the source-andsink method, only one temperature equation will be derived; whether it is in the solid or liquid region depends on the position that is assigned in the equation. This is in sharp contrast to the conventional methods, in which several temperature equations must be derived, again one for each phase region. The time and effort saved in the source-and-sink method can thus be readily appreciated. With the inclusion of the three separate problems in this dissertation, the presentation of materials follows a nontradi tional approach. Thus one chapter will be devoted for the presentation of each problem; and the chapters will be self -complete with the coverage of the statement of the problem, general solution methodology, analysis, examples, and results and discussion. It is felt that, only through such an organization, the material can be presented in a coherent manner without severe fragmentation. In what follows in the next chapter, a literature review will be given first. It is then followed by the presentation of the three PAGE 19 7 problems. Finally, conclusions and recommendations will be given to conclude this dissertation. As usual, all computer programs developed in this work are compiled in the appendix. It is hoped that the insight gained from this research will be useful in the analysis of the phase change in the future. PAGE 20 CHAPTER II LITERATURE REVIEW Heat conduction with phase change constitutes a large class of moving boundary problems. These problems are commonly referred to as Stefan problems and were first solved in 1831 by Lame and Clapeyron [6], who determined the thickness of the solid crust by freezing of water under a constant temperature condition. They were able to find the thickness to be proportional to the square root of time but did not find the constant of proportionality. Some thirty years later, Neumannn presented an exact solution to the problem in his unpublished lecture notes, and for the first time in history, the Stefan problem was solved in its entirety [7-8]. NeumannÂ’s solution was for the phase change in a semi inf ini te space imposed with a constant temperature condition. His work is important because, irrespective of the numerous efforts made by many others in search for additional exact solutions, Stefan problems that can be solved exactly today remain those that were originally analyzed by Neumann some one hundred and thirty years ago. While the exact solution to the Stefan problem was first obtained by Neumann, the problem has been named after Stefan in recognition of his published work appearing in the open literature in 1889 [6-9]. Practically the same problem was solved by Stefan, who expressed the temperature of the phase-change medium in terms of 8 PAGE 21 9 a similarity variable, x/^ ? and the position of the moving boundary to be proportional to the square root of time. Thus, in essence, the moving boundary was solved by means of a similarity transformation. Since the transformation can not be used for general problems, it has now been firmly established that, for a Stefan problem that can be solved exactly, it must be in an unbounded domain, of medium of constant properties, and imposed with a constant temperature condition. With the severe limitation imposed by the similarity transformation, it is not surprising to see a wide variety of approximate solutions developed in the literature. Books and monographs have been prepared for presentation of these methods. In what follows, some popular solution methods will be briefly reviewed; interested readers are referred to 10-14 for details. One of the early methods developed for the solution of the Stefan problems imposed with heat flux and convection conditions is the power-series expansion method. In this method, the interface position and the temperature in the medium are expanded in terms of power series [10] or complimentary error functions [11-14]. These series are then substituted into the governing equation and boundary conditions to determine the coefficients in the series expansion. Power series method works well for the short-time solutions. At that time, the series converge rapidly, and only a few terms of the series are sufficient to yield accurate results. The method also works for the solution in the neighborhood of singularities that may occur as a phase degenerates. For large time, more terms are needed in the series, and the method rapidly loses its appeal. PAGE 22 10 Stefan problems can also be solved by reducing them to integrodif f erential equations . Lightf oot [15] took the solidification as a moving heat source front. With the use of GreenÂ’s function, he was able to derive two integrodif f erential equations, one for the temperature and the other for the interface position. Then by assuming the position a function of the square root of time enabled him to retrieve the Stefan-Nuemann solution. LightfootÂ’s method has been extended to the solution of two dimensional Stefan problems involving phase change in an infinite wedge [16]. To avoid a tedious iteration in two dimensions, the interface position was taken to be a hyperbolic function and a correction term was introduced into the solution to account for the property variation in different phase regions. Indeed, LightfootÂ’s method has been limited to phase change of equal properties. This limitation has been lifted by Kolodner [17] by using double source and sink at the phase-change interface. Integrodif f erential equations of Volterra or Fredholm type can also be derived in the solution of the Stefan problems by a Fourier transform [18-19]. This occurs when the transformed equations are inverted in the solution. In general, the integral equation method is useful in providing an exact solution in integral form. However, the integral equations must still be solved numerically for a solution . One of the unique features of the Stefan problems is the occurrence of multiple phase regions whose domain changes continuously with time. The boundary position of the domain is also sought as a part of the solution. Boley [20] introduced an PAGE 23 11 embedding method, in which the time-variant domains are embedded in a large fixed domain for the solution of the problem. This introduces an unknown heat flux at the fixed boundary, and this flux must be solved together with the moving boundary position to complete the solution. BoleyÂ’s embedding method has been extended to the solution of multidimensional problems without internal heat generation [21]. It should be noted that all the moving boundary problems solved in this dissertation also work in a fixed domain; however, a source-and-sink method will be employed which is more convenient to use in the solution of the combination of Stefan and ablation problems. Stefan problems can also be solved by an asymptotic expansion [22-29]. In this method, a quasi-steady solution is derived by dropping the unsteady temperature terms in the governing equations to find a long-time solution. A quasi -stationary solution is also derived by dropping the interface velocity terms in one of the governing equations and the interface flux condition to establish a short-time solution. Asymptotic expansions of the temperature and the moving boundary position are then constructed by using these solutions for limits [26]. Asymptotic expansions work best in the solution of Stefan problems with nonlinear conditions and convective motions of the melt [22]. Basically a perturbation technique, it offers insight into the physics of the problem. However, the method is very tedious mathematically; even the determination of the first-order terms in the temperature expansion has proven to be an intractable task [26]. PAGE 24 12 Higher-order terms become increasingly difficult to obtain, a distinct drawback of the method. Coordinate transformation offers another approach to the solution of the Stefan problems [30]. Somewhat similar to the embedding method, the transformation works in a fixed domain without the need of introduction of any new unknowns such as the boundary heat flux in the embedding method. In fact, the time-variant domain is mapped into an invariant domain in the transformation method; the moving boundary is thus immobilized in the solution. It is noted that use of the transformation does not by itself solve the Stefan problem. It only provides the convenience of working in a fixed domain. The coordinate transformation has thus been used together with other methods for the improvement of accuracy and also in the numerical methods for facilitation of solution [31-37]. Like the Stefan problems briefly reviewed above, ablation problems also fall into the general category of the moving boundary problems. However, the ablation problem is more complex than the Stefan problem because in the latter the problem can be solved in a fixed domain, whereas in the former the domain is continuously changing with time. There have been numerous efforts developed for the solution of the Stefan problems; only very limited work, however, has been devoted to the solution of the ablation problems, which will now be reviewed as follows. Solution of the ablation problems dates back to Landau [38] who first used the coordinate transformation to change the variable domain encountered in the ablation problem to a fixed domain. An PAGE 25 13 exact, quasi-steady state solution for the case of a semi inf ini te medium imposed with the constant heat flux condition has been given in his paper. Also studied is the ablation for the semi inf ini te medium after dropping the steady-state assumption. Using Laplace transformation, Landau was able to derive the pre-melt solution for a time-variant heat flux imposed on the fixed boundary. However, the ablation solution was only derived for the constant heat flux condition. Specifically, two limiting cases were examined that include (i) Stefan number approaching zero (negligible heat capacity) and (ii) Stefan number approaching infinity (negligible latent heat). For the case of the Stefan number of the order of unity, Landau employed a finite difference method to solve the ablation problem with results presented graphically. LandauÂ’s method has also been used by Rogerson and Chayt [39] to find the exact melt-through time for ablation of a slab imposed with a constant heat flux condition on one side and an insulated condition on the other side. Rogerson integrated the heat conduction equation and showed the results to be independent of the thermal properties of the material engaged in phase change. The integral approach has been used in the solution of the phase-change problems [9,40,41]. Essentially a method of weighted residuals, the method was first introduced by von Karman and Pohlhausen in the approximate solution of boundary layer equations. The heat integral approach is simple to use; it also provides reasonably accurate results. However, it is somewhat handicapped in a detailed analysis of the temperature field. Specifically, the accuracy of the temperature is limited by the form of the profile PAGE 26 14 initially chosen for analysis. Also the approximation can not be systematically improved for accuracy. The moment method proposed by Zien [41] carries promise of improvement over the classical heat integral method. The moment method has been used to solve one-dimensional ablation imposed with time-dependent heat-flux conditions. Both pre-melt and ablation solutions were derived. Again, the heat balance integral was used, and the integral of the original heat equation was carried out after multiplying the integrand by powers of temperature. A temperature profile was chosen and substituted into the integrated version of the heat conduction equation. The heat balance integral based on the approximate temperature profile was then used as the expression for the boundary heat flux. Although the method appears to work for the general case of the time-variant heat flux condition, it is expected that the choice of an approximate temperature profile that works for the nonmonotonic heat flux may not be applicable to the case of more general time-variant conditions. The ablation problems can also be solved by a variational approach. Biot and Agrawal applied the variational analysis and Lagrangian thermodynamics to the solution of ablation problems with variable thermal properties [42]. They considered one-dimensional heat transfer in a semi-infinite cylinder imposed with a constant heat flux condition. Both pre-ablation and ablation stages were solved. In the procedure, the governing equations were transformed and the Lagrangian heat-flow equation was derived which provided a relationship between the surface velocity and the heat penetration depth. Another relation between these two quantities was also PAGE 27 15 found by using the energy equation. These equations were then solved simultaneously for the heat penetration depth and the velocity of the ablated surface. The variational method has been applied for approximate analysis of ablation of a semi inf ini te solid subject to convective and radiative heating [43,44]. Solutions were obtained in closed form for both the pre-melt and melt-removal heating regimes. In these studies, a cubic temperature profile was taken. The surface temperature, the thermal penetration depth, and the depth of the melt removal were treated as unknown and were determined as functions of time. The ablation problem has been solved numerically. Recently Blackwell [45] has employed the exponential differencing method to solve an ablation problem in one-dimension. He proposed the use of a moving coordinate system which was attached to the ablated surface. Then, by invoking the use of a finite control-volume approach, the element matrices could be defined for conduction, convection of the moving grid, and energy storage. In this method, all elements and control volumes were moving at a uniform velocity with the exception of the last ablating element and control volume. The node point at the moving interface was fixed in space. As such, the last ablating element had one moving boundary and one fixed boundary. The method has been applied to the solution of a steadystate ablation problem for which an exact solution was available. Comparisons were also made with central differencing for the conduction terms and upwind differencing for the convective terms; the exponential schemes appear to be better numerically. PAGE 28 16 Ablation problems have also been solved by a finite element method with deforming spatial grids [46]. In this method, the classical finite-element equations are transformed to account for the continuous deformation of the grid for a precise localization of the ablated surface. This is done at the expense of chores of construction of an additional convective matrix. Recently, ablation has been applied to model intense heating such as laser beam. Masters [47] used the finite difference method to solve the ablation of a one dimensional slab imposed with an intense uniform heat flux and analyzed the effect of melt on the temperature distribution during the heat pulse. The steady-state solutions were found for the velocity of the surface recession and the temperature history during ablation. Abakian and Modest [48] studied the ablation due to a continuous-wave laser beam irradiating on a moving semi inf ini te and semitransparent solid. Using an integral method, they were able to derive a set of nonlinear partial differential equations which were solved numerically for the groove depth and shape due to ablation. They also considered the ablation of a moving slab caused by irradiation from continuous-wave and pulsed laser beams and derived a solution for the temperature distribution [49]. The laser has also been used in material processing as investigated by Dabby and Paek [50]. In their work, vaporization occurred at the surface; however, below the surface, the material was heated by absorption of the laser radiation, which might reach a temperature higher than that for vaporization. Explosion may thus take place, which provides a means for material removal in the drilling process. PAGE 29 17 The review given above is strictly for the solution of regular heat conduction problems. In these problems, the system geometry, governing equation, initial and boundary conditions are fully specified, and the problems are solved mainly for the determination of the temperature field. In the inverse problems, however, the roles of known and unknown quantities are exchanged. There have been abundant studies documented in the literature for the solution of the regular problems; not much work, however, has been devoted to the solution of inverse problems, and for those that have been solved and documented in the literature, nearly all of them are for heat conduction without phase change [51-60]. There is a lack of study for inverse problems with phase change. Most inverse heat conduction problems deal with a situation where an extra temperature is available at one point in the domain, and this temperature is used together with others to find the condition imposed on the boundary. Stolz [51] first solved such a problem numerically. In his work, the inverse problem was formulated as though it were a direct problem. Since a linear heat conduction problem was solved, he was able to use the superposition principle. An integral equation was derived for the unknown surface condition and was solved by numerical inversion. The method was found to be inefficient if the time steps used were too small. Yet, a small time step must still be used for an accurate solution. Beck [52] was able to improve this method by using a procedure that involved minimization of the sum of the squared difference between the actual and the calculated temperatures at the location where the temperature data were given. Burggraf [53] developed an exact PAGE 30 18 series solution to a one-dimensional inverse problem with the lump capacitance approximation serving as the leading term. The temperature and heat flux histories were provided at an interior point. Approximate results were found if discrete or experimental data were used for input in solution. Beck [54] moved one step further by solving an inverse problem in which the material properties were treated as functions of temperature; thus the problem solved becomes nonlinear. A two-dimensional inverse finite element solver has also been reported in the literature [55]. It can be used for the solution of the heat flux imposed on the surface of nuclear fuel rod with the use of the interior temperature measurements for input. Only a handful of studies are found for the solution of inverse heat transfer with phase change. Macqueene et al . [61] proposed an inverse finite element method to determine the efficiency of an arc welding process. In their method, the latent heat of fusion is taken into consideration by the variation of the elemental specific heat when the average temperature of the element reaches the melting point. Conduction in both the solid and liquid regions was accounted for. Katz and Rubinsky [62] proposed a front-tracking finite-element method for the solution of one-dimensional inverse Stefan problems. His method was applied for the determination of the position of the solid-liquid interface and the transient temperature distribution in the solid region during stationary arc welding . An inverse method was also used by Landram [63] for the analysis of the interface position and the energy transport PAGE 31 19 mechanisms during welding. The vaporization energy loss was found to be important during the motion of the solid-liquid interface. The interface shape, being close to hemispherical, gives clear indication of a nearly one-dimensional radial symmetry for the heat transfer . An analytic solution to inverse Stefan problems in Cartesian and spherical geometries was provided by Rubinsky and Shitzer [64]. In their analysis, the inverse Stefan problem was characterized by two boundary conditions at the moving front. This is in sharp contrast to the technique developed in this dissertation (see Chapter 3) which needs only one condition at the interface to solve for the boundary condition. In their work, the medium was initially at the phase change temperature. An integral equation was then derived by integrating the governing equation, which was, in turn, solved by the method of analytic iteration in which the first guessed solution was taken to be the long-time solution to the problem. Series solution was then developed by induction following a number of iterations. A boundary element analysis with constant elements has been developed by Zabaras et al . for the solution of a one-dimensional inverse solidification problem [1]. They used the sensitivity analysis developed by Beck [52-54] and Burggraf [53] for inverse heat transfer solution. Using an integral formulation, they were able to solve two separate inverse Stefan problems, one in the solid region and the other in the liquid region. Thus similar to the ones given by Rubinsky and Shitzer, two conditions must be provided at the freezing front, also an inefficient method. PAGE 32 20 It should be mentioned that a source-and-s ink method has recently been developed in the heat transfer literature, which is particularly suited for the solution of regular and inverse Stefan problems. In a series of papers, Hsieh and his associates were able to apply this method to the solution of regular and inverse oneand two-phase melting and solidification problems for medium with and without subcooling and superheating and imposed with constant and monotonic temperature and heat flux conditions [2,3]. The method has been applied to the solution of phase change imposed with cyclic conditions [4,5,65]. Two inverse solution techniques have also been developed with the problem formulated with a source-and-s ink method as shown in reference 3. The method has shown to be closely related to the boundary element method as reported by Hsieh et al . [66,67]. The present study is a further extension of these works. PAGE 33 CHAPTER III SOLUTION OF INVERSE STEFAN PROBLEMS BY A SOURCE-AND-SINK METHOD This chapter presents the development of a source-and-sink method to solve inverse Stefan problems. In these problems, the conditions are specified at the moving rather than the fixed boundary. Typically, the temporal location of the interface between the phases is given and this location is used together with others to determine the boundary temperature or heat flux that is required to provide for this interface motion. In what follows in this chapter, the motivation for this study will be given first. It is followed by a general analysis that is designed for the solution of the inverse Stefan problems in three dimensions. This analysis is then used in the solution of onedimensional inverse problem examples. Finally, results for these examples are provided and discussed in detail and possible extensions of the method are included to conclude this chapter. 3 . 1 Motivation Heat diffusion in a medium with constant properties is governed by the partial differential equation V 2 T(f,t) + _ 1 5T(r,t) k Â“ 01 dt Â’ r G R t > 0 ( 3 . 1 ) 21 PAGE 34 22 where all notations have their usual meaning. For this medium, the conditions imposed on the boundary are usually one of the following types T(r,-,t) = F,(r,-,t), r,G BdT(f,,t) G,.(f ,,t) ~s^r Â— r ' Â€B ' GB, (3-2) (3.3) (3.4) which represent the familiar Dirichlet, Neumann, and Robin conditions, respectively. In (3.3) and (3.4), n t denotes an outward drawn normal. Then, with the additional initial condition given as T(r,0) = T-(r) (3.5) the temperature solution can be expressed by means of GreenÂ’s function as [68] t T(f,t) = G(f,t | f',0)T.(f / )dV / + f [ G(f,t | f',r )dV , dr + E{ } (3.6) R' Or' Here, the braced term is used to account for the three boundary conditions given earlier. Their expressions are listed in Table 3.1. A distinct feature is found in the GreenÂ’s function method above--the effects of the initial condition, heat generation (or destruction), and boundary conditions are embodied respectively in PAGE 35 23 Table 3.1 Expressions to account for effects of boundary conditions BOUNDARY CONDITION { } EXPRESSION T(T' / t) -F i (T' , t) -aff 3G PAGE 36 24 the first, second, and third terms on the right of equation (3.6). Then, in the case of regular problems, once these conditions are fully specified, the temperature can be easily found. In this effort, the GreenÂ’s function can be obtained by using the concept of point charges [68] or by solving auxiliary problems [69]. The format of (3.6) turns out to be particularly suited for the solution of the inverse problems. For the inverse problems, the left hand side of this equation can be used to represent the extra temperature that is provided either at the boundary or at interior points. Then this temperature can be used to determine the missing information, such as the initial condition, heat generation, or boundary conditions. In these efforts, the missing quantities can be expanded either in a polynomial or an infinite series, which is substituted into their respective integrals in (3.6). The resulting equation is then solved numerically for the coefficients in the series to complete the solution. This method clearly works for the initial condition and the heat generation. However, for the boundary condition, since it is unknown a priori, one must first assume a particular 4 type Â’ of condition that is imposed on the boundary. For convenience, one could use equation (1.2) or (1.3) for this condition. GreenÂ’s function is then found with this assumed condition. Notice that the boundary condition can be exchanged if necessary as shown in by Hsieh and Shang [70]. That is to say, a Dirichlet condition can be accomplished by means of a Neumann condition and vice versa. The search for the condition is thus not restricted by the type of the conditions assumed. Better yet, an incremental solution approach can also be developed to track PAGE 37 25 the boundary conditions accurately as will be shown later. The concept above will now be applied for the solution of inverse Stefan problems in which the interface positions are given, and these positions are used to find the boundary conditions that must be imposed to cause the interface motions. 3.2 General Analysis For the sake of illustration in what follows, the Stefan problems consist of two stages: a pre-melt stage, when heat is added to the surface of a subcooled medium to raise its temperature to the phase-change temperature; and a melting stage, when the medium changes phase and the melting starts at the surface and the interface moves inward with time. It is assumed that the properties for different phases are constant and of equal value. The medium has a distinct melting temperature; that is, no mushy zone in the medium. Convection is negligible. For generality, the analysis will be developed for the solution of both melting and solidification problems. The analysis can also be extended for the solution of Stefan problems in multiple phases and in a medium with unequal phase properties as will be discussed later. For the moment, the simplified problems will be solved by considering the medium shown in Figure 3.1. The formulation of these problems follows below. PAGE 38 Figure 3.1 System analyzed PAGE 39 27 PAGE 40 28 Pre-melt Stage Governing equation-V *r 0 (r,t) = i 1 ST 0 (r,t) dt r G (LUS) t 0 > t > 0 Initial condition-T o (r,0) = T-(r) Melting Stage Liquid Region: Governing equation-V *r L (r,t) = i 1 0T L (r,t) dt r G L t > t Q Solid Region: Governing equation-V 2 T 5 (i,t) 1 9T s (r,t) a dt fÂ€S t > t 0 (3.7) (3.8) (3.9) (3.10) PAGE 41 29 Initial condition-Ts( r >t 0 ) Â— T 0 (r,t 0 ) (3.11) Interface Conditions: T L (T f ,t) = T m = T 5 (r /)t ) (3.12) dT 5 (f/,t) dT L (t f ,t) _ pL dn On ~ k 6} v Â„(t)=V.n (3-14) Here ry denotes the interface position and v n the history of the interface motion. For the inverse Stefan problems of interest in this study, this history is used for the determination of the missing boundary conditions. The problems as posed can be solved by use of the GreenÂ’s function method described in the preceding section. However, a direct use of this method would require (3.6) to be applied to two separate regions, liquid and solid, and the solution so obtained may not be as efficient as one desires. A source-and-s ink method is thus used [2,4,15]. In this method, the melting interface is taken to be a moving heat-sink front and a freezing interface is taken to be a moving heat-source front. Then, in sharp contrast to conventional methods in which different equations are used to represent the temperatures in different regions, only one equation will be derived. Whether it is in the solid or liquid region is PAGE 42 30 determined by the position that is assigned in the temperature equation. The solution of the inverse problems can then be simplified with this method. Following this approach, the melting stage is solved by considering an equivalent problem as follows: Governing equation for the equivalent problem: V 2 T(r,t)+^ v n (t)5(r-r / ) = Â£ P L _ _ , i <9T(r,t) dt f G (LUS) (3.15) t > t 0 Initial condition for the equivalent problem: T(f,t 0 ) = T 0 (r,t 0 ) Interface conditions for the equivalent problem: (3.16) T(f / ,t) = T m , vÂ„(t) = V.n (3.17a, b) where S( rÂ— Fy) denotes a Dirac delta function. The signs preceding this function are used for freezing (+) and melting (-). It can be shown readily that (3.15) reduces to (3.9) and (3.10). Furthermore, by integrating (3.15) across the interface from Fy Â— e to Fy + e and forcing e to be zero in a limiting process, equation (3.15) reduces to (3.13). This can be proved by using the pill box at the interface in Figure 3.1. Other equivalences for this stage are apparent. PAGE 43 31 The equivalent problem can be solved by referring to (3.6) in which the heat generation term is changed to the interface motion term as T(r,t)= G(r,t | r , > < 0 )T t (f / )dV / LuS G(r,t |r , ,r)v n (r)i(r'-r / )dV , dr+ E{ } t Q L U S (3.18) where the plus sign is used for freezing and minus sign for melting. Finally, the boundary conditions can be found by setting r in this equation to the interface position, Fy, and T(ry,t) to the melting temperature, T m , as T = m G(r / ,t|f',< 0 )T t .(f / )dV , Â±^[ [ G(f / ,t|r / ,r)v n (r)5(r'-r / )dV / dr+ E{ } LuS t Q Lu5 (3.19) The missing boundary conditions can then be found by solving them implicitly. In this effort, the time when melting starts (t 0 ) can be determined by solving the pre-melt problem, whose solution can again be taken to be the special case of (3.6) in which the heat generation term is zero. Solution in this stage is thus elementary. 3.3 Example Problems The analysis above is now used to solve example problems which are in semi inf ini te domain in which the interface motion is given. For the sake of generality, the analysis will be developed for determining either the Dirichlet condition or the Neumann condition that is imposed on the boundary at x=0 , and the interface motion may PAGE 44 32 be the result of either melting or solidification in a oneand twophase medium. Also for illustration purposes, the heat flow is one dimensional, the initial temperature being uniform. Extensions to more dimensions are provided in Appendix A. For the problem given, GreenÂ’s function can be found to be G(x,t | x',r)= 1 2^7ro(t Â— r) [ ex P( Â“ (x-x') 2 N ^ , (x + x') 2 3S(t^) )Â±exp< -4^(t3T) )] (3.20) where the plus and minus signs are to be used when the flux and temperature condition is assumed to appear at the boundary, respectively. The temperature can then be obtained by using (3.18), which is recast in a general format as G(x,t | R(r-Ft 0 ),r)dr (3.21) where all temperatures, including T m , are measured in excess of the initial temperature, and T 0 ( X ,t) Â— > 7T t 0 E ( s ) r x 2 ex p[ Â— -jTu ~ ]ds -si 1 / 2 4a(t Â— s) (t-s) ( 3 . 22 ) in which jlZM 2a t Â— s E(s) = ^ (3.23a,b) r G ( s ) Here F(s) and G(s) denote the assumed temperature and heat flux condition, respectively. Also PAGE 45 33 1 t > t 0 for 0 t < t 0 (3.24) (3.25) where Ste is known as the Stefan number. With the use of the circumflexed Heaviside function given by (3.24), equation (3.21) holds for all time and for both oneand two-phase problems. Equation (3.21) can be used to determine the unknown boundary condition by invoking use of the condition at the interface as t-t Â±[iT Â° ( y t),t) ] 0 H(t 1 0 ) f dR(r+t 0 ) m Ste dr G(R(t),t t 0 | R(r-ft 0 ),r)dr (3.26) 0 where the plus and minus signs on the left hand side are to be used for freezing and melting, respectively. 3 . 4 Numerical Solution A local linearization can be used to solve (3.26) numerically. In this effort, the entire time range is divided into small increments in which the interface position is taken to be linear; see Figure 3.2. Then the dR/dt can be taken out of the integral for each increment, and the convolution integral written as a summation as PAGE 46 Figure 3.2 Linearization of solid-liquid interface position curves for the numerical solution PAGE 47 35 0 T PAGE 48 36 i T 0 (R(t Ar ),t 7V -) l H (t N t 0 ) N, dR(t n ) T m ~ Â±1 Ste *Â»Â» t 0 dt G(R(t N ),t N t 0 | R(r+t 0 ),r)dr] ^n-l Â“ ^0 (3.27) Here the signs are to be selected following the statement below (3.26). For the inverse problems, the right hand side can be evaluated with the input of the interface motion data. As for the left hand side, if the boundary conditions are represented by a power series, then the number of terms on this side must be equal to the number of the terms that are taken in the series. In practice, equation (3.27) can be written by means of matrix elements for a melting problem for the series method as N N 53 ^mn a n C ra 53 ^mn e n n=l n=l (3.28) where m=N ; t^>t 0 ; N=1,2,...,N; and * 0 for m < n N for m>n; F(s) = ^ a n 5 n n=l exp[ R(U 2 Mt ra -s) for m > n ; G(s) = 5Z a n s " 1 n=l c m ~ N 7 T rp a 1 m (3.29a,b,c) (3.30) PAGE 49 37 r 0 for m < n = < mn r 7 T t -t n-1 G(R(t m ),t m -t 0 | R(r + t 0 ),r)dr for m>n 0 (3.31a,b) dR (.S) (3.32) n dt v ' Notice that the integrals in these equations can be performed in closed form if the missing conditions are expressed in the power series as shown in (3.29) or in terms of a Fourier series. The coefficient matrices B and D are of a lower triangular structure, a characteristic permitting the solution of (3.28) to be carried out simply by using forward substitution, a simple numerical procedure. The boundary conditions can also be determined by use of an incremental approach. In this method, (3.27) is recast as t N t N-1 2L T / 41 _ Ip) yv r A m\ Ste n 4=l dR(f n ) Â“ t 0 dt G(R(t yv ),t N -t 0 | R(r+t 0 ),r)dr]jt N-1 n { jv Â— j Â±Â£ n=l E(s) _.^/2 ^n-l Â“ ^0 R{t N y \ T72Â“ p1 -iÂ»( 1 Â„-s) 1< ' s ) n Â— 1 (t N Â— s) 4a(t N -s) (3.33) where t 0 is set to zero all the time in the integral on the left and PAGE 50 38 the last integral on the right. In practice, equation (3.33) can be recast as F ( t Af) / P N,F P 7V,G -1 7T rr\ >Ja 1 m N + E d 11=1 Nn c n eÂ„ Â— N -1 E n=l F ( t n) P n,F 1 G(t n )P n , G J (3.34 a,b) where t N > t 0 ; N=1,2,...,N; d Nn can be obtained by referring to (3.31); and (3.35 a,b) In (3.34), the summation vanishes when N-1=0. Then starting from N=1 , 2, and so on, the conditions can be evaluated incrementally. In this effort, the F and G values found from the previous time steps are used as the input in the computation of the right hand side of (3.34), which immediately gives the F and G values at the succeeding new time step. This continues until the desired time is reached. The algorithms can be developed for a rapid solution of the conditions. Computer programs developed for the solution of inverse Stefan problems given in this chapter are provided in Appendices B through E . PAGE 51 39 3.5 Critique of the Method Using GreenÂ’s functions, the present method is limited to the solution of problems whose GreenÂ’s function can be obtained analytically. This excludes those problems whose boundary conditions and governing equations are nonlinear and such nonlinearity can not be resolved by using transformation (e.g., Kirchhoff transformation). Yet with the use of GreenÂ’s function, it projects an impression of being related to the boundary element method that has been undergoing through development in recent years [64]. Yet, they are functionally different. In boundary element method, the boundary element equations are written separately for the liquid and solid regions, whereas in the present method, only one equation (3.27) is derived. The number of the equations to be solved in the present method is thus reduced by half, which is particularly true in the solution by the incremental approach described earlier. More important, in the boundary element method, the problems cannot be solved without the information for the heat fluxes that appear on both sides of the interface [64,1]. Such fluxes, however, are unnecessary in the present work. Instead, equation (3.26) embodies conditions (3.12) and (3.13) in the form of a single integrodif f erent ial equation. There is no need for the satisfaction of the flux conditions; as a result, the present method is more effective because it requires less information for input. As will be shown in the next section, the present method is also accurate. In fact, such accuracy is not unexpected because equation (3.26) is exact. The only approximation in the analysis is PAGE 52 40 the local linearization which has been applied to the interface motion and the boundary condition; see (3.33) and (3.35). In fact, in the present analysis, they have been approximated by using constant elements. The analysis can thus be improved for accuracy by using higher order elements, such as linear and quadratic elements as in references [59,71]. 3 . 6 Results and Discussion For the numerical experiments performed in this study, the interface motion data are taken from reference 5. Aluminum will be used for tests; its properties are given in Table 3.2. The inverse solution techniques developed in this chapter are used to solve eight example problems of which four having exact solutions. The retrieved conditions for these examples can thus be compared with the exact solutions for error. In these four examples, the interfaces move as a function of square root of time, and the interface motion data are used to retrieve the constant temperature conditions that appear on the surface. This problem is known as the Stefan-Neumann problem, and Table 3.2 provides a summary of the conditions tested in the examples. The first example deals with a general two-phase StefanNeumann problem; the medium is initially subcooled to 300 K, which is lower than the melting temperature (see Table 3.3). For this example, the temperature at the boundary is unknown; an unknown temperature is thus assumed to appear at the boundary, and the interface motion data are used to find this temperature. The results are listed in Table 3.4, where two sets of results are PAGE 53 41 Table 3.2 Properties of aluminum Melting temperature, T m (K) 932 Vaporization temperature, T v (K) 2,543 Heat of fusion, L f (J/kg) 389,600 Heat of vaporization, L* (J/kg) 9,462,000 Thermal conductivity, k (W/m K) 200 Density, p (kg/m 3 ) 2,710 Specific heat, c (J/kg K) 1200 PAGE 54 42 Table 3.3 Conditions tested in eight examples PROBLEM DESCRIPTION INPUT DATA TRUE CONDITION F(t):K; G(t):W/m 2 TYPE OF BOUNDARY CONDITION ASSUMED SOLVED 1 two-phase Tj=300K PAGE 55 43 Table 3. 4 Comparison between true and retrieved conditions for first Stefan-Neumann example solved for boundary temperature Boundary Conditions, T (K) Time (sec) Imposed True Condition Retrieved Series Method N=3 N=4 Incremental Method 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1000.00 1003.22 1002.97 1002.72 1002.48 1002.25 1002.03 1001.82 1001.61 1001.42 100123 1001.06 1000.89 1000.73 1000.58 1000.43 1000.30 1000.17 1000.06 999.95 999.85 1003.17 1002.87 1002.58 1002.30 1002.04 1001.79 1001.56 1001.34 1001.14 1000.95 1000.78 1000.63 1000.50 1000.38 100028 100020 1000.14 1000.09 1000.07 1000.07 1003.49 1000.41 100020 1000.11 1000.03 1000.05 1000.04 1000.02 1000.68 999.88 999.97 999.99 999.99 999.99 999.99 999.98 999.98 999.98 999.90 999.94 N F(S) a n Sn ' 1 n m 1 For N=3, a,=1003.490405 a 2 = -0.268732 a 3 = 4.350448xl(T 3 For N=4, a, =1003.490405 aj= -0.322478 a 3 = 6.264645xia 3 aÂ«= 6.552981xia s PAGE 56 44 given--one for the series solution method and the other for the incremental solution method. In the series solution method, a power series of degree N-l is used to represent the temperature (see (3.29b)), and two values of N are tested. In this method, the interface position data at equal time intervals are used for input. Thus, for example, for N=3 , R data at times equal to 0, 6, 12, and 18 seconds are used to determine the series, which is, in turn, used to generate the temperature values listed for 20 time steps in the table. Comparing the retrieved temperatures at the boundary by both series with the true condition (exact temperature) listed to the left indicates they are in good agreement . In this case , the coefficients found for these series are listed at the bottom of the table. From the temperature values tabulated, the results for the low power series appear to be as good as those of the high power series, and such trend persists even with the test of a higher power series of degree 10 (results not shown). This gives the indication that the temperatures have been converged. As for the incremental solution method, the results are also good; errors are of the order of 10 % at large time. For the incremental method, the boundary conditions are found at exactly the same times when the interface positions are given. The time step for the solution is thus identical to that for the interface data input. Also notice that the convergence and stability that are normally encountered in the conventional finite difference methods are nonexistent in the present incremental solution of integral equations. Also as in the case of the series solution method, the incremental solution results are generally PAGE 57 45 better at large time than small time, which will be further discussed later. In the example above, a temperature condition is imposed, and the same Â‘typeÂ’ of condition is assumed in the process of the inverse solution. Since the type of the condition that is imposed on the boundary is unknown a priori, it may well be the heat flux condition that one assumes, and the question to be addressed now is whether it is still possible to retrieve the temperature condition via the assumed flux condition. Use will now be made of equation (3.21) to determine this temperature condition once the boundary flux condition is found, and the results are listed in Table 3.5. As shown in the table, the incremental solution results are still good but the series solution results are not as accurate as those listed in Table 3.4. This is certainly a result of the errors being accumulated first in the evaluation of the heat flux next in the evaluation of the temperature using the previously determined heat flux. While this example serves we ii to illustrate that the conditions are sti 1 1 exchangeable, such a twostep solution of the condition may lead to large errors, part i cu larly in the series solution method, and should thus be avoided in practice. In fact, according to experience, the computer time saved in this two-step approach of solving flux then temperature is insignificant as compared with that in the separate, one-step approach of direct evaluation of the flux and temperature. A slight modification is made in the next two examples: this time the medium is not subcooled, the initial temperature being equal to the melting temperature of the medium (see examples 3 and 4 PAGE 58 46 Table 3. 5 Comparison between true and retrieved conditions for second Stefan-Neumann example solved for boundary flux then temperature Boundary Conditions, T (K) Retrieved Time (sec) Imposed True Series Method Incremental Condition N=4 N=5 Method i 1000.00 585.89 620.85 1023.06 2 1000.00 715.58 758.46 1007.71 3 1000.00 799.99 844.90 1003.60 4 1000.00 862.06 905.50 1000.88 5 1000.00 909.11 948.70 1000.04 6 1000.00 945.07 979.11 999.68 7 1000.00 972.34 999.67 999.52 8 1000.00 992.57 1012.44 999.44 9 1000.00 1006.96 1019.03 998.1 1 10 1000.00 1016.46 1020.76 998.30 11 1000.00 1021.85 1018.76 998.37 12 1000.00 1023.76 1014.02 998.44 13 1000.00 1022.76 1007.48 998.51 14 1000.00 1019.34 999.99 998.59 15 1000.00 1013.95 992.40 998.66 16 1000.00 1007.01 985.51 999.45 17 1000.00 998.89 980.13 999.16 18 1000.00 989.05 977.04 999.11 19 1000.00 980.54 977.04 999.07 20 1000.00 971.00 980.93 997.15 N G(s)=Â£a n s J7-1 71*1 For N=4, a,=8001 187.956002 & 2 = -387201.858263 a 3 = 1581.702473 a 4 = 148.390536 For N=5, a,=8945600.085189 aj= -541131.047558 a 3 = 2763.120508 a 4 = 324.034498 aj= 3.817008 PAGE 59 47 in Table 3.3). The Stefan-Neumann problems solved thus become a one-phase problem. It will be tested that the results are unaffected by the change to the one-phase problem. Again, the same series of tests are made and the results are listed in Tables 3.6 and 3.7. Again, the one-step solution results are good (Table 3.6). The twostep solution results are poorer (Table 3.7), further reinforcing the recommendation made earlier in the testing of the two-phase problems. Having satisfactorily completed testing of the Stefan-Neumann problems, attention is now directed to the solution of inverse Stefan problems whose interface motions must be met by imposing the time-variant temperature and flux conditions. There are no exact solutions for these problems, and the interface motion data are taken from Choi [5] who have solved the regular (forward) version of the problems with great accuracy. The interface position data are then used to retrieve the boundary conditions and the results are listed in Tables 3.8 through 3.11. Tables 3.8 and 3.9 give the results for the direct retrieval of the linear temperature and heat flux conditions F(t) = 1000 + 5t (3.36) G(t) = 6.388034 x 10 6 + 2.82752 x 10 5 t (3.37) while Tables 3.10 and 3.11 give the results for the direct retrieval of the linear temperature and quadratic heat flux conditions PAGE 60 48 Table 3.6 Comparison between true and retrieved conditions for third Stefan-Neumann example solved for boundary temperature Boundary Conditions, T (K) Retrieved Time (sec) Imposed True Series Method Incremental Condition N=3 N=4 Method 1 1000.00 1033.83 1033.12 1037.49 2 1000.00 1030.31 1028.96 998.22 3 1000.00 1026.93 1025.01 998.58 4 1000.00 1023.69 1021.29 998.85 5 1000.00 1020.59 1017.80 999.05 6 1000.00 1017.63 1014.55 999.20 7 1000.00 1014.81 1011.55 999.24 8 1000.00 1012.13 1008.81 999.38 9 1000.00 1009.59 1006.33 999.44 10 1000.00 1007.19 1004.13 999.49 11 1000.00 1004.93 1002.21 999.53 12 1000.00 1002.81 1000.57 999.59 13 1000.00 1000.82 999.23 999.60 14 1000.00 998.98 998.20 999.65 15 1000.00 997.28 997.47 999.65 16 1000.00 995.71 997.07 999.69 17 1000.00 994.29 997.00 999.69 18 1000.00 993.00 997.26 999.73 19 1000.00 991.86 997.87 999.71 20 1000.00 990.85 998.82 999.78 For N=3, a,=1037 .496339 N F ( s) =52 a n s n ' 1 i n ii 3.728416 6.982696xl0" 2 For N=4, a,=1037 .496324 a,= -4.474102 a 3 = 1.005513x1(7' a<= 1.324365xia 3 PAGE 61 49 Table 3 .7 Comparison between true and retrieved conditions for fourth Stef an-Neumann example solved for boundary flux then temperature Boundary Conditions, T (K) Retrieved Time (sec) Imposed True Series Method Incremental Condition N=4 N=5 Method 1 1000.00 906.32 917.43 1056.39 2 1000.00 953.87 966.25 996.60 3 1000.00 980.48 991.90 1002.26 4 1000.00 997.67 1006.84 999.33 5 1000.00 1008.95 1015.10 999.70 6 1000.00 1016.02 1018.67 999.45 7 1000.00 1019.89 1018.86 999.16 8 1000.00 1021.28 1016.61 999.22 9 1000.00 1020.73 1012.66 999.09 10 1000.00 1018.66 1007.67 999.23 11 1000.00 1015.44 1002.18 999.05 12 1000.00 1011.39 996.73 999.32 13 1000.00 1006.79 991.80 998.89 14 1000.00 1001.89 987.87 999.72 15 1000.00 996.94 985.38 998.27 16 1000.00 992.16 984.80 1001.19 17 1000.00 987.77 986.56 994.04 18 1000.00 983.98 991.12 1014.19 19 1000.00 980.99 998.92 951.97 20 1000.00 979.00 1010.44 1162.46 N G(s)=Â£a n s n-l J7*l For N=4, a,=2805963.1 14624 aj= -204888.890714 a,= 2268.882336 99.892255 For N=5, a,=3 137 161 .977 178 a,= -286340.839757 a 3 = 3963.574559 a 4 = 218.130426 a,= 1.467676 PAGE 62 50 Table 3.8 Comparison between true and retrieved conditions for fifth inverse example problem solved for temperature Boundary Conditions, T (K) Retrieved Time Imposed Series Method Incremental Error (sec) true Condition N=3 Error N=4 Error Method % % % i 1005.00 1073.05 6.77 1078.13 7.28 1003.61 -0.13 2 1010.00 1064.27 5.37 1061.16 5.07 1007.83 -0.21 3 1015.00 1056.80 4.11 1048.37 3.29 1012.50 -0.24 4 1020.00 1050.64 3.00 1039.36 1.90 1017.33 -0.26 5 1025.00 1045.81 2.03 1033.71 0.85 1022.27 -0.26 6 1030.00 1042.29 1.19 1031.04 0.10 1027.17 -0.27 7 1035.00 1040.09 0.49 1030.94 -0.39 1032.24 -0.26 8 1040.00 1039.21 -0.07 1033.01 -0.67 1037.26 -0.26 9 1045.00 1039.64 -0.51 1036.86 -0.78 1042.29 -0.25 10 1050.00 1041.40 -0.81 1042.07 -0.75 1047.31 -0.25 11 1055.00 1044.47 -0.99 1048.26 -0.64 1052.35 -0.25 12 1060.00 1048.86 -1.05 1055.02 -0.47 1057.36 -0.24 13 1065.00 1054.56 -0.97 1061.95 -0.29 1062.38 -0.24 14 1070.00 1061.59 -0.78 1068.65 -0.13 1067.38 -0.24 15 1075.00 1069.93 -0.47 1074.72 -0.02 1072.39 -0.24 16 1080.00 1079.58 -0.03 1079.77 -0.02 1077.40 -0.23 17 1085.00 1090.56 0.51 1083.38 -0.15 1082.37 -0.24 18 1090.00 1 102.85 1.17 1085.17 -0.44 1087.50 -0.22 19 1095.00 1116.47 1.96 1084.72 -0.94 1092.06 -0.26 20 1100.00 1131.39 2.85 1081.65 -1.67 1098.80 -0.10 a,=1083. 1636095 %= -10.7640380 a 3 = 0.6587919 a,=1099 .6680063 aj= -23.9563496 a 3 = 2.4866623 34 = -0.0666941 For N=3, F(S) =52 a n sn Â‘ 1 n 1 For N=4, PAGE 63 51 Table 3.9 Comparison between true and retrieved conditions for sixth inverse example problem solved for heat flux Boundary Conditions, q (W/m 2 ) || Retrieved Time Imposed Series Method Incremental Error (sec) irue Method % Condition N=4 Error N=5 Error % % 4.25 7589730.94 8076200.72 6.41 7558675.50 -0.40 7184872.69 -5.33 II 4.50 7660418.94 8147001.48 6.35 7694176.58 0.44 7371675.43 -3.77 4.75 7731106.94 8206716.05 6.15 7815212.59 1.08 7493717.81 -3.07 5.00 7801794.94 8257180.74 5.83 7920261.67 1.52 759645424 -2.63 5.25 7872482.94 8300231.83 5.43 8008672.53 1.73 7689904.72 -2.3 i 5.50 7943170.94 8337705.60 4.96 8080664.49 1.73 7775533.53 2.11 5.75 8013858.94 8371438.35 4.46 8137327.48 1.54 7859641.15 -L92 6.00 8084546.94 8403266.35 3.94 8180622.00 1.19 7940899.93 -1.77 6.25 8155234.94 8435025.90 3.43 8213379.15 0.71 8020420.49 -1.65 6.50 8225922.94 8468553.29 2.94 8239300.64 0.16 8098888.64 -1.54 6.75 8296610.94 8505684.79 2.52 8262958.75 -0.40 8176244.61 -i.45 7.00 8367298.94 8548256.71 2.16 8289796.37 -0.93 8252784.49 -1.36 7.25 8437986.94 8598105.32 1.89 8326126.99 -1.33 8329871.72 -L28 7.50 8508674.94 8657066.91 1.74 8379134.69 -0.61 8402584.80 -1.25 7.75 8579362.94 8726977.77 1.72 8456874.13 -1.42 8480061.68 -1.15 8.00 8650050.94 8809674.19 1.84 8568270.59 -0.94 8553180.08 1.12 8.25 8720738.94 8906992.45 2.13 8723119.93 0.02 8628203.50 -1.06 8.50 8791426.94 9020768.84 2.60 8932088.60 1.60 8701314.04 1.02 8.75 8862114.94 9152839.65 3.28 9206713.65 3.88 8775467.98 -0.98 9.00 8932802.94 9305041.17 4.16 9559402.73 7.01 8847997.32 -0.95 For N=4, a,=3397040.84 N ^=2247925.15 G(S) =^^n S a3=-353 114.747 n 1 a 4 =19587.0623 . For N= =5. a,=8734725.30 a 2 =-3092029.40 83 = 1293821.95 a^-188030.644 aj= 9286.34410 1 PAGE 64 52 Table 3.10 Comparison between true and retrieved conditions for seventh inverse example problem solved for temperature Boundary Conditions, T (K) Retrieved Time Imposed Series Method Incremental Error (sec) true Method % Condition N=3 Error N=4 Error % % 1 972.00 976.35 0.44 1020.27 4.97 969.93 -0.21 2 1012.00 1014.88 028 1044.41 3.20 1006.51 -0.54 3 1052.00 1053.63 0.15 1071.70 1.87 1042.85 -0.87 4 1092.00 1092.58 0.05 1101.84 0.90 1079.43 -1.15 5 1132.00 1131.74 -0.02 1134.51 0.22 1116.09 -1.40 6 1172.00 1171.11 -0.07 1169.40 -0.22 1153.87 -1.55 7 1212.00 1210.69 -0.10 1206.22 -0.47 1192.15 -1.64 8 1252.00 1250.48 -0.12 1244.64 -0.58 1230.89 -1.69 9 1292.00 1290.48 -0.11 1284.37 -0.59 1270.11 -1.69 10 1332.00 1330.70 -0.09 1325.09 -0.51 1310.73 -1.60 11 1372.00 1371.12 -0.06 1366.49 -0.40 1348.84 -1.69 12 1412.00 1411.75 -0.01 1408.27 -0.26 1389.93 -1.56 13 1452.00 1452.59 0.04 1450.12 -0.12 1429.95 -1.52 14 1492.00 1493.64 0.11 1491.73 -0.01 1470.27 -1.46 15 1532.00 1534.91 0.19 1532.80 0.05 1510.54 -1.40 16 1572.00 1576.38 0.27 1573.01 0.06 1550.77 -1.35 17 1612.00 1618.06 0.37 1612.05 0.003 1590.93 -1.31 18 1652.00 1659.95 0.48 1649.63 -0.14 1631.04 -1.27 19 1692.00 1702.06 0.59 1685.42 -0.39 1671.39 -1.22 20 1732.00 1744.37 0.71 1719.12 -0.74 1710.13 -1.26 For N=3, a, =938.0379452 N a?= 38.2157864 F(s)=^a n s n 1 a 3= 0.1050509 /j* 1 For N=4, a,=999.5873340 %= 18.8517864 a 3 = 1.8834608 a 4 = -0.0513597 PAGE 65 53 Table 3.11 Comparison between true and retrieved conditions for eight inverse example problem solved for heat flux Boundary Conditions, q (W/m 2 ) Retrieved Time Imposed Series Method Incremental Error (sec) i rue Condition N=4 Error N=5 Error Method % % % 0.25 3003125.00 3007778.76 0.15 3033589.06 1.01 2995616.25 -0.25 0.50 3012500.00 3021875.02 0.31 3035992.94 0.78 3004320.50 -0.27 0.75 3028125.00 3040169.58 0.40 3043865.79 0.52 3019590.53 -0.28 1.00 3050000.00 3062934.56 0.42 3057990.24 0.26 3036905.38 -0.43 1.25 3078125.00 3090442.06 0.40 3078999.09 0.02 3061522.20 -0.54 1.50 3112500.00 3122964.20 0.34 3107375.31 -0.16 3092582.42 -0.64 1.75 3153125.00 3160773.08 0.24 3143452.06 -0.31 3130006.55 -0.73 2.00 3200000.00 3204140.84 0.13 3187412.63 -0.39 3173845.90 -0.82 2.25 3253125.00 3253339.57 0.01 3239290.49 -0.43 3223966.45 -0.90 2.50 3312500.00 3308641.38 0.12 3298969.30 -0.41 3280396.20 -0.97 2.75 3378125.00 3370318.40 -0.23 3366182.86 -0.35 3343088.04 -1.04 3.00 3450000.00 3438642.74 -0.33 3440515.14 -0.27 3411906.45 1.10 3.25 3528125.00 3513886.50 -0.40 3521400.30 -0.19 3487162.49 -1.16 3.50 3612500.00 3596321.81 -0.45 3608122.64 0.12 3568166.70 -1.23 3.75 3703125.00 3686220.76 -0.46 3699816.64 -0.08 3656322.04 -1.26 4.00 3800000.00 3783855.48 -0.42 3795466.96 0.12 3748732.19 -1.35 4.25 3903125.00 3889498.09 -0.35 3893908.40 -0.23 3851068.35 -1.33 4.50 4012500.00 4003420.68 0.22 3993825.95 -0.46 3951711.24 -1.51 4.75 4128125.00 4125895.38 -0.05 4093754.75 -0.83 4084808.11 -1.05 5.00 4250000.00 4257194.29 0.16 4192080.13 -1.36 4096803.88 -3.60 N G(s)=Â£aÂ„s n-1 71*1 For N=4, a, =2997608.684 aj=33 190.75084 a 3 =29232.59062 a 4 =2902.536744 For N=5, a, =303572 1.67 a^-16210.2870 aj= 27733.4034 a 4 = 12343.6674 35 =1598.2 1379 PAGE 66 54 F(t) = 932 + 40t (3.38) G(t) = 3xl0 6 + 5xl0 4 < 2 (3.39) In these tables, the computational errors are calculated using the following definition: respectively . As described in Table 3.3, those in Tables 3.8 and 3.11 are for the medium initially at phase change temperature, while those in Tables 3.9 and 3.10 are for the medium initially subcooled at 300 K . The former are thus one phase problems while the latter are twophase problems. Again both series and incremental solution methods are used for solution and their results are good. For example in Table 3.8, the series solution results converge even with a value of N that is as low as 3, whereas in seeking the flux condition in Table 3.9, the series converges rapidly from N=4 to N=5 (higher degree results not shown). In Tables 3.8 and 3.11 the medium melts as soon as the boundary conditions are imposed, whereas in Tables 3.9 the medium starts to melt at time greater than 4 seconds. In all cases, the accuracy of the results appears to be unaffected by the time when melting takes place. The results in Tables 3.10 and 3.11 follow the same trends as the ones in Table 3.8 and 3.9. Tests Error=l Â—2 (3.40) where p and q represent retrieved and imposed true conditions, for the time variant conditions are thus successful. PAGE 67 55 The inverse solution techniques developed in this study are expected to be accurate as mentioned earlier that is close to the end of the previous section. Yet for the Stefan problems solved in this work, the accuracy of the techniques is better at large time than small time. This can be attributed to the curvature of the interface position curve, which is always large at small time [2]; see Figure 3.2. Then, in the numerical solution of (3.27), there will be a slight error associated with the linearization of the position at small time. At large time, however, the position curve tends to be linear; the linearization error will be diminished. In fact, as time progresses, the accurate terms under the summation in (3.27) rapidly outnumber the inaccurate terms to the effect that the boundary conditions can always be evaluated accurately with the present method at large time; see the results in Tables 3.4 through 3.11. This is a distinct departure from the trends of other time marching schemes reported in the literature in which the errors tend to grow with time. It should also be pointed out that, for the Stef an-Neumann problem chosen for comparison in the present study, there is a singularity of the temperature at zero time. This also contributes to the large discrepancy of the results at small time, which must not be overlooked. It has been firmly established that, in the solution of the Stefan problems, only the Stefan-Neumann problems can be solved exactly. Yet, it is also possible to develop an exact solution for an exponential condition imposed on the boundary; such condition, however, has been considered as physically untenable in the literature [7]. Worse yet, such a condition gives rise to a PAGE 68 56 constant velocity of the interface, a situation making the present linearization scheme exact in the solution of the inverse problems [5]. Perfect results will be obtained, rendering the test of the exponential conditions meaningless. 3 . 7 Extension and Concluding Remarks The analysis developed in this chapter can be readily extended for the solution of inverse problems with multiple phases. For such problems, times for re-melt and re-freeze of the medium must be closely accounted for, and the analyses [4,65] can be readily adapted for the development of the inverse solution techniques presented in this study. The present analysis can also be extended for the solution of inverse problems in which the properties are unequal for different phases of the medium. For such problems, double source and sink fronts must be used as given in the solution of the regular problems in reference 17. Finally, it is noted that although problems in one-dimensional, semi-infinite domain have been solved for examples in this work, problems in finite domains (e.g., plane wall) can also be solved with the present methods. For these problems, there are two boundaries and two boundary conditions are imposed. Two unknowns are thus sought simultaneously, and this requires the input of one additional condition in the form of either temperature or heat flux at any interior point close to the boundary where no phase change takes place. On the other hand, for problems with two phase change interfaces caused by separate heat input simultaneously from both sides of the boundaries, the interface motion data for the second interface will serve as this additional PAGE 69 57 condition . In any event, no flux information is needed at both sides of the interfaces. Furthermore, the present method can also be applied to the solution of problems in multiple dimensions. Again the inverse solution for these problems can be developed on the basis of the solution of regular versions of these problems. PAGE 70 CHAPTER IV SOLUTION OF ABLATION PROBLEMS WITH ONE MOVING BOUNDARY BY A SOURCE-AND-SINK METHOD It is the purpose of this chapter to present a source-and-s ink method for the solution of ablation problems with one moving boundary. As will be shown, the essential feature of this method is that the solution will be sought in a fixed domain which does not change with time. In fact, the ablated region is to be treated as a fictitious domain in which the flux condition at the ablated boundary will match that of the imposed condition at the moving boundary. Then, with the additional vaporization temperature given for the moving boundary, the conditions at this boundary are overspecified for the fictitious domain. The problem can thus be taken as an inverse problem in the sense that the condition on the fixed boundary is sought during the ablation period. This condition will provide for the necessary conditions at the moving boundary for the solution of the problem. 4.1 Solution Methodology For a subcooled medium, the problem can be divided into two stages; namely, pre-ablation stage and ablation stage. During the pre-ablation stage, heat is added to the surface of the medium in order to raise its temperature to the phase-change temperature. Then with continuous heating, ablation takes place. The ablated 58 PAGE 71 59 surface moves inward with time and the imposed heat flux follows this boundary motion. As shown in Figure 4.1, it is assumed that the medium is homogeneous and isotropic. The thermophysical properties are constant. For the ablation problem, the ablated region is immediately removed upon formation. Radiation is taken to be a surface phenomenon. It is also assumed that the medium changes phase at a distinct temperature; that is, no mushy zone in the medium. Moreover, there is no volumetric heat generation in the medium. Under these assumptions, the ablation problem can be formulated as follows: Pre-ablation Stage Governing equation-V 2 T ,. t x_iÂ«W) * l 0 ( r Â»t) Â— a f G (A US) t 0 > t > 0 (4.1) Initial condition-T o (r,0) = T,(r) (4-2) Boundary condition-3T 0 (f 0 ,t) dn t Â’ (4.3) PAGE 72 Figure 4.1 System analyzed PAGE 73 61 PAGE 74 62 Ablation Stage Solid Region: Governing equation-V 2 T 5 (r,t) _ 1 ST s (r,t) Q dt r Â€ S t > t 0 ( 4 . 4 ) Initial condition-T s ( r ,t 0 ) Â— T 0 (r,to) ( 4 . 5 ) Boundary conditions-T s (r s ,t) = tÂ„ ( 4 . 6 ) G(t) = k + pLvSl ( 4 . 7 ) where all notations have their usual meaning; see nomenclature. Here r PAGE 75 63 The ablation problem above is solved by using a source-andsink method. In this method, the domain for investigation is extended to cover AUS; the ablation front becomes an interior phasechange interface. Solution is thus rendered in a fixed domain, in which the ablated region is fictitious in the sense that it is physically nonexistent; however, the unablated (solid) region is real. The conditions imposed at the moving boundary in the original ablation problem are taken to be those imposed on the interior phase-change interface. Thus, for the new problem in the fixed domain, the interface position and the fictitious condition on the fixed boundary are the unknowns to be determined. Once they are found, they can be used for input in the determination of the temperature in the solid region. The problem is thus solved totally. An equivalent problem is formulated with a source-and-s ink method as follows: Equivalent Problem Governing equation-V 2 T(f,t) pL Â— r/ Â— Â— \ 1 v.nÂ«5(r r 5 ) = s dT(r,t) dt FG(AuS) (4.8) t > t 0 Initial condition-T(r, t 0 ) = T 0 (r,t 0 ) (4.9) Interface conditions-- PAGE 76 64 T(? 5 ,t) = TÂ„ (4.10) dT(fg,t) <9n + pLv.n (4.11) where 8(rÂ—r s ) denotes a Dirac delta function. As has been shown in Chapter 3, equation (4.8) reduces to (4.4). Also integrating (4.8) across the interface at r^ enables the G(t) to represent the heat flux imposed on the interface at the fictitious side. Equation (4.8) can be solved by using (3.18) in which the domains for integration for the first and second terms on the right are changed to AUS. Also minus sign is taken for the Dirac delta function term in which fy is changed to r s as T(r,t)= G(f,t | r',t 0 )T .(f')dV -/\ A U 5 t G(r,t | r',r)v.n<5(r' Â— r^)dV'dr + | t Q Au5 1 (4.12) This temperature equation is forced to satisfy the interface temperature equation (4.10) by setting r to the ablated surface position, r 5 , and T(r 5 ,t) to the phase-change temperature, T^, as t V G(r 5 ,t | r / ,t 0 )T-(r / )dV / Â— ^ G(r 5 ,t | r',r)v.7t6(r' r 5 )dV'dr + ^( j (4.13) A U 5 t Q ,4 U 5 Equation (4.12) is also differentiated to satisfy (4.11) as PAGE 77 65 G(t)= k< Â™^T j(f Oav' t L c 'A u S % 3G(r 5 ,t| Â•/ % dn -i -v.n6( v' Â— r^)dV / dr t 0 AuS + }_ _ } i r=r s ) + pLv.n (4.14) where differentiation is to be effected at the solid side of the interface. Notice that the boundary condition is embodied in the summation term, which, according to Table 3.1, should be expressed in the form Â»= a r( t ,-/ r ^T(r 0 ,t) G(r,t I r 0) ^) Qn arpor (4.15) o r o where the following condition is taken at the fixed boundary: dT(rp,t) _ G(fQ,t) dn : k for t < t 0 and Â— ^ for t > t Q (4.16) where g(rg,t) represents the imposed fictitious heat flux. Clearly, there are two unknowns in (4.13) and (4.14): g(Fg,t) and v; there are two equations to solve them. In this effort, t 0 and T 0 (r,* 0 ) are to be found from the solution of a pre-ablation problem as described previously . 4.2 Illustrative Examples Use is now made of the solution methodology given in the preceding section to solve examples in semi-infinite medium. Four examples will be provided and they include medium with or without subcooling with the boundary imposed with constant or variable heat flux conditions. PAGE 78 66 For a subcooled medium, the pre-ablation stage solution has been given in Chapter 3. The pre-ablation temperature can be taken from (3.22) as T 0 (x,t) = T, + i N a 7 r t 0 G(s) (t-s) 1/2 exp[4a(t Â— s) ]ch (4.17) Thus the time when ablation starts (t 0 ) can be found by solving the following equation implicitly T = T-h1 1 ; 1 * ' V kN a 7 r (4.18) Also at the moment when ablation starts, the temperature in the solid region is given by the relation To( x ?lo) Â” o r a 7 r 0 G(s) 1/2 exp[X 4a(< 0 -s) ]ds (4.19) which serves as the initial condition for the ablation stage. For the problem at hand, (4.12) can be used to derive temperature as J 0 T(x,t) = T,+ i a * 1 T Â— 0 t exp[ Â—77 v ]dr -f r (t -r ) 1/2 4a(t Â— r) k t a g(' r ) kNir 172 ex P[ X T Â— t (t -r ) 1/2 4a(t Â— r) ]d0 L c TÂ—t d#i( r ) ^J47ra(t Â— r) 1 |exp( ft Â— t) ^ (x-^i(r)) : 4o(t Â— r) ) + exp( (x+ii 1 (r)) 4a(t Â— r) 0 (4.20) Also (4.13) is used to derive PAGE 79 67 t o 1 T = T + A v ^ v a G(r) ex p[R i(t) kN7r J (t-r) 1/2 4a(t Â— r) r=o v Â’ ]dr + r a 1 k N7T g( r ) r u_ exp[ Â— T777 Z\ J dr TÂ— t 0 t L c r= t dT ^47ra(t Â— r) L_ {exp( ft Â— t) ^ (^(t)-^)): 4a(t Â— r) Â•) + exp( 4a(t Â— r) (^(tj+^r)) 4a(t Â— r) 0 (4.21) Equation (4.14) follows as t o G(t) = G(r) r R\{ t) nj , R,( t) ex P[ ~ T7Tl ~ l dr + 2^f!r fil(t) 2>f7ra J ( t _ r \3/2 4a(t Â— r) r=o v ' t g( r ) r u_ ex P[ Â— J dr T Â— t (t-r) 3/2 4a(t Â— t) 0 Lp 4>j7ra d ^i ( r ) dr r=t * 3/2 {(fli(t) ^ 1 (r))exp( (t-r) 4a(t Â— r) ) 0 + (fl 1 (t)+ii 1 (r))exp((R,(t)+R,(r)) 4a(t Â— r) + * dt (4.22) Solution of (4.21) and (4.22) subject to the initial condition, f?i(t 0 )=0, yields the results for two unknown functions i?j(t) and SO)4.3 Numerical Solution of Ablated Front Equations (4.21) and (4.22) are coupled nonlinear integrodif f erential equations, which will be solved numerically. In this effort, a local linearization scheme is employed as described in Section 3.4. For the present problem, the entire time range is dR,(t) and g(t) are divided into small increments in which both dt treated as constants. They can thus be taken out of their respective integrals, and the convolution integrals in these equations are changed to summations as PAGE 80 68 t o 1 T V = T t + fc a G(r) exp[k ri 0 (tjV -^ /2 4 a (^-r) ]d r t _ N -f Â£Â«Jf E g(tÂ„) K n=l J TÂ—t 1 r ^1 ( t 7 v) :j ' xp[ 53(i^) ]dr n1 _L gdiii(tÂ„) n /* C ^ nÂ— l ai . ^ 47 ra(t N -r) r=t n-l (^l(t n) ~ ^i( r )) 4 a(t n -t) ) + exp( (fi^t N )+R 1 (r))\ \ A _ Aa(t N -r) >!* (4.23) and G(t j\f) Â— ^l(tyv) { J 0 GM Â„ p[ _agg) |dT 2 ^ U (t w -r ) 3/2 4 o (t N~ t ) TÂ— 0 V N 7 t N + E g(t n ) n=i r=t n 1 exp[ R i(t n) ( t N -rf /2 n-1 t Lp n 4 MWa n ^ dt J (tjv _ r) 1 f/p /, \ n / \\ / (Â•^l( t 7v) _ -^l( r ))\ 73372 4a(tjv _ r ) Â— ) n-l , ,r> u \ . n / \\ t (^l( t Af) + -^l( T )) 2 \\j_ , + (Ri(tN)+ R i( T )) e M 4a ( t N -7j Â— ) / tlr + pL dt (4.24) Notice that the removal of diJ 1 (t) dt and g(t) will cause a slight error in the numerical solution of the ablated front position. However, R 1 is a gradual function of t, and so is g(t) during the early stage of ablation. The error associated with the linearization is expected to be small. In practice, such an error can always be reduced by taking small time increments. PAGE 81 69 It is also noted that the . ablated position iZj(t) in the GreenÂ’s function can be related to by means of the interface velocity, Â’ as s ^ own figure 3.2. The position R a at any time t should thus be discounted as unknown. Numerically, the R 1 ( t) and g(t) in these equations will be solved incrementally starting from t 0 until the final time tjy is reached. This corresponds to a time marching scheme in the numerical solution. Computer programs developed for the solution of the ablation problems given in this chapter are compiled in Appendix F. 4.4 Results and Discussion Four examples are tested and they include semi inf ini te medium with and without subcooling with the moving boundary imposed with constant and time-variant conditions (see Table 4.1 for summary). Notice that, in this table, under the column of problem description, the ablated region is always counted as one phase. Thus, for a medium initially at phase-change temperature imposed with a constant heat-flux condition, such as Example 1, there is one ablated phase. It is thus called one -phase problem in the table. Example 1 deals with an ablation problem that can be solved exactly. As shown in Carslaw and Jaeger [7], for a semi inf ini te medium ablated with a constant velocity U, the heat flux imposed on the boundary is given by the relation G = [L + c(T ( ,-T i )]pU (4-25) where T t is constant and the velocity is related to the ablation position as PAGE 82 70 Table 4.1 Conditions tested in four examples Problem Description Material Heat Flux Condition Imposed G(t) (W/m 2 );! (S) Remarks 1. One-phase T=932 K=T m =T v Aluminum G(t)=5xl0 5 Exact solution is available for this case. 2. Two-phase TÂ—882 K PAGE 83 71 U _ (4.26) It is expected that, for a medium initially at the phase-change temperature, G is equal to /?LU product. It is also noted that equation (4.25) is strictly valid for the medium in a quasi-steady state. At that time, the temperature in the solid is stabilized given as T(x,t) = T v e ~ ~ Ut ) (4.27) which occurs for the medium having ablated for a long period of time. The first two examples are thus chosen to substantiate these points . For the first example, the ablation is due to melting. The medium is initially at the melting temperature so that the surface ablates as soon as the heat is applied. The position of the ablated surface is shown in Figure 4.2, where the curve appears to be linear. Using the exact solution (4.25) for comparison yields the error curves shown in Figure 4.3. Clearly, the results are stable and converge satisfactorily, errors being less than one tenth of one percent for all the time increments tested. Also the errors are slightly larger at small time, a result of the singularity associated with the heat flux that is abruptly appl ied on the surface at time zero . In what follows in the numerical computation of all the examples to be presented in this chapter, the time increment is taken as 0.5 sec, unless otherwise noted. Example 1 provides an excellent test for the algorithms PAGE 84 0) P c 3 0 P aj P o3 Sh * rH 0) ~C U a c 0 E 0 P a) u p c X 0 0 3 Â• pH bO f-H P G P Â• pH o3 i to J3 P 0 0 o3 a l 0) PAGE 85 73 o GO a CO A CM o # w o w o w o o o (ujo) Â‘uoi^isod ooDj-ins p9}D|qy Time, t (s) PAGE 86 U U P P p U cfi O O In 0) > p U U Bj (1) CH box 0) P P P b0 cb cb r * H U P! P 0) P U I > O I p c U O P Cfl U U O P P5 Â•rH yJZj -o p a p cp c flj I-H p flj O b3 P >>0) co p >Â» p Â• rH > c O -P 0 C c$ 0 D rÂ“H bOX> 03 P 0) -fi > Sh P> fi 0 Â• pH 0 S u c 0 T3 Â•"O Â»ri 0) C 4-> (0 O H a >> 0 E P CO Â• pH iÂ“H Q) E Â•
> 0 E S3 P CO Â• pH 0 Â• pH Â• pH r-H 0 E P Â• r-l P *-H rD P Â«$ HD P P fi CD 0 CD **h E U Oi P P bO Â• pH Uh cb E CO *Â»-i 0) p 'POE S-< G G >> Â•Â•* O U G T3 ,H G O O 'P Â« > P O O P P Â•* Â• iÂ—* fÂ— H Â£ P p-h O O Â•fH (/} I-H P X! X T3 C$ O O C$ P X P 3 cx p a cr 05 X o p hO Â• pH &U |