Session 5
PHASE CHANGE I
Chairman: T. Kowalewski
TRANSPORT PROCESSES IN POLYMER EXTRUSION AND OPTICAL FIBER DRAWING
Yogesh Jaluria Department of Mechanical and Aerospace Engineering Rutgers, the State University of New Jersey New Brunswick, NJ 08903, USA
The numerical modeling of the transport processes in materials processing, particularly for polymers
and glass, is discussed in this presentation. Schematics for two important manufacturing processes
concerned with these materials, the optical glass fiber drawing and the plastic screw extrusion
processes, are shown in Fig.1. Heat transfer and fluid flow considerations play a very important role
in determining the properties of the final product. In this paper, the governing equations, the relevant
boundary conditions, and the basic approach to numerical study of these processes are discussed.
The importance of material characterization in an accurate mathematical and numerical modeling is
brought out, along with the coupling between the process and the resulting properties. The various
concerns and computational difficulties that arise due to the nonNewtonian behavior of plastics,
exponential temperature dependence of glass viscosity, complex geometry, free surface flow,
conjugate transport, and many other such effects encountered in practical processes are considered.
The material undergoing the thermal transport process is moving, with the velocity level changing
dramatically for optical fiber drawing due to the large change in diameter. Additional mechanisms
such as surface tension effects and chemical reactions are also important. Some of the important
methods to treat these problems and challenges are presented. A few characteristic numerical results
are discussed and validated by comparisons with experimental results. The implications of such
results in improving practical systems is also considered.
Polymer Extrusion
The variation of viscosity , requires special consideration for polymers, such as plastics, rubbers,
and food materials, that are of interest in a variety of manufacturing processes. Most of these
materials are nonNewtonian in behavior, implying that the shear stress is not proportional to the
shear rate. The viscosity is a function of the shear rate and, therefore, of the velocity field. The
viscosity is independent of the shear rate for Newtonian fluids like air and water, but increases or decreases
with the shear rate for shear thickening or thinning fluids, respectively. These are
viscoinelastic (purely viscous) fluids, which may be timedependent or timeindependent, the shear
rate being a function of both the magnitude and the duration of shear in the former case. Viscoelastic
fluids show partial elastic recovery on the removal of a deforming shear stress. Food materials are
often viscoelastic in nature.
Various models are employed to represent the viscous or rheological behavior of polymers of
practical interest. For instance, timeindependent viscoinelastic fluids without a yield stress are often
represented by the powerlaw model, given by^{l}
where K_{c} is the consistency index and n the power law fluid index. Note that n = 1 represents a
Newtonian fluid. For n < l, the behavior is pseudoplastic (shear thinning) and for n > 1, it is dilatant
(shear thickening). The viscosity variation may be written as^{l, 2}
for a twodimensional flow, with the velocity components u and w varying only with coordinate y.
Similarly, expressions for other two and threedimensional flows may be written. Here is the shear
strain rate, is the shear stress, the subscript o denotes reference conditions and b is the temperature
coefficient of viscosity.
The fluid viscosity is very large, being more than a million times that of water. This usually gives
rise to very small Reynolds numbers, for which the creeping flow approximation is often employed.
For instance, the Reynolds number Re is generally smaller than 1.0 for plastic flow in a single screw
extruder and the inertia terms are usually dropped. The coordinate system is generally fixed to the
rotating screw and the channel straightened out mathematically, ignoring the effects of curvature.
Then the complicated flow in the extruder is replaced by a pressure and shear driven channel flow,
with shear arising due to the barrel moving at the pitch angle over a stationary screw. This
substantially simplifies the numerical model^{2}. Figure 2 shows the typical computed temperature
distribution in an extruder channel. Large temperature differences are seen to arise across the
channel height because of the relatively small thermal conductivity of plastics. The flow is welllayered,
with little bulk mixing, due to the high viscosity of these fluids. Viscous dissipation causes
the temperature to rise beyond the imposed barrel temperature^{2}.
Optical Fiber Drawing
Glass is another very important, though complicated, material. It is a supercooled liquid at room
temperature. The viscosity varies very strongly with temperature. In optical fiber drawing, for
instance, the viscosity changes through several orders of magnitude in a relatively short distance.
This makes it necessary to employ very fine grids and substantial underrelaxation to achieve
convergence. Even a change of a few degrees in temperature in the vicinity of the softening point,
which is around 1600 ^{o}C for fused silica, can cause substantial changes in viscosity and thus in the
flow field and the neckdown profile in optical fiber drawing^{3, 4}. At the free surface, between the
glass and the inert environment, the shear stress is often specified as zero. If the shear stress exerted
by the ambient fluid is significant, it is taken into account. Basically, a balance of all the forces
acting at the surface is used to obtain the interface. The free surface may be determined numerically
by iterating from an initial profile and using the imbalance of the forces for correcting the profile at
intermediate steps^{4}. Figure 3 shows such an iterative convergence for the neckdown profile of the
glass preform in the optical fiber drawing process.
At its softening point, the viscosity is still very high, being of the same order as that of polymer
melts. Thus, viscous dissipation is important and the energy and momentum equations are coupled.
However, glass flow may be treated as Newtonian. In optical fiber drawing, the diameter of the
cylindrical rod, known as a preform, typically changes from 25 cm to about 125 m in a distance of
only a few centimeters. This places stringent demands on the grid as well as on the numerical
scheme because of the large change in the surface velocity. The flow in the neckdown region is
smooth and welllayered because of the high viscosity. A typical temperature difference of 50100
^{o}C arises across the fiber. Even this small difference is an important factor in fiber quality and
characteristics. Viscous dissipation, though relatively small, is mainly concentrated near the end of
the neckdown, in the small diameter region, and plays an important role in maintaining the
temperatures above the softening point.
REFERNCES
 Tadmor, Z. and Gogos, C., Principles of Polymer Processing, Wiley, New York,1979.
 Jaluria, Y., Heat and Mass Transfer in the Extrusion of NonNewtonian Materials, Advances in Heat Transfer, Vol. 28, pp.145230,1996
 Lee, S.H.K. and Jaluria, Y., Effect of Variable Properties and Viscous Dissipation During Optical Fiber Drawing, ASME J. Heat Transfer, Vol. 118, pp. 350358, 1996. Also, Int. J. Heat Mass Transfer, Vol. 40, pp. 843856,1996.
 Roy Choudhury, S., Jaluria, Y. and Lee, S.H.K., Generation of NeckDown Profile for Furnace Drawing of Optical Fiber, ASME Heat Transfer Div., Vol. 306, pp. 2332,1995.
