AN APPROACH FOR INTERFACE CONDITION OF PHASE-CHANGE HEAT CONDUCTION IN CURVILINEAR COORDINATES

Phase change materials are vastly used in thermal engineering applications. The model studies reduce the experimental time and cost and gives insight into the physical process and and provides relation between the process outcomes and the influencing parameters on the process. One of the challenges in the model study related to the phase change problem is setting the appropriate boundary conditions across the phases. This is because of the fictitious definition of the mush zone across the phases. This situation becomes complicated when setting the boundary conditions across the odd geometric shapes. In this study, mathematical formulation of the condition for energy-balance at the interface of the phase changing is investigated using the curvilinear coordinate system without requiring the coordinate system. The proposed arrangement enables to create a curvilinear system via transformation equations from another curvilinear coordinate system. It also provides mathematical formulation of the interfacial boundary conditions across the phases.


INTRODUCTION
Interfacial conditions for the phase changing materials remain important in terms of computational efforts and correct description of the physical process involved during the phase change. Considerable research studies were carried out to examine phase change problem and the interfacial conditions. The turbulenceinduced interfacial instability in two-phase flow with moving interface was studied by Balabel [1]. He introduced the transition from one phase to another through incorporating a consistent balance of kinematic and dynamic conditions on the interface separating the two phases. However, the topological changes of the interface were formulated via adopting the level set approach. In this case, through using the interfacial markers on the intersection points, the interfacial stresses and the interfacial driving forces were estimated. This arrangement allowed prodicting the normal interface velocity, which could be extended to the higher dimensional level set function and used for the interface advection process. The interfacial conditions incorporating Stefan boundary at solid-liquid interface was examined by Turkyilmazoglu [2] after using the single and double phase models. The findings revealed that the physical phase transition process took place at a constant speed due to the imposed movement of the material along or reverse directions. In addition, the coefficients determining the movement of the phase change interface could be presented analytically. The interfacial moving boundary problem occurred in chemical processes, such as combustion. The conservation equations and constitutive equations, within the framework of moving interface, could be incorporated to formulate the interfacial conditions. Coordinate transformation should be incorporated towards achieving a fixed interface formulation from the moving interface problem, which resulted in a fixed domain of each phase ranging from 0 to 1. The interfacial boundary conditions and residual trapping in relation to wetting phase flow was investigated by Heshmati and Piri [3]. They considered the effect of changes in invading wetting phase flow rate and injection of a non-wetting droplet on pore fluid configuration. They demonstrated the local perturbations of non-wetting phase in the flow system. A sharp-interface level-set method for the phase change interfacial conditions was presented by Lee and Son [4] towards simulating growth and collapse of a compressible vapor bubble. However, the interface tracking method was extended including the influence of bubble compressibility and liquid-vapor phase change. They used the ghost fluid method implementing the matching conditions of velocity, stress and temperature at the interface. In phase change problems, such as evaporation or condensation, there were two-way coupling of momentum, heat and species transfer took place. In formulating the interfacial conditions for phase change problem, such as cavitation problem, a sharp interface approach incorporating the volume-of-fluid and ghost-fluid methods was demonstrated to be very effective Michael et al. [5]. However, the level set function should be constructed from the volume fraction function. A marching cubes method could be used determining the interface area at the interface grid cells. A parametric study exploring the approximation of the exact analytical solution of the Stefan problem in a finite phase changing layer was carried out by Mazzeo and Oliveti [6]. They have undertaken the exact analytical solution determining the sensible and latent stored energy in the layer through incorporating the implicit transcendental equation with complex thermal parameters. In this case, the dimensionless solution was demonstrated to be a function of the Fourier and Stefan numbers for the two phases. The analytical and numerical analysis for the solidification of the phase change material inside a rectangular finned container with time-dependent boundary condition was untaken by Taghilou and Talati [7]. They introduced two special cases of linear and sinusoidal time-dependent boundary conditions and presented analytical results for non-dimensional temperature distribution and position of the solidification front. The findings were compared with those of the lattice Boltzmann method with respect to the Fourier number at which the Stefan number was fixed. The gradient augmented level set method for phase change simulations was introduced by Anumolu and Trujillo [8]. They demonstrated that the gradient augmented level set method enabled to obtain sharp capturing of the vaporization process. This arrangement provided i) identification of the vapor-liquid interface at the subgrid level, ii) discontinuous treatment of thermal physical properties (except for viscosity), and iii) enforcement of mass, momentum, and energy jump conditions, where the gradients of the dependent variables were obtained at interface. The solution of the two-phase Stefan problem using analytical method was presented by Khaled et al. [9]. The Eigen-conditions pertinent to the governing equations were formulated using the separation of variable technique. The eigenvalues were obtained by applying the boundary conditions for liquid and solid phases. They noted that the radial eigenvalues were free from imaginary values and interface equation were solved and analyzed through varying the Stefan number.
On the other hand, heat conduction problems related to a phase change from liquid to solid phases or vice versa require an energy balance equation to be satisfied at the liquid-solid interface with the condition that the temperature at the interface remains the melting temperature. Moreover, to introduce the interface conditions in computations while satisfying the energy-balance, it is necessary to modify the conditions appropriate to the particular applications. Some of these modifications were reported in the previous studies [1 -12]. Although the interface conditions were provided separately either in the Cartesian or in the Curvilinear coordinate system in the previous studies [1 -12], the conditions describing the interface of the phase changing materials were not generalized for any coordinate system. In the present study, the condition of energy-balance at the phase change interface is presented comprehensively without referencing to any particular coordinate system. In addition, the proposed arrangement allows to generate a curvilinear coordinate system through transformation equations from another curvilinear coordinate system.

