PHASE CHANGE I
Chairman: T. Kowalewski
TRANSPORT PROCESSES IN POLYMER EXTRUSION AND OPTICAL FIBER DRAWING
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 non-Newtonian 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.
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 non-Newtonian 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 time-dependent or time-independent, 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, time-independent viscoinelastic fluids without a yield stress are often
represented by the power-law model, given byl
where Kc 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 asl, 2
for a two-dimensional flow, with the velocity components u and w varying only with coordinate y.
Similarly, expressions for other two- and three-dimensional 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 model2. 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 well-layered,
with little bulk mixing, due to the high viscosity of these fluids. Viscous dissipation causes
the temperature to rise beyond the imposed barrel temperature2.
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 under-relaxation to achieve
convergence. Even a change of a few degrees in temperature in the vicinity of the softening point,
which is around 1600 oC for fused silica, can cause substantial changes in viscosity and thus in the
flow field and the neck-down profile in optical fiber drawing3, 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 steps4. Figure 3 shows such an iterative convergence for the neck-down 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 2-5 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 neck-down region is
smooth and well-layered because of the high viscosity. A typical temperature difference of 50-100
oC 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 neck-down, in the small diameter region, and plays an important role in maintaining the
temperatures above the softening point.
- Tadmor, Z. and Gogos, C., Principles of Polymer Processing, Wiley, New York,1979.
- Jaluria, Y., Heat and Mass Transfer in the Extrusion of Non-Newtonian Materials, Advances in Heat Transfer, Vol. 28, pp.145-230,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. 350-358, 1996. Also, Int. J. Heat Mass Transfer, Vol. 40, pp. 843-856,1996.
- Roy Choudhury, S., Jaluria, Y. and Lee, S.H.-K., Generation of Neck-Down Profile for Furnace Drawing of Optical Fiber, ASME Heat Transfer Div., Vol. 306, pp. 23-32,1995.
EFFECT OF BUOYANCY-DRIVEN 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 36849-5341, USA
+Technical University of Brno, College of Mechanical Engineering, Technická 2, 616 69 Brno, CZECH REPUBLIC
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, latent-heat-of-fusion 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 diffusion-controlled 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 diffusion-controlled 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 latent-heat-of-fusion thermal storage system. In high-power space systems, large amount of waste heat will be generated. As the heat loads are usually time-dependent, 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 latent-heat-of-fusion storage system, a thorough study of melting within spheres is needed.
The present work will utilize the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) procedure of Patankar11 to solve the unsteady coupled fluid flow and heat transfer equations governing the melting process. The fixed-grid 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 well-suited 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 phase-change 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 buoyancy-driven convection during melting within spherical containers has been undertaken. Emphasis is placed on low-Prandtl 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 solid-liquid interface will move into the specimen. As the extent of the liquid state is increased, buoyancy-driven convection will be promoted and a skewed solid-liquid 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 diffusion-controlled melting. This is due to the active role of buoyancy-driven 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 diffusion-controlled melting.
- Tao, L. C., Generalized Numerical Solutions of Freezing a Saturated Liquid in Cylinders and Spheres, AIChE Journal, Vol. 13, No. 1, pp. 165-169, 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. 1507-1516, 1974.
- Pedroso, R. I. and Domoto, G. A., Inward Spherical Solidification--Solution by the Method of Strained Coordinates, Int. J. Heat Mass Transfer, Vol. 16, pp. 1037-1043, 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. 2611-2616, 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. 397-406, 1996.
- Moore, F. E. and Bayazitoglu, Y., Melting within a Spherical Enclosure, Journal of Heat Transfer, Vol. 104, pp. 19-23, 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. 4174-4181, 1991.
- Saitoh, T., Kato, H. and Maruhara, K., Numerical Computation on Combined Natural Convection and Close-Contact Melting in a Cylindrical Energy Storage Capsule, Proc. Int. Conf. Comput. Sci., pp. 1-5, 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. 273-281, 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. 1-8, 1993.
- Patankar, S. V., Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, NY, 1980.
The last two authors acknowledge the support provided by the US-Czech 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 so-called 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 two-dimensional ice accretion computer program which was initially developed at the University of Dayton Research Institute and later modified by Ruff and Berkowitz1. 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 ad-hoc formulas based on Reynold’s analogy are replaced. A finite-difference method is used to solve the conservation equations with a turbulence model based on the eddy-viscosity concept. The iced surface is treated as an arbitrary rough surface and its characteristics are related to an equivalent sand grain roughness, ks. The roughness correlation between ks and iced surfaces is established from experimental data. In this way, ks is expressed as a function of freestream speed V¥
, static freestream temperature Ts, 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 boundary-layer 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.
- Ruff, G.A. and Berkowitz, B.M.: “User’s Manual for the NASA Lewice Ice Accretion Code (LEWICE),” NASA CR-185129, May 1990.
NON-LINEAR INTERACTIONS AND TEMPERATURE OSCILLATIONS IN LOW-PRANDTL
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
Non-linear interactions and temperature oscillations of gravity-driven fluid flow with rotation as typical
example of non-linear 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 non-linear
coupling fluid flows in Czohralski growth (hydrodynamic model) using the axisymmetric model on the
basis of unsteady Navier-Stokes equations (Boussinesq approach).
Two types of numerical approaches: stream function-vorticity and velocity-pressure 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  (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 ground-based
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 ground-based case that temperature
oscillations due to gravity-driven convection in Czohralski growth melt flow may be of two types:
Rayleigh-Benard 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  and in the water . Non-symmetry
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 PC-based
system "COMGA" as a part of computer laboratory . 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 gravity-driven regime with counter rotation.
This work was supported partly by joint NASA-RSA program (project No TM-11)
- 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.682-694.
- Kosushkin V.G., Polezhaev V.I., Zakharov B.G. Ground-based 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.141-146.
- Jones A.D.W. T'he temperature field of a model Czochralski melt. J. Crystal Growth 69,1984, p.165-
- 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, IAF-95
J.3.11, October 2-6,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
E-mail : firstname.lastname@example.org
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,
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 T2 > Tf (see fig.1 ). At time t = 0, the left face is reduced to a constant temperature
T1 < Tf. 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.
T1(t=O,x)=T2 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, 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
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 semi-infinite 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 Tf1 and TF2, at which latent
heats LH1 and LH2 are liberated. With T1 < Tf1 < Tf2 < T2.Suffixes 1,2 and 3 respectively are used
here to describe the phases (T1, Tf1) , (Tf1 , Tf2) and (Tf2 , T2).
Fig.3 Two-phase change problem in rectangular cavity
Using the same analysis as in the one-phase change problem, the expressions of the position of
the two fronts, denoted respectively by X1(t) and X2(t) are:
a1, a2, a3 and a4 are calculated before, l1 and l2 are solutions to the equations (see, p.290)
the times of connection t1 and t2 are solutions to two equations of t1 and t2
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.
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 eLH.
c specific heat, K thermal conductivity, k thermal diffusivity, LH Latent heat, e porosity
- C. Beckermann and R. Viskanta, Natural convection solid/liquid phase change in porous media, Int. J. Heat Transfer, Vol. 31 (1988), N0.1, PP. 35-46.
- A. M. Cames - Pintaux and M. Nguyen - Lamba, Finite-element enthalpy method for discrete phase change, Numer. Heat Transfer, 9 ( 1986), pp.403-417.
- Y. Cao and A. Faghri, A numerical analysis of Stefan problems for generalised multidimensional phase-change structures using the enthalpy transforming model, Int.J Heat Transfer, Vo1.32 (1989), N0.7, PP. 1289-1298
- H. S. Carslaw and J. C. Jaeger, Conduction of Heats in solids(2nd edn.), Oxford University Press, Oxford ( 1959) ,pp282-296.
- 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, 29-30 septembre, I er octobre 1992.
- R.W. Lyczkowski and Y.-T. Chao, Comparaison of Stefan model with two-phase model of cool drying, Int. J. Heat Transfer, Vol.27 ( 1984), N0.8,pp.l 157-1969.
- M. Picasso, An adaptative finite element algorithm for a two dimensional stationary Stefan - like problem, Comput. Meth. appl Mech. Engrg. 124 ( 1995), pp. 213-230.