EFFECT OF BUOYANCYDRIVEN CONVECTION ON MELTING WITHIN SPHERES
Y. Zhang*, J. M. Khodadadi* and F. Kavi ka+
*Department of Mechanical Engineering, Auburn University, 201 Ross Hall
Auburn, AL 368495341, USA
+Technical University of Brno, College of Mechanical Engineering, Technická 2, 616 69 Brno, CZECH REPUBLIC
INTRODUCTION
Melting and solidification of materials are important processes which are encountered in nature and numerous industrial applications. Due to the recent upsurge of interest in such processes as spray casting, latentheatoffusion thermal storage systems, space radiators, spacecraft energy conversion systems, rapid solidification, purification of materials, containerless levitation technique, etc., there is an urgent need to analyze the transport phenomena of melting and solidification in spherical configurations. Only a few authors have studied the transport phenomena with phase change in a spherical shape. Tao1 presented some numerical results for the diffusioncontrolled solidification process inside spheres with a highly idealized boundary condition. Few years later, Riley et al.2 found an approximate solution for this particular problem by employing the perturbation method. Pedroso and Domoto3 and Gupta4 also obtained some results for the diffusioncontrolled inward solidification problem in the spherical geometry, all with idealized boundary conditions. Gammon and Howarth5, still employing the perturbation method, obtained some results for the inward solidification problems with a little more complicated boundary conditions. Moore and Bayazitoglu6 were the first authors to present some experimental results for the melting process within a spherical enclosure, where the convective effects were first mentioned.
Transport phenomena can play a significant role during the processes of melting and solidification. For example, consider the latentheatoffusion thermal storage system. In highpower space systems, large amount of waste heat will be generated. As the heat loads are usually timedependent, heat sinks and heat sources have to be used. In the inputing and extracting processes of thermal storage operations, it has been reported that natural convection plays a significant role in affecting the heat inputing and extracting rate for cylinder containers (Ishikawa et al.7, Saitoh et al.8, Saitoh and Kato9), as the high thermal storage density Phase Change Materials (PCMs) usually have unacceptably low thermal conductivities. Experimentally, the significance of natural convection in the operation processes of thermal storage system have been confirmed in spherical configurations (Saitoh and Moon10). However, no analytical investigations have been conducted for the spherical configurations to the best knowledge of the authors. Therefore, a thorough investigation of transport phenomena during the melting process inside spherical configurations is in need.
In addition, when convection is considered in phase change problems, all the research papers published were conducted with very idealized boundary conditions, and most of them are conducted in the Cartesian coordinate system, actually in rectangular cavity. In order to understand the significance of natural convection in such practical fields as spray casting, rapid solidification, and quantitative information in the latentheatoffusion storage system, a thorough study of melting within spheres is needed.
METHODOLOGY
The present work will utilize the SemiImplicit Method for PressureLinked Equations (SIMPLE) procedure of Patankar11 to solve the unsteady coupled fluid flow and heat transfer equations governing the melting process. The fixedgrid method which utilizes single formulations for different phases is used and the coupling condition at the phase change interface is formulated into the governing equations. Thus, there is no need to track the moving interface and a fixed grid can be used throughout the computations. This approach is wellsuited for treating the continuous transition between solid and liquid phases, as well as the evolution of latent heat over a finite temperature range. An enthalpy formulation which can handle phasechange problems occurring both at a single temperature and over a temperature range will be utilized.
RESULTS AND DISCUSSION
A parametric study to assess the role of buoyancydriven convection during melting within spherical containers has been undertaken. Emphasis is placed on lowPrandtl number liquid metals undergoing melting. The surface temperature of a solid spherical metal specimen is raised to a value greater than its melting temperature. Melting will initiate at the surface and the solidliquid interface will move into the specimen. As the extent of the liquid state is increased, buoyancydriven convection will be promoted and a skewed solidliquid interface will be observed. A typical plot of the streamlines and temperature contours at a given instant (dimensionless time = t* = 0.095) are given in Figure 1 for the Rayleigh and Stephen numbers of 1.2x105 and 0.4, respectively. Clearly, the temperature contours shown on the right half of Figure 1 are skewed and different than concentric contours expected for pure diffusioncontrolled melting. This is due to the active role of buoyancydriven convection as evidenced on the left half of Figure 1. Also note that at this particular time instant, a second recirculation zone is observed at the top portion of the sphere.
More detailed unsteady flow field information including velocity vectors, migration of the interface and temperature contours within the sphere will be plotted and discussed in the paper. Results will be compared to the limiting case of diffusioncontrolled melting.
REFERENCES
 Tao, L. C., Generalized Numerical Solutions of Freezing a Saturated Liquid in Cylinders and Spheres, AIChE Journal, Vol. 13, No. 1, pp. 165169, 1967.
 Riley, D. S., Smith, F. T. and Poots, G., The Inward Solidification of Spheres and Circular Cylinders, Int. J. Heat Mass Transfer, Vol. 17, pp. 15071516, 1974.
 Pedroso, R. I. and Domoto, G. A., Inward Spherical SolidificationSolution by the Method of Strained Coordinates, Int. J. Heat Mass Transfer, Vol. 16, pp. 10371043, 1973.
 Gupta, S. C., Analytical and Numerical Solutions of Radially Symmetric Inward Solidification Problems in Spherical Geometry, Int. J. Heat Mass Transfer, Vol. 30, No. 12, pp. 26112616, 1987.
 Gammon, J., and Howarth, J. A., The Inward Solidification of Spheres with a Slight Perturbed Temperature Distribution at the Boundary, Int. Comm. Heat Mass Transfer, Vol. 23, No.3, pp. 397406, 1996.
 Moore, F. E. and Bayazitoglu, Y., Melting within a Spherical Enclosure, Journal of Heat Transfer, Vol. 104, pp. 1923, 1982.
 Ishikawa, M., Hirata, T. and Tamaki, H., Melting and Solidification Phenomena of Two Insoluble Phase Change Materials In a Cylindrical Capsule, Trans. JSME, Vol. 57 B, pp. 41744181, 1991.
 Saitoh, T., Kato, H. and Maruhara, K., Numerical Computation on Combined Natural Convection and CloseContact Melting in a Cylindrical Energy Storage Capsule, Proc. Int. Conf. Comput. Sci., pp. 15, 1992.
 Saitoh, T. and Kato, H., Experiment on Melting in Heat Storage Capsule with Close Contact and Natural Convection, Exp. Therm. Fluid Sci., Vol. 6, pp. 273281, 1993.
 Saitoh, T. and Moon, J. H., Experimental Performance of Latent Heat Thermal Storage Unit Packed with Spherical Capsule, Proc. of the 5th Int. Energy Conf., pp. 18, 1993.
 Patankar, S. V., Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, NY, 1980.