MATHEMATICAL ANALYSIS
Consider the heat conduction problem in which there is a phase change taking place from liquid to solid or vice versa at a single temperature in a pure substance. At the liquid-solid interface two conditions have to be satisfied, where is the liquid/solid temperature, is the melting temperature, ″ is the liquid/solid heat flux vector evaluated at the solid/liquid interface, is the density of the solid phase, is the latent heat of melting, is the 'coordinate velocity' at the interface [5], is a vector normal to the interface, but not necessarily of unit magnitude, and the subscripts and denote the liquid and solid phases, respectively.
The energy-balance at the interface condition (Equation 2) is symbolic and it is not suitable for use in computations; therefore, various researchers have modified this condition so as to be of the practical applicability [10 -12]. One such result was presented by Patel [10] for the Cartesian coordinate system. Ozisik [12] describes the interface condition in the cylindrical coordinate system. The concern is that the interface conditions described in [10,12] solely depend on the variables of the particular coordinate system, the Cartesian or cylindrical, without referencing to any other coordinate system. The energy-balance interface condition, Equation 2, for a general, non-orthogonal, curvilinear coordinate system is derived earlier [12]. Although the derivation is general, it has a shortcoming of incorporating transformation equations explicitly. In this case, the energy-balance interface condition involves the transformation equations from the Cartesian coordinate system to the non-orthogonal, curvilinear coordinate system. This implicitly demonstrates that their interface condition has the short come of the generality of those described in [10,12]. This paper is an attempt to address this shortcoming by restating Equation 2 in a mathematical form that only makes reference to the curvilinear coordinate system for which it has been derived. A further benefit of restating the energy-balance interface condition in this form is that the curvilinear coordinate system can be generated from another curvilinear coordinate system and not solely from the Cartesian coordinate system. In the following analysis Tensor calculus is to be used and the final result incorporates the metric tensor. The range and summation conventions are to be employed as used in Tensor Calculus.
Let ( , , ) be the coordinates in a general, non-orthogonal, curvilinear coordinate system that is generated from the Cartesian system by means of the transformation equations, One can use the following terminology so as to be able to employ the methods of Tensor Calculus, ( , , ) → ( 1 , 2 , 3 ) or simply (4) ( , , ) → ( 1 , 2 , 3 ) or simply (5) The liquid-solid interface is in general a three-dimensional surface and in a curvilinear coordinate system it can be defined as, The interface can also be defined parametrically as, where and are the surface coordinates. In the present analysis we will first define the interface by means of Equation (6b) and later we will switch to the form defined by Equation (6a). Let = ( ( , , )) be the position vector of a point on the interface in the curvilinear coordinate system. Then two, linearly independent tangent vectors on the interface are given as, where are the covariant bases vectors in the general curvilinear coordinate system. The normal vector to the surface is written as [13], where is the Levi-Civita symbol, = | | is the determinant of the metric tensor, is the e-permutation symbol and are the contravariant bases vectors in the general, curvilinear, coordinate system. The coordinate velocity of the surface is given as [14], Also, one can note that for a pure substance, the solid-liquid interface is an isotherm. This implies that as one moves from one point to another on the interface, the temperature remains constant and is equal to the melting temperature of the substance. Since the position on the interface is defined by the surface coordinates and therefore, the partial derivatives of the solid/liquid temperatures with respect to the coordinates and vanishes. Therefore, from Equation (10) we get, and, One can now restrict the analysis to interfaces that can be described by equation (6a), 3 = 3 ( 1 , 2 , ). This can be achieved by letting = 1 and = 2 in equations (12), (13). Furthermore, since 3 = 3 ( , , ), therefore, now the coordinate 3 on the surface will be a function of the coordinates 1 and 2 as well as that of time , i.e. 3 = 3 ( 1 , 2 , ). From Eq. (12) one can obtain, and from Eq. (13) we obtain, Carrying out the implied summations in Eq. (11b); noting that now Equation (16) has been derived by using Equation (3) as the transformation equations and these involve the Cartesian coordinates, however, the final mathematical form of Equation (16) is completely independent of the Cartesian coordinate system and involves no reference to the transformation Equation (3). The components of the metric tensor appearing in Equation (16) are themselves functions of the coordinate variables ( , , ). This mathematical form of the interface condition is in contrast to the one derived in [12] and has the character of those described in [11].
An important implication of the particular mathematical form of Equation (16) is that it is now not necessary that the curvilinear coordinate system with coordinates ( , , ), be generated from a Cartesian coordinate system by means of Equation (3). It can now be generated from another curvilinear coordinate system with the coordinates (̄,,̄) and the mathematical form of the interface condition Equation (16) remains unchanged. This is again in contrast to the interface condition described in [12]. To further elaborate this point, consider the transformation equation between two curvilinear coordinate systems, ̄=̄( , , ), ̄=̄( , , ), ̄=̄( , , ) The metric tensor so obtained can be substituted into Equation (16) to obtain the interface condition. Finally, if the curvilinear coordinate system is orthogonal then the interface condition simplifies to,

