Old Meetings: CHT97

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 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.

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 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 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 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 model². 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 temperature².

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 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⁴. 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.

REFERNCES

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
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, 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.
METHODOLOGY

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.
REFERENCES
ACKNOWLEDGEMENT

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 Berkowitz¹. 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, 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 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.

REFERENCE
1. 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 [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 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 [3] and in the water [4]. 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 [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 gravity-driven regime with counter rotation.

This work was supported partly by joint NASA-RSA program (project No TM-11)

References
1. Wheeler A.A. Four Test Problems for the Numerical Simulation of Flow in Czochralski Crystal Growth. J. Crystal Growth,1990, v.102, p.691.
2. Buckle U., Schafer M. Benchmark results for the numerical simulation of flow in Czochralski crystal growth. / J. Crystal Growth 126,1993, p.682-694.
3. 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.
4. Jones A.D.W. T'he temperature field of a model Czochralski melt. J. Crystal Growth 69,1984, p.165- 172.
5. 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 : 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₂ > T_f (see fig.1 ). At time t = 0, the left face is reduced to a constant temperature T₁ < 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₁(t=O,x)=T₂ 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 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 T_f1 and T_F2, at which latent heats L_H1 and L_H2 are liberated. With T₁ < T_f1 < T_f2 < T₂.Suffixes 1,2 and 3 respectively are used here to describe the phases (T₁, T_f1) , (T_f1 , T_f2) and (T_f2 , T₂).

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:

a₁, a₂, a₃ and a₄ are calculated before, l₁ and l₂ are solutions to the equations (see[4], p.290) the times of connection t₁ and t₂ are solutions to two equations of t₁ and t₂
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
1. 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.
2. A. M. Cames - Pintaux and M. Nguyen - Lamba, Finite-element enthalpy method for discrete phase change, Numer. Heat Transfer, 9 ( 1986), pp.403-417.
3. 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
4. H. S. Carslaw and J. C. Jaeger, Conduction of Heats in solids(2nd edn.), Oxford University Press, Oxford ( 1959) ,pp282-296.
5. 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.
6. 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.
7. 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.

Future Meetings | Past Meetings | Proceedings on Sale | Related Links

ICHMT World Wide Web Administrator