ACKNOWLEDGEMENT
The last two authors acknowledge the support provided by the USCzech Republic Science and Technology Program
Figure 1. Streamlines and temperature contours at t^{*} = 0.095
A NEW HEAT TRANSFER MODEL FOR PREDICTING ICE
ACCRETION ON AIRFOILS
Orhan Kural* and Tuncer Cebeci*
*Department of Aerospace Engineering
California State University, Long Beach
Long Beach, California 90840, USA
Aircraft icing presents a serious hazard for flight at subsonic speeds in visible moisture and at temperatures near or below freezing. Many aircraft have been lost due to ice accumulation. In the absence of thermal ice protection, ice on wings, control surfaces and engine intakes can reduce the aerodynamic performance of the aircraft. Therefore, the Federal Aviation Administration (FAA) requires an airplane manufacturer to demonstrate that its aircraft can fly safely in icing conditions as defined by the socalled icing envelopes in the FAA’s Federal Airworthiness Regulations. For this reason, it is necessary to predict ice accretion on lifting bodies and be able to determine the aircraft performance degradation due to icing.
Currently the most widely used method in the United States for predicting ice accretion on airfoils is LEWICE, a twodimensional ice accretion computer program which was initially developed at the University of Dayton Research Institute and later modified by Ruff and Berkowitz^{1}. This program provides a basis for the determination of the ice buildup on the leading edges of airfoils by specifying the atmospheric parameters of temperature and pressure, the flight velocity and the meteorological parameters of liquid water content (LWC), droplet diameter, and relative humidity. The three major elements of this code are: (1) the flowfield calculation, (2) the particle trajectory calculation, and (3) the energy balance. While the first two elements of this code are relatively rigorous, the third element is not. This is (1) partly because the local heat transfer coefficients used to calculate the convective heat transfer are computed by empirical formulas based on Reynolds analogy and (2) partly because the physical model for the ice accretion process is limited only to rime ice and is not applicable to glaze ice.
This paper addresses the above shortcomings by developing an accurate method for computing the heat transfer coefficients for both laminar and turbulent flows by using boundary layer theory and describes work in progress for developing a physical model for the ice accretion process.
To improve the accuracy of predicting the heat transfer coefficient, the conservation equations for mass, momentum and energy are solved. In this way the adhoc formulas based on Reynold’s analogy are replaced. A finitedifference method is used to solve the conservation equations with a turbulence model based on the eddyviscosity concept. The iced surface is treated as an arbitrary rough surface and its characteristics are related to an equivalent sand grain roughness, k_{s}. The roughness correlation between k_{s} and iced surfaces is established from experimental data. In this way, k_{s} is expressed as a function of freestream speed V¥
, static freestream temperature T_{s}, and liquid water content, LWC, of the air.
Our physical model for the ice accretion process, which is in the development phase, takes into account the presence of supercooled droplets entrained in the flow. These droplets follow trajectories that cause them to either be carried past a body or impinge upon the body. Upon impact with a clean surface, the droplets coalesce into larger surface droplets under the effects of surface tension and flow along the surface as dictated by the airflow. In our approach, we assume that these droplets then either freeze on the surface or are shed off from the surface because of the aerodynamic forces acting on the droplets. The ice accretion resulting from this initial freezing forms a rough surface. The heat transfer coefficient is then calculated for both laminar and turbulent flows again by using boundarylayer theory. The turbulence model is modified in order to account for roughness effects. Depending on our progress, ice shapes will be calculated with this approach and comparisons will be made with the prediction of LEWICE and experimental data for both glaze and rime ice shapes.
REFERENCE
 Ruff, G.A. and Berkowitz, B.M.: “User’s Manual for the NASA Lewice Ice Accretion Code (LEWICE),” NASA CR185129, May 1990.