APPLICATION TO AN AXISYMMETRIC LASER HEATING OF A SOLID SUBSTRATE
Consider the problem of the axisymmetric, laser heating and subsequent melting of a solid substrate in line with the previous study [15]. Since the problem is axisymmetric therefore, = ( , , ) , = ( , , ) (23)

Figure 1. Schematic of an axi-symmetric laser melting process
The numerical solution of this problem may require the generation of body-fitted grids in the irregular liquid and solid regions as shown in figure 1. Now, for the convenience of the numerical solution it may be necessary to preserve the axisymmetric nature of the problem. This problem is then essentially equivalent to generating a curvilinear coordinate system from another curvilinear (cylindrical) coordinate system. In other words, the grid generation involves Equation (18) We make the following designations, The components of the metric tensor in the cylindrical coordinates are, 1 1 =3 3 = 1,2 2 = 2 and ̄= 0, ≠ The components of the metric tensor in the curvilinear coordinates are calculated from Equation (20) as, The contravariant components of the metric tensor can then be determined from the matrix inverse of , Finally, the interface condition for the axisymmetric, laser heating problem can be formulated from Equation (16).
Since the heating problem considered is associated with the laser interaction of metallic substrate (figure (1)), the output power intensity distribution at the workpiece remains important. In this case, it is assumed to be a Gaussian and its spot centre is at the centre of the co-ordinate system is considered ( figure (1)). This arrangement results in an axisymmetric heating of the substrate material. In the initial stage of the laser heating, the conduction heating of the solid substrate with insulated boundary condition at the surface can be considered. This is due to that the pulse duration is short and the convective and radiate losses from the surface are negligibly smaller than the internal energy gain of the substrate material during heating [16]. Consequently, the heat transfer equation for a solid phase heating due to a laser irradiation pulse with a Gaussian intensity profile can be written as: , , and are the laser peak power intensity, absorption coefficient, reflectivity and the Gaussian parameter, respectively. Since the heating problem is transient, the initial condition should be defined, i.e. initially it is assumed that the substrate material is at a uniform temperature, which can be specified. Therefore: at time zero ⇒ = 0: ( , , 0) = (specified). In order to solve Equation 30, two boundary conditions for each principle axis should be specified. Due to the short duration of laser pulse, insulated boundary is assumed at the surface and at a distance considerably away from the surface (at infinity), it is also assumed that the heating has no effect on the temperature of the substrate material; consequently, at a depth of infinity, temperature is assumed to be constant and equals to the initial temperature of the substrate material. The boundary conditions, therefore, are: at infinity is = ∞: ( , ∞, ) = (specified) and at infinity is = ∞: (∞, , ) = (specified). However, at symmetry axis, it is = 0: (0, , ) = 0 and at the substrate surface it is = 0: ( ,0, ) = 0.
In the case of melting and evaporation, phase change takes place. Moreover, the functional relation of pressure dependence of boiling temperature is not known for the metals, which are used in industry; therefore, it is assumed that the substrate material has single melting and boiling temperatures. To accommodate the phase change process, heat transfer equations for solid heating should be modified. In this case, an energy or enthalpy methods can be used. In the enthalpy method, the governing equation of energy transport can be written in terms of enthalpy equation [16] and interfacial condition introduced in Equation 22 can be incorporated.
In order to solve Equation 30 with the appropriate boundary conditions, the numerical code is developed to incorporate the governing equations of heat transfer and the interface condition (Equation 28) to obtain the interface velocity of the phase chaining substrate material. The short description of the numerical method is as follows: To discretize the governing equations, a finite difference scheme is introduced. Figure 2 shows the grid used in the numerical analysis. The details of the numerical scheme are given in [17]. To compute the equations discretized for temperature field and relative positions of solid-liquid and liquid-vapor interface, an implicit scheme is used, i.e. using the initial conditions, the temperature in the whole domain is calculated for following time steps with the respective conditions. The calculation domain is divided into grids and grid independence test is being performed for different grid size and orientation. A in home computer program based on implicit scheme is developed to compute the temperature field.

