Modeling of Laser Keyhole Welding: Part I. Mathematical Modeling, Numerical Methodology, Role of Recoil Pressure, Multiple Reflections, and Free Surface Evolution HYUNGSON KI, PRAVANSU S. MOHANTY, and JYOTI MAZUMDER A three-dimensional laser-keyhole welding model is developed, featuring the self-consistent evolution of the liquid/vapor (L/V) interface together with full simulation of fluid flow and heat transfer. Important interfacial phenomena, such as free surface evolution, evaporation, kinetic Knudsen layer, homogeneous boiling, and multiple reflections, are considered and applied to the model. The level set approach is adopted to incorporate the L/V interface boundary conditions in the Navier–Stokes equation and energy equation. Both thermocapillary force and recoil pressure, which are the major driving forces for the melt flow, are incorporated in the formulation. For melting and solidification processes at the solid/liquid (S/L) interface, the mixture continuum model has been employed. The article consists of two parts. This article (Part I) presents the model formulation and discusses the effects of evaporation, free surface evolution, and multiple reflections on a steady molten pool to demonstrate the relevance of these interfacial phenomena. The results of the full keyhole simulation and the experimental verification will be provided in the companion article (Part II).
I. INTRODUCTION
LASERS, since their inception, have become critical in manufacturing processes involving joining, precision machining, and surface modification. Conventional manufacturing technologies, such as welding, cutting, and drilling, are continuously being replaced by new laser technologies. However, without a thorough understanding of the associated physics, the potential of lasers cannot be fully realized. Laser-aided materials processing is based on light-matter interaction and may result in four phases, such as solid, liquid, vapor, and plasma, over a wide range of temperatures. Hence, the involved physics has a wide spectrum ranging from heat conduction, melting, and fluid flow to vaporization, plasma formation, and laser-plasma interaction. Mathematical analysis of the process is computationally demanding and requires multidisciplinary approaches. One of the physics, the “keyholing” phenomenon, which occurs at high laser intensities, is critical to deep penetration welding. Upon irradiation with a high-intensity laser beam, the target material melts first, and subsequently, vaporization occurs. The vapor flux generates a recoil pressure on the evaporating surface. Also, there exists a huge temperature gradient on the liquid/vapor (L/V) interface due to the spatial distribution of laser-beam energy, which generates a thermocapillary force. Recoil pressure and thermocapillary force together provide the driving force for liquid ejection forming a vapor-filled cavity called the “keyhole.” As the keyhole deepens, the energy deposition pattern changes considerably due to the multiple reflections of the laser beam inside the hole. The major advantage of keyholing lies in its ability to HYUNGSON KI, Research Fellow, and JYOTI MAZUMDER, Professor, are with the Center for Laser Aided Intelligent Manufacturing, Mechanical Engineering Department, University of Michigan, Ann Arbor, MI 481092125. Contact e-mail:
[email protected] PRAVANSU S. MOHANTY, Assistant Professor, is with the Mechanical Engineering Department, University of Michigan–Dearborn, Dearborn, MI 48126-1409. Manuscript submitted November 8, 2001. METALLURGICAL AND MATERIALS TRANSACTIONS A
improve energy utilization and propagation tremendously. The existence of a keyhole, however, makes the process physics very complicated. The dramatic influence of the keyhole on the molten pool characteristics is illustrated in Figure 1. Figure 1(a) presents a high-speed charge coupled device (CCD) camera picture of the molten pool in the absence of a keyhole, i.e., in low-laser-intensity conduction-mode welding. Here, the melt pool is small with a shape close to a circle. The flow field, which is solely generated by the thermocapillary force, is very stable, and the free surface deformation is very small. In contrast, the presence of the keyhole dramatically changes the molten pool characteristics (Figure 1(b)). In keyhole mode (Figure 2), the melt pool is much longer and deeper, and the flow pattern is very complicated. Here, the difference and complexity of the physics come from the existence of the keyhole, i.e., when the laser power density is high enough to generate vaporization processes. Many modeling studies exist on conduction-mode welding.[11,12,37,62] In recent years, the focus has been to develop better computational models for the keyhole-mode welding process. Dowden et al.[17] obtained approximate solutions of laser-keyhole welding for low welding speeds. Kroos et al.[32] investigated the stability of a keyhole by prescribing the geometry of the keyhole to be cylindrical and concentric to the laser beam right through the workpiece. They only considered the energy and pressure balance in their model and neglected fluid flow. Kaplan[24] proposed a model by calculating the keyhole profile using a point-by-point determination of the energy balance at the keyhole wall. The model, however, neglected fluid flow. Dowden and Kapadia[18] suggested a mathematical model based on the instability resulting from the variable absorption at the workpiece due to reflections at the keyhole wall. Lankalapalli et al.[33] estimated penetration depth based on two-dimensional heat conduction on an assumed conical keyhole. Matsunawa and Semak[36] simulated the front keyhole-wall VOLUME 33A, JUNE 2002—1817
(a)
(b) Fig. 1—(a) Conduction-mode laser welding (D.D. Voelkel: Ph.D. Thesis, University of Illinois at Urbana-Champaign, IL, 1992). (b) Keyhole-mode laser welding.
behavior assuming that only the front part of the keyhole wall is exposed to the high-intensity laser beam. Semak et al.[46] also developed a model for the transient behavior of the front keyhole, assuming that the keyhole wall is dominated by evaporation-recoil-pressure-driven melt expulsion from the beam interaction zone. They showed that the front part of the keyhole has a steplike shape, which occurs when the component of the keyhole velocity along the sample surface exceeds the welding speed. Solana and Ocana[52] developed a model to determine the three-dimensional weld pool and keyhole geometry by considering heat conduction, evaporation losses, and evaporation effects at the keyhole surfaces, neglecting the fluid flow. Sudnik et al.[54,55] found a linear correlation between the depth and the length of the weld pool with the laser power intensity. In their model, they used a semiempirical approach, i.e., assuming a simplified pattern of flow trajectories between plane sections and determining the flow velocities analytically on the basis of continuity considerations. Fabbro and Chouf [19] modeled the keyhole, based on a drilling or penetration velocity whose combination with the welding speed causes the inclination of the front keyhole wall. They analyzed the problem by considering front and rear keyhole walls separately. From the preceding discussion, it is evident that the existing modeling studies are oversimplified and do not present the complete picture of the associated physical phenomena. The interface shape and the interaction physics are intimately coupled. In other words, the surface shape affects the interaction physics, which, in turn, affects surface shape. Modeling these phenomena necessitates an interdisciplinary approach involving heat transfer, fluid mechanics, and phase transformation while tracking the L/V interface self-consistently. As far as the authors know, none of the existing models fully solved the thermal and the flow fields together with the evaporation process and the transient evolution of the keyhole. This article presents a three-dimensional laserkeyhole welding model with the focus on self-consistent simulation of keyhole evolution coupled with complete fluid flow and heat transfer. As a heat source, this study assumes CO2 lasers even though the developed model is capable of simulating other types of lasers, such as high-power yttrium aluminum garnet (YAG) lasers. II. KEYHOLE EVOLUTION A. Free Surface Evolution: Narrow-Band Level Set Method
Fig. 2—Schematic of laser keyhole welding process. 1818—VOLUME 33A, JUNE 2002
It is no exaggeration to say that keyhole welding is a L/V interface process. In order to account for the interaction of interface shape and the process physics, a method of tracking the free surface is needed. The level set method, which was recently developed by Osher and Sethian (1988[40]) has been very successful in tracking complex, free surface movements. Indeed, this method is today considered one of the best methods to capture moving interfaces.[1–3,39,47,50] Since its inception, the method has been used to solve many problems involving moving boundaries, such as fluid mechanics,[39,56] solidification,[13,48] flame propagation,[50] and stereolithography.[1,2,3] The method has many advantages, since it transforms an interface tracking problem into a partial differential equation. First, the implementation is relatively easy, since the solution METALLURGICAL AND MATERIALS TRANSACTIONS A
(a)
the x ⫺ y plane, x2 ⫹ y2 ⫺ 1 ⫽ 0, is a zero level set of a function, (x, y) ⫽ x2 ⫹ y2 ⫺ 1 ⫽ C when the constant, C, is zero. Thus, instead of directly dealing with the circle, we can define a function, (x, y), of which the zero level set is the surface of interest. Since (x, y) is defined in the entire x ⫺ y plane, it is much easier to solve it numerically. The level set equation is a novel way of viewing the welldefined kinematic boundary condition of the interface. The equation, which has been thought to be valid only for the interface, can be defined in the entire space that contains the interface by considering it as the one level set of a wellcreated function, . To extend the zero level set (the physical interface) to the entire domain, the level set function should be well-defined, and the signed distance function is chosen as the standard: (x, t) ⫽ ⫾d [1] where d is the actual distance from the zero level set, 0; and the plus (minus) sign denotes the outside (inside) of the 0. The general form of the level set equation is ⭸ [2] ⫹ F0앚앚 ⫹ U(x, y, z, t) ⭈ ⫽ ⑀앚앚 ⭸t where ⑀ is a number, and is curvature. Here, F0 is the speed function, U(x, y, z, t) is convection velocity, and the term on the right-hand side accounts for the surface movement due to interface curvature. Equation [2] must be modified to be used in this problem. The surface curvature effect will be taken into account when formulating the Navier– Stokes equation with the level set function in Section III–A. If the surface recession speed due to the evaporating mass flux, is Fe , then, Eq. [2] becomes ⭸ [3] ⫹ Flv앚앚 ⫽ 0 ⭸t where Flv ⫽ Fe ⫹ n ⭈ U(x, y, z, t) [4] and [5] 앚앚 An advancement to this technique, the “narrow-band level set method,”[47] can reduce the computational cost significantly by taking the computational domain as a tiny band around the interface we are interested in (Figure 3(b)). In this study, the narrow-band level set method has been used to simulate the free surface evolution. In the modeling of laser-keyhole welding, surface normal and curvature values are used to calculate capillary and thermocapillary terms, which are dependent on the interface geometry. Especially when it comes to the modeling of multiple internal reflections using the ray tracing method, it is imperative to compute the surface normal accurately and easily, since the accuracy of the energy deposition profile determines the accuracy and reliability of the whole solution, such as fluid flow and heat transfer. The surface normal (Eq. [5]) can be simply discretized using the central difference scheme. Since the surface geometry may undergo a jump at a corner, Sethian et al.[48] have suggested a technique to approximate the normal in which the one-sided difference approximations to the unit normal in each possible direction are all obtained and then averaged. Their method is adopted in this study. n⫽
(b) Fig. 3—(a) An example of the level set function. (b) Narrow band method.
methods for partial differential equations are generally well known. Second, geometric complexities, such as merging and splitting of multiple surfaces, can be handled. Those geometric complexities are very troublesome for other front tracking methods, such as the volume of fluid method, the particle or node tracking method, the cell and marker method, etc. Finally, geometric properties, such as normal and curvature, can be easily calculated from the level set values. This is very important in the modeling of laser welding, where normal and curvature have physical importance. In implementing this method, the interface of interest is embedded as the zero level set of one-order-higher dimensional space (Figure 3). This zero level set is the surface where all the physics should be implemented to evaluate the force that moves this surface. For example, a circle in METALLURGICAL AND MATERIALS TRANSACTIONS A
VOLUME 33A, JUNE 2002—1819
p s pv ⫽ pa pa
冒
pv ps
[10]
sTs ps ⫽ pv vTv
[11]
av ⫽ 冪␥vRvTv / 冪␥aRaTa aa
[12]
where a is the sound speed. For the supersonic case, the Mach number outside of the Knudsen layer is always 1, and the rarefaction fan accelerates the flow to supersonic. The corresponding relations are ps pok pv ⫽ pa pv pa
Fig. 4—Flow structure outside the evaporating surface.
冢
uv ⫽ 冪␥v /2 Mv, 冪2RvTv
[6]
where u, R, ␥, and M denote velocity, gas constant, specific heat ratio, and Mach number, respectively. The corresponding temperature, density, and pressure jump conditions are[31] 2
␥v ⫺ 1
m 冢冪1 ⫹ 冢␥␥ ⫹⫺ 11 m2 冣 ⫺ 冪 ␥ ⫹ 1 2 冣 2
v
v ⫽ s
冪T 冋冢m Ts
2
v
v
⫹
v
冣
1 m2 m e erfc (m) ⫺ 2 冪
册
[7]
[8]
1 Ts ⫹ [1 ⫺ 冪 m exp (m2) erfc (m)] 2 Tv
冢
av ␥a ⫹ 1 av pv Mv ⫽ 1 ⫹ ␥aMv pa aa 4 aa
冋
⫹ 1⫹
冢
␥a ⫹ 1 av Mv 4 aa
冣册
2 1/2
[9]
冣
respectively. If the Mach number outside the Knudsen layer reaches 1, a rarefaction fan is produced. Therefore, subsonic and supersonic cases should be dealt with separately. When the Mach number is less than 1, 1820—VOLUME 33A, JUNE 2002
[13]
2␥ v
When the L/V interface temperature reaches the boiling point, evaporation begins to occur. There exists a very thin layer of several mean free paths, called a kinetic Knudsen layer, just outside the L/V interface. Across this layer, the continuum hypothesis fails, and a steep change in temperature, pressure, and density occurs. Therefore, it is treated as a mathematical discontinuity at the interface. Jump conditions with back pressure were derived by Knight[31] and have been adopted in the authors’ previous work on laser drilling.[29] The method will be briefly described here. Figure 4 shows the flow structure generated by the evaporation process. The circled letters are subscripts used in the equations below. Defining
Tv ⫽ Ts
pok pok ⫽ 0.206 ps ps
pok ␥v ⫺ 1 2 ⫽ ⫹ M pv ␥v ⫹ 1 ␥v ⫹ 1 v
B. Evaporation at the L/V Interface
m⫽
冒
av ␥vRvTsTok ⫽ aa ␥aRaTaTs
冢
1/2
冣冢
冣
␥ v⫺1
[14]
␥v ⫺ 1 2 ⫹ M ␥v ⫹ 1 ␥v ⫹ 1 v
⫺1
冣
Tok ⫽ 0.669. Ts
[15]
[16]
The net mass loss due to evaporation can be calculated as[61] m ˙ evap ⫽ s
冪 2 ⫺ 冪 2  F RTs
RTv
v
⫺
(m) ⫽ lFe
[17]
where
⫽
2(2m2 ⫹ 1) 冪Tv /Ts ⫺ 2冪 m F ⫺ ⫹ 冪Tv /TsG⫺
F ⫺ ⫽ 冪 m (⫺1 ⫹ erf (m)) ⫹ exp (⫺m2) G⫺ ⫽ (2m2 ⫹ 1) (1 ⫺ erf (m)) ⫺
2 m exp (⫺m2) 冪
In the preceding expression, the factor, , is the modification factor, which accounts for the back-scattered flux,[61] and Fe is the speed function of the L/V interface due to evaporation. To complete the Knudsen-layer jump relations, the Clausius–Clapeyron relation is used in almost all research including Knight’s original work,[31] assuming that the L/V interface is always saturated. In reality, however, it is far from the actual heating pattern of material caused by the intense laser beam, and the interface is highly superheated. Therefore, the degree of superheat is approximated using the momentum exchange balance across the Knudsen layer. In the state of thermodynamic equilibrium between gas and liquid, the evaporating and back-scattered fluxes are equal and contribute equally to the surface pressure. In a vacuum, however, since there is no back-scattered flux, the interface pressure is around 50 pct of the saturation pressure. Rigorous calculation reveals that the pressure is around 0.56 psat(Ts)[5]. Thus, in a general situation, back-scattered flux is small compared to the evaporating flux, and the actual surface pressure can be expressed as METALLURGICAL AND MATERIALS TRANSACTIONS A
Fig. 6—Surface tension value for liquid iron.
Fig. 5—P–T relationship[26].
ps ⬵
svs ⫹ vvv psat(Ts) 2svs
[18]
J ⫽ Nl
C. Homogeneous Boiling The multiple reflections may increase the energy to such a level that it can cause homogeneous boiling in some locations and can accelerate the movement of the L/V interface. However, there exists a maximum temperature, the “critical point” temperature, that a liquid metal can attain. At the critical point, the distinction between liquid and gas vanishes, and only the fluid state exists.[44,35,26,60,53,34,21,22,5] It is reported that intense heating of the substrate can cause homogeneous boiling right beneath the L/V interface when the interface temperature reaches around 80 pct of the critical point due to a significant drop in surface tension and a huge amount of liquid superheat.[29] Craciun et al.[14] examined the role of subsurface boiling for micrometer-sized droplet emission during pulsed-laser ablation of Ge. Using scanning electron microscopy, they observed circular micron-deep cavities and suggested that those cavities were formed by the release of gas caused by the subsurface boiling. The phenomenon of homogeneous nucleation is summarized in References 10, 15, and 49. During boiling, the critical radius of a bubble, re , can be computed from the following equation: re ⫽
2 Psat(Tl) exp [vl(Pl ⫺ Psat(Tl))/RTl] ⫺ Pl
[19]
Here, Tl , Pl , Psat(Tl), vl , and R are liquid temperature, liquid pressure, saturation pressure at Tl , specific volume of liquid, and the gas constant, respectively. The pressure inside the bubble, Pv , and the specific volume of the bubble, vv , can be calculated from the Young–Laplace equation and modified Van der Waals equations, respectively: Pv ⫽ Pl ⫹ Pr ⫽
2 re
8Tr 3 ⫺ 3vr ⫺ 1 v2r
METALLURGICAL AND MATERIALS TRANSACTIONS A
Here, Pr ⫽ P/Pcr , Tr ⫽ T /Tcr , and vr ⫽ v/vcr . Since surface tension values are not known for the higher temperatures, they are extrapolated to the critical point using the experimental data up to the normal boiling temperature and the theoretical asymptotic behavior[16] near the critical point (Figure 6). The homogeneous evaporation rate, which is the number of bubbles with critical radius per unit time per unit volume of liquid, can be computed by
[20] [21]
冋
6 m(2 ⫺ Pl /Pve)
册
1/2
exp
冋
册
⫺4 re2 3kBTl
[22]
where Nl , m, Pve , kB are the number of liquid molecules per unit volume, the mass of one molecule, the pressure in the vapor bubble, and the Boltzmann constant (1.3807 ⫻ 10⫺23 J/K), respectively. The number of molecules in the bubble, ne , is now ne ⫽
Vd NA vv M
[23]
where Vd is the volume of the bubble, and NA is Avogadro’s number (6.022 ⫻ 1023). Vd is expressed as: Vd ⫽
4 r 3e 3
[24]
M is molecular mass (55.85 for iron). Finally, mass flux and energy flux due to homogeneous bubble generation near the critical point are m ˙ ho ⫽ J ne mm
[25]
冢2 k T ⫹ m L ⫹ 3 r 冣
q˙ho ⬵ J ne
3
4
B
m v
2
e
[26]
which complete the mass and energy balances. It is assumed that the vapor molecule inside the bubble is monoatomic, and mm is the atomic mass. Once a bubble of critical size is generated in the superheated liquid pool, it will grow or collapse depending on the surrounding conditions. And, boiling, whether the distance scale is atomically small or much larger, has prohibitive kinetic obstacles because it requires bubble diffusion if the bubbles are formed other than at the outer surface.[27] In the case of intense laser energy, bubbles generally grow due to the high energy influx. This homogeneous bubble growth, when combined with instability, will possibly lead to the so-called explosive boiling. Explosive-type boiling has been observed in a laser-matter interaction experiment.[20] Here, the general bubble-growth relation is used, ignoring the explosive boiling phenomena. VOLUME 33A, JUNE 2002—1821
D. Multiple Internal Reflections In keyhole welding, the energy is transferred into the substrate by multiple internal reflections inside that cavity. With this mechanism, generally, a metal surface with a very low absorptivity, such as 0.10, can behave similar to a blackbody, since a laser beam transfers most of its energy to the keyhole surface. Thus, multiple reflection phenomena determine how the energy is transferred from the laser beam to the workpiece, and, most importantly, all other physics occurring in the laser materials processing, such as fluid flow and heat transfer, are highly dependent on it. In short, to better understand laser-keyhole welding, an accurate model of multiple internal reflections is crucial. Several studies on multiple reflections have been done as part of laser-keyhole welding/drilling modeling.[19,25,38,51,59] Kar et al.,[25] in their laser drilling research, used the ray tracing technique for their two-dimensional axisymmetric drilling model. They assumed a parabolic hole shape for the simulation. Milewski and Sklar[38] simulated three-dimensional multiple internal reflections for V-groove weld joints, using the ray tracing technique and assuming the conical beam shape. Solana and Negro[51] investigated the effect of multiple reflections on keyhole welding, also using the ray tracing model and assuming an initial conical keyhole shape. Wei and Ho[59] used a paraboloid of revolution as an assumed profile of the drill hole for their laser drilling research. Fabbro and Chouf [19] employed the ray tracing method in their laser-keyhole welding study. In this study, a multiple reflection model has been developed based on the level set method. The ray tracing technique has been adopted, considering the specular mode of reflection. As emphasized earlier, use of the level set information to calculate the surface normals is the key idea of this model. This model has an ability to deal with arbitrary hole shapes if accurate level set information is provided. The normal vector can be calculated at each point using the level set values. In finding the next reflection point, the bisection method[9] is employed. At each time of reflection, the laser energy is absorbed into the wall, according to the material’s absorptivity. The whole process can be simulated by simply repeating these steps until the ray escapes from the computational domain. The grids for multiple reflection calculations are much finer than the main grid system for the sake of a good resolution of solutions. III. FLUID FLOW AND HEAT TRANSFER A. Momentum Equation with L/V Interface Boundary Conditions Boundary conditions for the moving L/V interface are very difficult to implement using conventional methods. The level set-based method, however, has the ability to incorporate the interface boundary conditions in the formulation, since it provides information on the interface shape, no matter how complicated it is. Besides, in laser-material processing, surface tension has a spatial distribution due to the huge temperature gradient. Thus, for accurate force balance at the interface, surface tension must be considered to vary spatially. Ki[28] incorporated the thermocapillary boundary condition in the momentum equation by modifying the level set approach proposed by Sussman et al.[56] The details of the method can be found in Reference 28. 1822—VOLUME 33A, JUNE 2002
The momentum equation can be described as follows:
Du d ⫽ ⭈ T ⫺ n* ⫺ sT ␦ ( ) Dt dT
冢
冣
[27]
Here, T is the stress tensor, and n* is the surface normal directing liquid to gas. The two terms inside the bracket on the right-hand side are the L/V interface boundary conditions. The n* ␦ ( ) term denotes the capillary effect that acts in the normal direction due to the interface curvature and surface tension. This capillary force is the major restoring force, which counteracts the recoil pressure. Note that the capillary and thermocapillary terms are successfully combined in the governing equation, not as boundary conditions. The recoil pressure will be incorporated as part of the implementation of the Symmetrically Coupled Gauss Seidel (SCGS) method, in Section III.E, which will complete the L/V interface boundary conditions. B. Energy Equation with L/V Interface Boundary Conditions In this section, the boundary conditions for the energy equation are incorporated in the governing equation. A Gaussian laser-beam distribution is considered, and it is assumed that the beam profile and radius do not change in the z direction: q⬙laser (x, y) ⫽
冢
2Plaser 2(x2 ⫹ y2) 2 exp ⫺ Rb R2b
冣
[28]
where Plaser is the laser power, and Rb is the beam radius. The actual energy deposition pattern on the L/V interface is modified by many factors. First, the multiple reflections inside the laser-created cavity greatly distort the original laser-beam distribution. Second, extra heat is lost due to the latent heat required by the L/V phase transformation when the surface temperature is above the normal boiling point. Third, radiation loss may not be ignored since the surface temperature can attain very large values. Therefore, the actual laser energy on the surface is q⬙ ⫽ q⬙L/V ⫺ lLvFe ⫺ ⑀ (T 4 ⫺ T 4⬁),
[29]
where q⬙L/V is the spatial laser-beam distribution in threedimensional space after multiple reflections, lLvFe is the energy loss due to evaporation, and ⑀ (T 4 ⫺ T 4⬁) is the energy loss due to radiation. The and ⑀ denote the Stefan– Boltzmann constant ( ⫽ 5.67 ⫻ 10⫺8W/m2 ⭈ K4) and emissivity, respectively. As seen in Eq. [29], convection heat loss is not included, since convection in the gas phase, as well as the liquid phase, are solved together. The final form of the energy equation for both liquid and gas phases, including the L/V interface boundary conditions is given below. ⭸(CpT ) ⫹ u ⭈ (CpT ) ⫽ kⵜ2T ⫹ q⬙ ␦ ( ) ⭸t
[30]
Here, ␦ ( ) is the delta function of the level set values. C. Melting and Solidification: S/L Interface Boundary Conditions In this section, the phase transformation at the S/L interface will be considered. Due to the absorption or release of METALLURGICAL AND MATERIALS TRANSACTIONS A
latent energy, phase change problems are nonlinear. Moreover, the interface between solid and liquid can be morphologically very complex, forming a mushy zone that contains both solid and liquid states. Due to absence of a discrete phase change, the level set method would be inappropriate for the S/L interface problem. In this study, a method developed by Bennon and Incropera[6] is adopted, which is an extension of the classical mixture theory. They developed a consistent set of continuum equations for the conservation of mass, momentum, energy, and species in a binary, S/L phase change system, assuming laminar, constant property flow in the liquid phase and characterizing the mushy region as a porous solid of isotropic permeability, K. The details and capabilities of the model have been demonstrated through many works.[7,8,23,42,43] In this section, the method will be summarized, closely following Bennon et al.[6,7] and Prakash et al.[42] with minor modifications. Development of conservation equations is based on the following principles. First, mixture components may be viewed as isolated subsystems if interactions with other mixture components are properly treated. Second, all properties of the mixture are mathematical consequences of the component properties. Finally, the mean collective-mixture behavior is governed by equations similar to those governing the individual components. Defining mass fraction, f, and volume fraction, g, the velocity, density, thermal conductivity, mass diffusion coefficient, concentration, and enthalpy for the liquid and solid mixture are
K ⫽ K0
g3l (1 ⫺ gl)2
[38]
The mixture transport equations reduce to the standard momentum equation as K → ⬁ and limits to an appropriate equation for the solid as K → 0. For implementation purposes, it is better to express the energy equation in a temperature form, and this can be achieved by linearizing the phase enthalpies. The phase enthalpies are hs ⫽ hl ⫽
兰
Te
0
兰
T
Cps dT
0
Cps dT ⫹ Lm ⫹
[39]
兰
T
Te
Cpl dT
[40]
where Te is the eutectic temperature. Defining average specific heats Cpl , Cps , and a constant, ᏸ, as Cps ⫽ Cpl ⫽ ⫽ ⌬ CpTe ⫹ Lm
1 T
兰
T
0
1 T ⫺ Te
Cps dT
兰
T
Te
Cpl dT
ᏸ ⫽ (Cps ⫺ Cpl)Te ⫹ Lm
[41] [42] [43]
we obtain
u ⫽ fsus ⫹ fl ul
[31]
hs ⫽ CpsT
[44]
⫽ gss ⫹ gll
[32]
hl ⫽ CplT ⫹ ᏸ
[45]
[33]
h ⫽ CpT ⫹ flᏸ
[46]
k⫽
⫺1
冢k ⫹ k 冣 gs
gl
s
l
The final form of the energy equation is
D ⫽ f lD l
[34]
c ⫽ f s cs ⫹ f l cl
[35]
h ⫽ fshs ⫹ flhl
[36]
respectively. Using the defined mixture variables, continuity equation, momentum equations in x, y, and z directions, an energy equation and solute conservation equations are given elsewhere[6] and are not repeated here. The momentum, energy, and species conservation equations contain the phase interaction terms. The phase interaction part is crucial and needs to be handled carefully. For a wide range of S/L phase change problems, the multiphase region is characterized by a fine permeable solid matrix, which can either be stationary or constrained to free body translation. In this study, it is assumed that the solid phase translates with the prescribed scanning velocity, Vs. For such systems, analogies can be drawn to flows through porous media. The Darcian damping force is written as F⫽
l gu K l r
[37]
where K is the isotropic permeability, and ur ⫽ ul ⫺ us is the superficial liquid velocity relative to the velocity of the porous solid. In the present analysis, permeability is assumed to vary with liquid volume fraction according to the Kozeny– Carman equation:[4] METALLURGICAL AND MATERIALS TRANSACTIONS A
⭸( flᏸ) ⭸( CplT ) ⫹ u ⭈ (CplT ) ⫽ ⭈ (kT ) ⫺ ⭸T ⭸t ⭸( fs⌬CpT) ⫹ ⫹ ⭈ [ fs(ᏸ ⫹ ⌬CpT )]Vs ⭸t
[47]
At this point, a scheme to update the mass fraction is necessary. Prakash and Voller[42] presented a stable and rapidly convergent algorithm for updating the mass fraction. The liquid mass fraction is related to the temperature and composition as T ⫽ F( fl, c)
[48]
which denotes that the temperature, T, is some function, F, of the liquid fraction, fl , and concentration, c.
D. Generalized Transport Equations Combining the results from Sections III–A, B, and C, the generalized transport equations with S/L and L/V boundary conditions can be written as follows. Continuity equation: ⭸ ⫹ ⭈ ( u) ⫽ 0 ⭸t
[49]
VOLUME 33A, JUNE 2002—1823
Table I. Parameters for the Governing Conservation Equations ⌽
⌫
1
0
0
x-momentum
u ⬅ u⬘ ⫺ Vs
⫺
y-momentum
v
z-momentum
w
Energy
T
k/Cp
Equation Mass
S⌽(x, y, z, t)
l ⭸p ⫺ ex ⭈ ( n* ⫺ s )␦ ( ) (u ⫺ Vs) ⫺ k l ⭸x l ⭸p ⫺ ey ⭈ ( n* ⫺ s )␦ ( ) ⫺ (v) ⫺ k l ⭸y l ⭸p ⫺ ez ⭈ ( n* ⫺ s )␦ ( ) ⫺ (w) ⫺ k l ⭸z ⭸( flᏸ) ⭸( fs⌬CpT ) ⫹ ⫹ ⭈ [ fs(ᏸ ⫹ ⌬CpT )]Vs ⫹ (q⬙L/V ⫺ lLvFe ⫺ ⑀ (T 4 ⫺ T 4⬁))␦ ( ) ⫺ ⭸t ⭸t
Momentum equations: ⭸(u) ⫹ ⭈ ( uu) ⫽ ⭈ (l u) ⭸t l
[50] l ⭸p (u ⫺ Vs) ⫺ ⫺ ⫺ ex ⭈ ( n* ⫺ s )␦ ( ) K l ⭸x
l ⭸(v) ⫹ ⭈ ( uv) ⫽ ⭈ (l v) ⫺ (v) ⭸t l K l ⭸p ⫺ ey ⭈ ( n* ⫺ s )␦ ( ) ⫺ ⭸y
[51] E. Solution Scheme and Incorporation of Recoil Pressure
⭸(w) ⫹ ⭈ ( uw) ⫽ ⭈ (l w ) ⭸t l
[52] l ⭸p ⫺ ez ⭈ ( n* ⫺ s )␦ ( ) ⫺ (w) ⫺ K l ⭸z
Energy equation: ⭸( CplT ) ⭸( flᏸ) ⫹ u ⭈ ( CplT ) ⫽ ⭈ (kT ) ⫺ ⭸t ⭸t ⫹
⭸( fs⌬CpT ) ⫹ ⵜ ⭈ [ fs(ᏸ ⫹ ⌬CpT )]Vs [53] ⭸t
⫹ (q⬙L/V ⫺ lLvFe ⫺ ⑀ (T 4 ⫺ T 4⬁))␦ ( ) where ex, ey, and ez are unit vectors in the x, y, and z directions, respectively. Since the energy source scans over the substrate, the Eulerian frame is the most desirable frame to develop a multidimensional problem. To relocate the Cartesian coordinate system from the reference frame of the substrate under consideration to that of the reference frame of the scanning energy source, the following coordinate transformation is used. x ⫽ x⬘ ⫺ Vst
[54]
where Vs is the translation velocity of the workpiece. Equations [49] through [53] can be cast into the general form in a Cartesian coordinate system as follows, using the Einstein summation convention.
冢 冣
⭸ ⭸ ⭸ ⭸⌽ ( ⌽) ⫹ (ui⌽) ⫽ ⌫ ⭸t ⭸xi ⭸xi ⭸xi ⫹ S⌽(x1, x2, x3, t) 1824—VOLUME 33A, JUNE 2002
where is the density of material, t represents time, ⌽ is the dependent variable, ⌫ is diffusivity, and ui is the i-directional Cartesian velocity component, respectively. The S⌽(x, y, z, t) term encapsulates any additional terms that do not follow this generalized formulation. The governing equations describing conservation of mass, momentum, energy, and species can be obtained by appropriate assignment of ⌽, ⌫, and the source term, S⌽(x, y, z, t). A listing of the terms is given in Table I.
[55]
The formation procedure for transport of a general property, ⌽, is based on Patankar’s method.[41] Scalar quantities, such as pressure, temperature, the level set function, and material properties, are stored at the nodes of an ordinary control volume. To evaluate velocity components, staggered grids, which are centered around the cell faces, are used.[41,58] The momentum and energy equations are discretized as follows. aP ⌽ ⫽ aE⌽E ⫹ aW⌽W ⫹ aN⌽N ⫹ aS⌽S a⌽ P
[56]
aP ⫹ aT⌽T ⫹ aB⌽B ⫹ b ⫹ (1 ⫺ ␣⌽) ⌽* ␣⌽ P Here, ␣⌽ is the relaxation factor for the dependent variable, ⌽; and ⌽*P is taken as the value of ⌽P from the previous iteration. In this study, temperature and velocity variables are all underrelaxed to avoid unpredicted divergence, and a relaxation factor of 0.8 is chosen for u, v, w, and T. In this analysis, an upwind-differencing scheme is used to accommodate advection of the dependent variable, ⌽. Likewise, the continuity equation is discretized by integrating it over a cell for scalar variables, and we have an equation as follows: (ae)cue ⫺ (aw)cuw ⫹ (an)cun ⫺ (as)cus ⫹ (at)cut ⫺ (ab)cub ⫹ bc ⫽ 0
[57]
The details are found in Reference 28. In this study, the energy equation is solved by an implicit method. Since the energy equation is solved for all phases, including the solid phase, the equation can be solved by any method. For velocity and pressure coupling equations, the coupled point-solution scheme is adopted since solid phase, as well as liquid and vapor phases, exist in the computational domain. In addition, the scheme offers easy extension to METALLURGICAL AND MATERIALS TRANSACTIONS A
three dimensions and addresses the velocity-pressure coupling in a more direct manner than segregated methods. The SCGS point-relaxation scheme developed by Vanka[57] is employed in this work. The SCGS method can also be vectorized with ease. The implementation details can be found in References 45 and 57. The SCGS method gives a matrix equation: Au ⫽ r
[58]
where
A⫽
冢
(ap)uijk 0 0 0 0 0 c ⫺(aw)ijk
0 0 0 (aP)ui⫹1jk v 0 (aP)ijk 0 0 0 0 0 0 c c (aE)ijk ⫺(aB)ijk
0 0 0 v (aP)ij⫹1k 0 0 c ⫹(aT)ijk
冣
0 0 (A)uijk 0 0 ⫺(A)ui⫹1jk v 0 0 (A)ijk v 0 0 ⫺(A)ij⫹1k w w (aP)ijk 0 (A)ijk w w 0 (ap)ijk⫹1 ⫺(A)ijk⫹1 c c ⫺(aS)ijk (aN)ijk 0
u⫽
r⫽
From the preceding discussion, it is evident that the laser welding process is extremely complex, and the present formulation incorporates the transport behavior of three interacting phases, i.e., vapor, liquid, and solid, having widely varying physical properties in a single computational domain. The following simplifying assumption are used in this study.
By forward substitution, the following expression for ␦p is obtained. [59]
c c c Rijk ⫺ (aE)ijk Rui⫹1jk /(aP)ui⫹1jk ⫹ (aW)ijk Ruijk /(aP)uijk c w w c w w ⫺(aN)ijk Rijk⫹1 /(aP)ijk⫹1 ⫹ (aS)ijk Rijk /(aP)ijk c v v c v v ⫺(aT)ijk Rij⫹1k /(aP)ij⫹1k ⫹ (aB)ijk Rijk /(aP)ijk
冦 冦
c c (aE)ijk (A)ui⫹1jk /(aP)ui⫹1jk ⫹ (aW)ijk (A)uijk /(aP)uijk c w w c w w ⫹(aN)ijk (A)ijk⫹1 /(aP)ijk⫹1 ⫹ (aS)ijk (A)ijk /(aP)ijk c v v c v v ⫹(aT)ijk (A)ij⫹1k /(aP)ij⫹1k ⫹ (aB)ijk (A)ijk /(aP)ijk
Then, ␦u, ␦v, and ␦w can be solved by back-substituting ␦p. Until this point, all L/V interface boundary conditions have been accounted for except the recoil pressure. The recoil pressure can be incorporated into the flow field as part of the solution scheme for velocity pressure coupling equations. When updating the pressure in Section III–E, Eq. [59] should be replaced by recoil ␦pijk ⬅ pijk ⫺ p*ijk
[60]
if evaporation occurs at the node point, and this completes the L/V interface boundary conditions. METALLURGICAL AND MATERIALS TRANSACTIONS A
1. The number of grid points around the keyhole must be large enough to capture the transient movement of a keyhole. Since the length of a melt pool (7 to 8 mm) is, in general, very long compared to the keyhole diameter (500 to 600 microns), the total number of grid points in the direction of scanning speed must be very large even though a nonuniform grid system is adopted. 2. In this study, fluid flow, heat transfer, melting/solidification, microstructure evolution, and free surface movement are all solved in three dimensions, together with multiple reflections, evaporation, and homogeneous boiling phenomena. 3. The time-step of this problem must be very small (⬃10⫺6 seconds) to satisfy the convergence criteria requested by the governing equations. In addition, the time-step is also restricted by movement of the L/V interface. The maximum displacement of the interface at one time-step should be less than the minimum grid spacing.
IV. RESULTS AND DISCUSSION
Ruijk Rui⫹1jk v Rijk v Rij⫹1k w Rijk w Rijk⫹1 c Rijk
␦pijk ⫽
As mentioned earlier, the amount of computation required to solve this problem is tremendous. The specific reasons are as follows.
Hence, the code has been parallelized and optimized on the Origin 2000 supercomputer platforms at the National Center for Supercomputing Applications at Urbana-Champaign, IL.
冢冣 冢冣 ␦uijk ␦ui ⫹1jk ␦vijk ␦vij⫹1k ␦wijk ␦wijk⫹1 ␦pijk
F. Vectorization
(1) In this study, material properties across the L/V interface are smoothed out using an interpolation function for numerical purposes with some degree of discontinuity preserved. (2) Material properties for the high-temperature range are extrapolated. (3) Plasma is not taken into account. (4) The gas flow is assumed incompressible for the sake of simplicity. (5) It is assumed that the evaporating mass flux does not interact with the hole surfaces, i.e., the recondensation phenomenon is neglected. In order to demonstrate the influence of each surface phenomenon on the melt pool, four different cases have been considered for simulation, as shown in Table III. Starting from a conduction-mode welding, which is the case A in Table III, evaporation, free surface evolution, and multiple reflections effects will be added to the simulation one-byone. In case A, evaporation, free surface evolution, and multiple reflections are all turned off. In case B, only the evaporation switch is turned on but the free surface is frozen in order to obtain a pressurized melt pool without a keyhole. In case C, the free surface evolves according to the recoil pressure and melt flow, but multiple reflection phenomena are not simulated. In case D, all three free surface phenomena VOLUME 33A, JUNE 2002—1825
Table II. Material Properties for Iron ( pcr is obtained by substituting Tcr into the CC relation; for simulation, material properties are extrapolated to the critical point) Property
Symbol
Value
Melting temperature Normal boiling temperature Critical point temperature (Beautl 1994) Critical point pressure Liquid density Solid density Kinematic viscosity Surface tension Latent heat of vaporization Latent heat of fusion Solid thermal conductivity Liquid thermal conductivity Liquid constant-pressure specific heat Solid constant-pressure specific heat Liquid thermal diffusivity Solid thermal diffusivity Laser absorptivity for flat surface Liquid enthalpy (Beautl 1994)
Tm Tb
1809.0 (K) 3133.0 (K)
Tcr
9250.0 (K)
pcr l s vl Lv
8.973 ⫻ 108 (Pa) 6518.5 (kg/m3) 7870.0 (kg/m3) 4.936 ⫻ 10⫺7 (m2/sec) surface tension curve 6.3639 ⫻ 106 (J/kg)
Lm ks
2.7196 ⫻ 105 (J/kg) 40.96 (W/m K)
kl
43.99 (W/m K)
Cpl
804.03 (J/kg K)
Cps
658.63 (J/kg K)
␣l ␣s A0
8.39 ⫻ 10⫺6 (m2/s) 7.90 ⫻ 10⫺6 (m2/s) 0.1
hl
(⫺0.18626 ⫹ 8.2516 ⫻ 10⫺4 * T ) * 1 ⫻ 106
Table III. Four Cases Studied in This Article Simulation Identification
Evaporation
Free Surface Evolution
Multiple Reflections
A B C D
no yes yes yes
no no yes yes
no no no yes
are simulated. In this way, we can have a better understanding of the role of each component systematically. The same parameters are used for all cases. As an energy source, a continuous-wave 4 kW CO2 Gaussian laser beam with a 500-m beam diameter is considered. Beam profile is assumed constant along the z direction, i.e., the depth-offocus effect is not taken into account. For the bean scanning speed, 60 ipm is used. For the target material, a steel coupon with a thickness of 1.214 mm is selected. Material properties are listed in Table II. A laser-beam absorptivity of 0.1 is assumed. A. Case A Figure 7 presents the flow field when evaporation, free surface evolution, and multiple reflections are all neglected even though the interface temperature is well beyond the normal boiling point. In this case, thermocapillary force is the only driving force for the flow. As seen, the melt flow structure is very similar to that of conduction-mode welding even though a high 1826—VOLUME 33A, JUNE 2002
Fig. 7—Case A: conduction mode welding (units in meters, and laser beam scanning to the right).
power laser is assumed for this simulation. Two recirculation zones are found along the symmetry line. In fact, in three dimensions, there exists only one recirculating pool of a torus shape. The melt pool is very small and shallow, and its length is less than 1 mm. Figure 7(2) presents the velocity field of liquid metal and gas. Just like the melt flow, the gas flow in this case is solely induced by the thermocapillary force at the interface. Overall, the flow is directed downward and then turns to the direction of the thermocapillary force near the L/V interface. At the solid/liquid/vapor interfacial point (in fact, this is a circle in three dimensions), the gas flow forms a recirculation zone at the front; at the back, however, the flow fails to form a zone and directs further backward due to the beam scanning effect. B. Effect of Evaporation: Case B In this case, the effect of evaporation is studied while the free surface is frozen. In Figure 7, recoil pressure is not incorporated in the flow field even though the L/V interface temperature is much higher than the normal boiling point. Figure 8 demonstrates how the flow field, both melt and gas flows, changes dramatically when the recoil pressure is incorporated. Now the directions of the melt and gas flows are completely reversed. Due to the strong recoil pressure, both the melt and gas flows are very intense and start from the L/V interface. The vaporizing gas speed is much higher than that of melt because of the huge density difference between the two fluids. It is also clearly seen that the melt pool becomes METALLURGICAL AND MATERIALS TRANSACTIONS A
Fig. 8—Case B: evaporation-generated recoil pressure is applied but the free surface is frozen (units in meters, and laser beam scanning to the right).
much longer and wider than that in case A, owing to the recoil pressure. It is apparent that recoil pressure is responsible for the big melt pool in high-power laser processings of metals. In other words, we can conclude that the thermocapillary force is relatively insignificant compared to the evaporation-generated recoil pressure. Since the L/V interface is artificially frozen, the flow structure shown in the figure is somewhat unrealistic in terms of mass balance. Especially at the interface where evaporation occurs, both the melt and gas head to opposite directions and conservation of mass is not satisfied.
Fig. 9—Case C: velocity field with self-evolving L/V interface shape. Multiple reflections are not simulated (units in meters, and laser beam scanning to the right).
C. Effect of Free Surface Evolution: Case C In this case, free surface evolves self-consistently as a result of fluid flow and evaporation, still neglecting multiple reflections. Since the interface is moving dynamically instead of being frozen intentionally, conservation of mass can now be fully achieved even at those regions of strong evaporation. Figure 9 shows the velocity field seen from various viewpoints. Figure 9(3) is a combined plot of melt and gas flows. Overall, the gas flow is dominated by the intensive evaporating flux, which is much stronger than the melt flow. This evaporating gas flux creates a keyhole in the melt pool. A big recirculating region is observed behind the keyhole. Figure 10 presents the intensity field and temperature field. As expected, the intensity profile (Figure 10(1)) is spatially Gaussian on a curved keyhole surface. The maximum intensity occurs right at the beam center location and is about 4 ⫻ 105 W/cm2, which is same as the value for a flat surface. Figure 10(2) is the temperature profile on the L/V interface. The maximum temperature is about 3600 K and is located almost right at the Gaussian laser-beam center. As shown in this section, the effect of free surface evolution on the melt flow is substantial. The melt flow field deviates significantly from the conduction-mode welding and falls into a totally different regime. METALLURGICAL AND MATERIALS TRANSACTIONS A
D. Effect of Multiple Reflections: Case D In this case, the multiple reflections effect is considered together with all the previous interface phenomena, as proposed in this article. Since the purpose of this section is to make a comparison between the simulations without and with multiple reflections, we present the plots of velocity, intensity, and temperature fields taken at the same time as in Section IV–C. More in-depth and detailed study of the full simulation will be dealt with in the Part II of this work.[30] Figure 11 shows the velocity field seen from various viewpoints. First of all, it is apparent that, with multiple reflections, the hole diameter has become larger in size and the keyholing or penetration speed has increased substantially. The gas flow is much more intense and occurs over a broader area. These drastic changes by multiple reflections can be better explained in intensity and temperature plots in Figure 12. As shown in Figure 12(1) and (2), the maximum intensity and maximum temperature on the L/V interface have risen to 106 W/cm2 and 4200 K from 4 ⫻ 105 W/cm2 and 3600 K, respectively, due to multiple reflections. In fact, the effective laser absorptivity of keyhole has increased to 0.227 from the initial value of 0.1, which is the our assumed initial value for a flat surface. The intensity profile starts to deviate VOLUME 33A, JUNE 2002—1827
Fig. 11—Case D: velocity field with self-evolving L/V interface shape. Multiple reflections are simulated (units in meters, and laser beam scanning to the right). Fig. 10—Case C: laser intensity profile and temperature profile on the selfevolving keyhole surface. Multiple reflections are not simulated.
from the original Gaussian profile. In Figure 12(1), it is observed that, with multiple reflections, the laser-beam energy is more highly concentrated near the center area. This is in line with authors’ study of laser drilling with multiple reflections.[29] In laser welding, however, it is only valid when the hole is shallow. The keyhole shown in this figure is not deep enough to demonstrate the more dynamic evolution of temperature and intensity distributions by multiple reflections. More thorough investigation of intensity, temperature, and velocity distributions with keyhole propagation can be found in the companion article.[30] V. CLOSURE A self-consistent three-dimensional laser-keyhole welding model is presented, which simulates movements of L/V and S/L interfaces using the level set method and mixture continuum model. The SCGS method[45] is adopted. At the evaporating node points, the recoil pressure values calculated in the evaporation model replace the flow-generated pressure values when updating the pressure field in the SCGS method. In this way, the recoil pressure could be incorporated into the flow solution. By adopting the level set method and 1828—VOLUME 33A, JUNE 2002
mixture continuum model, the authors have derived the momentum and energy equations, which are valid for gas, liquid, and solid phases. All the L/V and S/L interface boundary conditions, such as the thermocapillary, capillary, recoil pressure boundary conditions, and the Stefan’s conditions, are successfully applied to the moving L/V and S/L interfaces. In this study, each individual process associated with L/V interface, such as free surface evolution, evaporation, homogeneous boiling, and multiple reflections has been modeled. The narrow-band level set method has been implemented to capture free surface movements consistently. In the evaporation model, Knight’s jump conditions[31] are introduced to model the kinetic Knudsen layer at the evaporating surface. The homogeneous boiling phenomenon has also been considered by estimating the liquid superheat and by extrapolating the surface tension curve to the critical point. A method to simulate multiple reflections for arbitrary three-dimensional surfaces has been developed using the level set values. In addition, this article presents simulation results regarding the effect of each L/V interface phenomenon and studies how and why the simple conduction-mode welding can transform to the complex keyhole-mode welding. It has been revealed that each and every interface phenomenon play a METALLURGICAL AND MATERIALS TRANSACTIONS A
Fig. 12—Case D: laser intensity profile and temperature profile on the selfevolving keyhole surface. Multiple reflections are simulated.
critical role in the laser-keyhole welding process, so any oversimplified models can result in erroneous and unrealistic predictions. Detailed simulation results with practical welding parameters and experimental verification to the model will be presented in Part II of this study.[30] ACKNOWLEDGMENTS This work was made possible by the continued support of the Office of Naval Research under Grant No: N0001497-1-0124. Dr. George Yoder is the program manager. The authors also acknowledge the National Center for Supercomputing Applications (NCSA). REFERENCES 1. D. Adalsteinsson and J.A. Sethian: J. Comput. Phys., 1995, vol. 118, pp. 269-77. 2. D. Adalsteinsson and J.A. Sethian: J. Comput. Phys., 1995, vol. 120, pp. 128-44. 3. D. Adalsteinsson and J.A. Sethian: J. Comput. Phys., 1995, vol. 122, pp. 348-66. METALLURGICAL AND MATERIALS TRANSACTIONS A
4. S. Asai and I. Muchi: Trans. Iron Steel Inst. Jpn., 1978, vol. 18, pp. 90-98. 5. V. Batanov, F. Bunkin, A. Prokhorov, and V. Fedorov: Sov. Phys. JETP, 1973, vol. 36 (2) pp. 311-22. 6. W. Bennon and F. Incropera: Int. J. Heat Mass Transfer, 1987, vol. 30 (10), pp. 2161-70. 7. W. Bennon and F. Incropera: Int. J. Heat Mass Transfer, 1987, vol. 30 (10), pp. 2171-87. 8. W. Bennon and F. Incropera: Num. Heat Transfer, 1988, vol. 13, pp. 277-96. 9. R.L. Burden and J.D. Faires: Num. Analysis. PWS Publishing Company, Boston, MA, 1993. 10. V.P. Carey: Liquid-Vapor Phase Change Phenomena, Hemisphere Publishing Corporation, Bristol, PA, 1992. 11. C. Chan, J. Mazumder, and M. Chen: Metall. Trans. A, 1984, vol. 15A, pp. 2175-84. 12. C. Chan, J. Mazumder, and M. Chen: J. Appl. Phys., 1988, vol. 64 (11), pp. 6166-74. 13. S. Chen, B. Merriman, S. Osher, and P. Smereka: J. Comput. Phys., 1997, vol. 135, pp. 8-29. 14. V. Craciun, D. Craciun, C. Boulmer-Leborgne, and J. Hermann: Phys. Rev. B, 1998, vol. 58, p. 6787. 15. P.G. Debenedetti: Metastable Liquids Concepts and Principles, Princeton University Press, Princeton, NJ, 1996. 16. C. Domb: The Critical Point, Taylor & Francis, London, 1996. 17. J. Dowden, M. Davis, and P. Kapadia: J. Fluid Mech., 1983, vol. 126, pp. 123-46. 18. J. Dowden and P. Kapadia: J. Phys. D: Appl. Phys., 1995, vol. 28, pp. 2252-61. 19. R. Fabbro and K. Chouf : J. Appl. Phys., 2000, vol. 87 (9), pp. 4075-83. 20. F. Gagiliano and U. Paek: Appl. Opt., 1974, vol. 13 (2), pp. 274-79. 21. F. Hensel: J. Phys.: Condens. Matter, 1990, vol. 2, pp. SA33-SA45. 22. F. Hensel and P. Edwards: Phys. Worlds, 1996, Apr., pp. 43-46. 23. F. Incropera, A. Engel, and W. Bennon: Num. Heat Transfer, Part A, 1989, vol. 16, pp. 407-27. 24. A. Kaplan: J. Phys. D: Appl. Phys., 1994, vol. 27, pp. 1805-14. 25. A. Kar, T. Rockstroh, and J. Mazumder: J. Appl. Phys., 1992, vol. 71 (6), pp. 2560-69. 26. R. Kelly and A. Miotello: Appl. Surf. Sci., 1996, vol. 96–98, pp. 205-15. 27. R. Kelly and A. Miotello: Phys. Rev. E, 1999, vol. 60, p. 2616. 28. H. Ki: Ph.D. Thesis, University of Michigan, Ann Arbor, MI, 2001. 29. H. Ki, P. Mohanty, and J. Mazumder: J. Phys. D: Appl. Phys., 2001, vol. 34, pp. 364-72. 30. H. Ki, P.S. Mohanty, and J. Mazumder: Metall. Mater. Trans. A, 2002, vol. 33A, pp. 1831-42. 31. C.J. Knight: AIAA J., 1979, vol. 17 (5), pp. 519-23. 32. J. Kroos, U. Gratzke, and G. Simon: J. Phys. D: Appl. Phys., 1993, vol. 26, pp. 474-80. 33. K.N. Lankalapalli, J.F. Tu, and M. Gartner: J. Phys. D: Appl. Phys., 1996, vol. 29, pp. 1831-41. 34. M. Martynyuk: Russ. J. Phys. Chem., 1975, vol. 49 (10), pp. 1545-47. 35. M. Martynyuk: Russ. J. Phys. Chem., 1983, vol. 57 (4), pp. 494-501. 36. A. Matsunswa and V. Semak: J. Phys. D: Appl. Phys., 1997, vol. 30, pp. 798-809. 37. J. Mazumder: Opt. Eng., 1991, vol. 30 (8), pp. 1208-19. 38. J. Milewski and E. Sklar: Modeling Simul. Mater. Sci., 1996, vol. 4, pp. 305-22. 39. W. Mulder, S. Osher, and J.A. Sethian: J. Comput. Phys., 1992, vol. 100, pp. 209-28. 40. S. Osher and J.A. Sethian: J. Comput. Phys., 1988, vol. 79, pp. 12-49. 41. S.V. Patankar: Numerical Heat Transfer and Fluid Flow, Taylor & Francis, 1980. 42. C. Prakash and V. Voller: Num. Heat Transfer, Part B, 1989, vol. 15, pp. 171-89. 43. P. Prescott, F. Incropera, and W. Bennon: Int. J. Heat Mass Transfer, 1991, vol. 34 (9), pp. 2351-59. 44. R. Ross and D. Greenwood: Progr. Mater. Sci., 1969, vol. 14, pp. 175-242. 45. P. Sathyamurthy: Ph.D. Thesis, University of Minnesota, Twin Cities, MN, 1991. 46. V.V. Semak, W.D. Bragg, B. Damkroger, and S. Kempka: J. Phys. D: Appl. Phys., 1999, vol. 32, pp. L61-L64. 47. J.A. Sethian: Level Set Methods and Fast Marching Methods, 2nd ed., Cambridge University Press, Cambridge, United Kingdom, 1999. 48. J.A. Sethian and J. Strain: J. Comput. Phys., 1992, vol. 98, pp. 231-53. VOLUME 33A, JUNE 2002—1829
49. V. Skripov: Metastable Liquids, John Wiley & Sons. New York, NY, 1974. 50. V. Smiljanovski, V. Moser, and R. Klein: Combust. Theory Modelling, 1997, vol. 1, pp. 183-215. 51. P. Solana and G. Negro: J. Phys. D: Appl. Phys., 1997, vol. 30, pp. 3216-22. 52. P. Solana and J. Ocana: J. Phys. D: Appl. Phys., 1997, vol. 30, pp. 1300-13. 53. K. H. Song and X. Xu: Appl. Surf. Sci., 1998, vol. 127-129, pp. 111-16. 54. W. Sudnik, D. Radaj. S. Breitschwerdt, and W. Erofeew: J. Phys. D: Appl. Phys., 2000, vol. 33, pp. 662-71. 55. W. Sudnik, D. Radaj, and W. Erofeew: J. Phys. D: Appl. Phys., 1996, vol. 29, pp. 2811-17.
1830—VOLUME 33A, JUNE 2002
56. M. Sussman, P. Smereka, and S. Osher: J. Comput. Phys., 1994, vol. 114, pp. 146-59. 57. S. Vanka: J. Comput. Phys., 1986, vol. 65, pp. 138-58. 58. H. Versteeg and W. Malalasekera: An Introduction to Computational Fluid Dynamics, Longman, England, 1999. 59. P. Wei and C. Ho: Int. J. Heat Mass Transfer, 1998, vol. 41, pp. 3299-3308. 60. B. S. Yilbas and M. Sami: J. Phys. D: Appl. Phys., 1997, vol. 30, pp. 1996-2005. 61. T. Ytrehus and S. Østmo: Int. J. Multiphase Flow, 1996, vol. 22 (1), pp. 133-55. 62. R.L. Zehr: Ph.D. Thesis, University of Illinois at Urbana–Champaign, Urbana, IL, 1991.
METALLURGICAL AND MATERIALS TRANSACTIONS A