NONLINEAR INTERACTIONS AND TEMPERATURE OSCILLATIONS IN LOWPRANDTL
MELT IN CZOHRALSKI MODEL: VALIDATION OF COMPUTATIONAL SOLUTIONS FOR GRAVITY DRIVEN AND ROTATORY FLOWS
V.I. Polezhaev, M.K Ezmakov, S.A. Nikitin
The Institute for Problems in Mechanics, Russian Academy of Sciences
Prospect Vernadskogo 101,117526 Moscow, Russia
Nonlinear interactions and temperature oscillations of gravitydriven fluid flow with rotation as typical
example of nonlinear coupling fluid system is studied as a benchmark of computer model of Czohralski
growth. A focus of the paper is concentrated on the new results of multiparametric analysis of nonlinear
coupling fluid flows in Czohralski growth (hydrodynamic model) using the axisymmetric model on the
basis of unsteady NavierStokes equations (Boussinesq approach).
Two types of numerical approaches: stream functionvorticity and velocitypressure formulations and two
types of numerical methods were used: "upwind" for qualitative pictures multiparametric simulation and
second order approximation for quantitative analysis of flow and temperature field structure, convective
instability and oscillations.
Parametric function of the characteristics in the melt can be written as follows:
A = (r, z, t, Rex, Rec, Gr, Pr, H/Rc, Rx/Rc, g, go)
here g indicates the boundary conditions and go  the initial values.
First part of the paper is devoted to the benchmark proposal [1] (H/Rx=1, Rx/Rc= 0.4, Pr=0.05) in the
case of prescribed temperatures Tx and Tc on the seed and crucible wall and isolation of the crucible
bottom and linear distribution on the melt surface. Results of tne calculations errors for different meshes
and for different values of Rex, Rec and Gr are presented. Fig.l shows selected pictures for general
gravity driven and rotation regimes. Critical Grashof number for onset of temperature oscillations in the
melt is also defined.
Second part of the paper related to the multiparametrical analysis of the general industrial groundbased
regime with counter rotation, without rotation as well as zero gravity regime with rotation .Dependencies
of oscillation amplitude on the Grashof number for the case of zero rotation as well as for the coupling
nonlinear case are presented. It is shown for above mentioned groundbased case that temperature
oscillations due to gravitydriven convection in Czohralski growth melt flow may be of two types:
RayleighBenard type due to cooling above of ctystal growth surface and side heating.
Experimental validation of computational solutions is discussed on the basis of comparisons with
experimental data for temperature oscillations in the GaAs melt [3] and in the water [4]. Nonsymmetry
effects and their impact on the above mentioned results are also discussed.
Computer simulation for second part of the paper is done on the basis of special version of PCbased
system "COMGA" as a part of computer laboratory [5]. Video visualisation pictures of unsteady
canvective interaction behavior of temperature and flow fields on the basis of this system are also
presented. Fig. 2 shows pictures of isotherms for general gravitydriven regime with counter rotation.
This work was supported partly by joint NASARSA program (project No TM11)
References
 Wheeler A.A. Four Test Problems for the Numerical Simulation of Flow in Czochralski Crystal