RESULTS AND DISCUSSION
Interfacial boundary condition for phase changing material is formulated and its application relevant to laser melting and evaporation under the pulse heating of solid substrate is presented. Stainless steel is considered as the phase changing material. The spatial and temporal behavior of vapor and liquid fractions are computed. This indicates that at the liquid phase is present on the liquid surface. However, as the heating due to laser pulse irradiation progresses, the fraction of the liquid-vapour ( ) increases sharply in the region close to surface (0 ≤ ≤ 9.7 × 10 −9 m). The rapid increase of ( ) is related to the laser beam irradiation, which increases rapidly with progressing time in the surface region, up to the depth ( = 1 , where  is the laser beam absorption depth) below the surface. However, as the depth below the surface increases further, the absorbed power reduces significantly and evaporation ceases. It should be noted the source term resembling the laser power in Equation 29 reduces exponentially with the depth according to the Lambert's Beer law. In this case, as the depth increases below the surface, beam energy absorbed by the substrate material reduces exponentially. This in turn modifies the interfacial conditions for liquid-vapor front; in which case, fraction of liquid-vapor changes drastically and size of the mush zone changes. Despite the fact that the energy gain via the absorption of the laser irradiation reduces at a depth = 1 , the energy transport via conduction (diffusion) from the surface region to the depth = 1 below the surface contributes to the phase change process while influencing the size of the mushy zone. The change of the radial location from the symmetry axis to 2 ( being the irradiated spot radius) alters the fraction of the liquid-vapor. Hence, the rise of the fraction of liquid-vapor becomes smaller at = 2 than that corresponding to the symmetry axis. The time shift for the increase of the fraction of the liquidvapour at = 2 is associated with the time required for temperature increasing the boiling temperature of the substrate material in the liquid phase. In this case, fast change of the fraction of liquid-vapor in the region of 1.75 ns from the laser pulse beginning at the locations z = 0 and r = ro/2 and z = 1/ and r = 0 is because of the rapid heating of the liquid -vapor mushy zone in the early pulse heating. However, when the heating period increases, the mushy zone size enlarges and the change of fraction (xb) becomes less. Since the laser heating pulse at the surface along the symmetry axis (r = 0) is higher than that of the radial location r = ro/2, the mush zone depth becomes larger along the symmetry axis as compared to that of r = ro/2. Figure 4 depicts the variation of fraction of solid-liquid ( ) phases in the mushy zone at locations along the x-axis in the substrate material for two different radial locations. The increase of fraction of solidliquid is fast, which is more pronounced along the symmetry axis ( = 0) in the region close to the surface. This behavior is related to the absorption of the incident laser radiation by the substrate material in this region, which is considerably high in accordance with the Lambert's Beer law. Increasing depth lowers the fraction of the solid-liquid interface; however, along the symmetry line this reduction is less as compared to other regions. This attributed to the heat diffusion from the surface towards the solid bulk along the symmetry axis, which remains high as compared to other regions. The increase of the fraction of solid-liquid phases becomes about same at various depths below the surface, i.e. high rate of thermal energy is carried through diffusional heat transfer into the solid bulk. However, for the radial location at = 2 , the fraction solid-liquid becomes similar to that of its counterpart along the symmetry axis; nevertheless, the increase of the fraction remains smaller for = 2 than that corresponding to the symmetry axis. This is due to the fact that the laser irradiated energy absorbed at location = 2 remains less than corresponding to the symmetry axis ( = 0) because laser intensity distribution at the surface, which is Gaussian. Figure 4. Solid-liquid quality with time at various locations below depth (z-axis) and radial location is the symmetry line (r = 0) Figure 5 depicts the surface recession velocity along the radial axis for various heating times. As the mush zone thickness approaches zero, the solid surface appears below the liquid layer; therefore, the recession of the solid surface with time is considered to be a recession velocity of the surface. During the initial irradiation of the laser pulse, the recession velocity becomes high, but the size of the surface recessed along the radial direction becomes small. When the laser pulse heating progresses, melt zone along the radial direction becomes large; this in turn causes rapid rise of the recession velocity in this period. When the laser heating progresses further ( ≥ 15.3ns), the surface recessed enhances along the radial direction despite the fact that the velocity of the recession decreases. The enlargement of the surface recessed along the radial direction is related to the absorbed irradiated energy from the laser pulse, i.e. the laser power intensity distribution at the surface follows the Gaussian profile. Figure 6 shows the behavior of the evaporation velocity of the surface at locations along the radial direction. Temporal behavior of the vapor front recession velocity becomes almost similar to that of the temporal behavior of the laser pulse. This behavior remains true for all the locations along the radial directions. However, increasing the radial location to = 2 , the recession velocity decay becomes different than that corresponding to other locations along the radial directions, particularly for the heating periods of ≥ 1.5 × 10 −8 s. This is attributed to the heat conduction in the radial direction, which remains small when the radial location increases away from the center of the irradiated spot, i.e. the irradiated energy absorbed from the laser power towards the phase change becomes less in this region, which in turn lowers the recession of the solid surface.

CONCLUSION
The general form of the energy-balance interface condition that is applicable in a non-orthogonal, curvilinear, coordinate system is derived in this paper. The mathematical form presented here is such that it involves no reference to any other coordinate system which might have been used in its generation. The formulation introduced for the interface condition can be used in situations where the curvilinear coordinate system is generated through transformation equations from another curvilinear system and not exclusively from a Cartesian system. The paper also presents the formulation of the energy-balance as an interface condition for the special case of an orthogonal coordinate system while demonstrating the equation satisfying the energy balance at the interface. The interfacial condition developed from the mathematical analysis is incorporated in the high power laser pulse heating of steel surface. The recession velocity of evaporation front and the fraction of the liquid-solid interface are predicted via using in the home developed computer code.