Growth. J. Crystal Growth,1990, v.102, p.691.
 Buckle U., Schafer M. Benchmark results for the numerical simulation of flow in Czochralski crystal
growth. / J. Crystal Growth 126,1993, p.682694.
 Kosushkin V.G., Polezhaev V.I., Zakharov B.G. Groundbased experiments and alternatives in GaAs
cryatal growth. In.: Proceedings of the Microgravity Science and Applications Session. International
Aerospace Congress, Moscow, ( Eds.) R.K. Crouch, V.I. Polezhaev,1995, p.141146.
 Jones A.D.W. T'he temperature field of a model Czochralski melt. J. Crystal Growth 69,1984, p.165
172.
 Polezhaev V.I., Etmakov M.K, Griaznov V.IL. et al. Computer Laboratory on Convective Processes in
Microgravity: Concepts, Current Results and Perspective. 46 Internatìonal Astronautical Congress, IAF95
J.3.11, October 26,1995, / Oslo, Notyvay.
AN ANALYTICAL MODEL FOR SOME CHANGE PHASE PROBLEMS IN BOUNDED DOMAINS COMPARED WITH NUMERICAL AND EXPERIMENTAL RESULTS
M.BENZADI  R.DUVAL  R.BENNACER  H.BEJI Laboratoire Mécanique et Sciences des Matériaux Rue d'Eragny 95031 Neuville Sur Oise France Phone : (33) Ol 34 25 68 70 Fax : (33) Ol 34 25 68 41 Email : boulmeza@mesopb.obspm.fr
Abstract  We present here an analytical model for some Stefan problems in a bounded domain.
This model, based on a connection between the unbounded domain approximation and the
quasistationnary approximation, predicts a priori what phase change front evolving look like for
small and large times,
Introduction
Problems which involve a phase change have received enormous attention in recent years,
especially in geophysics, geology, ice manifacture, enginering, mathematics and many other
technical fields. The main feature of such problems is the moving interface at which, the phase
change occurs, the way in which this interface moves has to be determined. Unfortunately, exact
solutions are available only for a few number of simplified problems those exact solutions for
small times. For greater times, no analytical solutions are available. therefore, numerical methods
appear to be the only practical methods to handle phase change problems providing that one can
successfully trace the moving interface. For that we want presenting an analytical model for some
Stefan problems in a bounded domain.
Solidification and melting problem in a rectangular cavity
We consider a rectangular cavity [O,L] x [O,H] which contain a fluid initially at a constant
temperature T_{2} > T_{f} (see fig.1 ). At time t = 0, the left face is reduced to a constant temperature
T_{1} < T_{f}. For t > 0, the two faces x = 0 and x = L are maintained at those temperatures respectively.
The upper and the lower faces are insulated. Hence, after some time, a layer of solid builds up
adjacent to the left face and moves to the right one.
The equations
Boundary conditions
Interface conditions
Initial conditions
T_{1}(t=O,x)=T_{2} and x(t = 0) = 0. (6)
For very early times , the solution is close to the analytical solution of the one
dimensional solidification problem in semi infinite region (see[4], p.283) of course
this solution is not valid when t becomes sufficiently large since when .
The asymptotic value of X(t) when is obtained very easily by solving equations (1) to (5)
in the study case:
by considering times this solution will connect the starting solution given above to the
asymptotic one.
The front evolving is given by:
Numerical simulation and validation
Fig.2 The position of the front as a function of time obtained by: (1) Numerical simulatïon (continue), (2) The analytic solution using the semiinfinite domain approximation (dashed), (3) The analytic solution obtained by the mode) exposed above.
The above analysis may easily be extended to many others change phase problems in a bounded
domain. As an illustration we consider this extension to two problems where we expose here just
one and the second it concern the problem of freezing with cylindrical symetry.
Multiphase problem in a cavity
We consider here that the material has two transitions temperatures T_{f1} and T_{F2}, at which latent
heats L_{H1} and L_{H2} are liberated. With T_{1} < T_{f1} < T_{f2} < T_{2}.Suffixes 1,2 and 3 respectively are used
here to describe the phases (T_{1}, T_{f1}) , (T_{f1} , T_{f2}) and (T_{f2} , T_{2}).
Fig.3 Twophase change problem in rectangular cavity
Using the same analysis as in the onephase change problem, the expressions of the position of
the two fronts, denoted respectively by X1(t) and X2(t) are:
a_{1}, a_{2}, a_{3} and a_{4} are calculated before, l_{1} and l_{2} are solutions to the equations (see[4], p.290)
the times of connection t_{1} and t_{2} are solutions to two equations of t_{1} and t_{2}
Experimental method
We compare also the results obtained by numerical and analytical methods to
experimental results wich are obtained from an experimental study based on an acoustic method
which has the advantage to be not destructive.
Conclusion
This study can be extended to solid/liquid phase change in porous media. In the
analytical part of the study, a model is developed which treats the entire domain as a single region
governed by one set of conservation equations. We have introduced equivalent conductivities in the
solid and liquid regions
and the latent heat become in the porous matrix eL_{H}.
Nomenclature
c specific heat, K thermal conductivity, k thermal diffusivity, L_{H} Latent heat, e porosity
References
 C. Beckermann and R. Viskanta, Natural convection solid/liquid phase change in porous media, Int. J. Heat Transfer, Vol. 31 (1988), N0.1, PP. 3546.
 A. M. Cames  Pintaux and M. Nguyen  Lamba, Finiteelement enthalpy method for discrete phase change, Numer. Heat Transfer, 9 ( 1986), pp.403417.
 Y. Cao and A. Faghri, A numerical analysis of Stefan problems for generalised multidimensional phasechange structures using the enthalpy transforming model, Int.J Heat Transfer, Vo1.32 (1989), N0.7, PP. 12891298
 H. S. Carslaw and J. C. Jaeger, Conduction of Heats in solids(2nd edn.), Oxford University Press, Oxford ( 1959) ,pp282296.
 J.M. Konrad, Anayse numérique des chaussées soumises à I'action du gel, du dégel et du trafic routier, Géotechnique et Informatique, Actes du colloque organisé par I'Ecole Nationale des ponts et chaussées, Paris, 2930 septembre, I er octobre 1992.
 R.W. Lyczkowski and Y.T. Chao, Comparaison of Stefan model with twophase model of cool drying, Int. J. Heat Transfer, Vol.27 ( 1984), N0.8,pp.l 1571969.
 M. Picasso, An adaptative finite element algorithm for a two dimensional stationary Stefan  like problem, Comput. Meth. appl Mech. Engrg. 124 ( 1995), pp. 213230.
