Numerical Analysis of Melting/Solidification Phenomena Using a Moving Boundary Analysis Method X-FEM

Abstract A numerical analysis method for melting/solidification phenomena has been developed to evaluate feasibility of the several candidate techniques in the nuclear fuel cycle. Our method is based on the extended finite element method, which has been used for moving boundary problems. The basic idea of the extended finite element method is to incorporate the signed distance function into the standard finite element interpolation to represent a discontinuous gradient of the temperature at a moving solid-liquid interface. This technique makes it possible to simulate movement of the solid-liquid interface without the use of a moving mesh. Construction of the finite element equation from the energy equation in the case of melting/solidification problems has been discussed and is reported here. The technique of quadrature and the method to solve the governing equations for the problem involving liquid flows have also been constructed in the present work. The numerical solutions of the basic problems - a one-dimensional Stefan problem, solidification in a two-dimensional square corner, and melting of pure gallium - were compared to the exact solutions or to the experimental data. Through these verifications, validity of the newly developed numerical analysis method has been demonstrated.


I. INTRODUCTION
Numerical analysis of melting0solidification or dissolution0precipitation phenomena is useful to evaluate feasibility of the several candidate techniques in the nuclear fuel cycle. Examples involving such phasechange processes are aqueous reprocessing and fabrication of a solidified glass from high-level radioactive wastes. [1][2][3] There is a difficulty in numerical analysis for the phase-change problems due to the existence of a moving solid-liquid interface. Tracking of the solid-liquid interface is a challenging problem in the numerical analysis.
Numerical analysis methods for melting0solidification problems can be classified into two groups: a moving mesh method and an enthalpy method. In the former group, a mesh is updated so that an element boundary conforms to the solid-liquid interface. Independent energy equations of each phase are coupled by the appropriate boundary condition at the solid-liquid interface. A moving mesh finite element method, which is typical of the former group, was developed by Lynch and O'Neill. 4 Sampath and Zabaras 5 applied the moving mesh finite element method to a three-dimensional directional solidification problem. A variable space network method developed by Murray and Landis 6 and a front-tracking arbitrary Lagrangian Eulerian method developed by Jaeger and Carin 7 are also classified as the former group. While the moving mesh methods give high accuracy in predicting an interface position, there is the drawback that the algorithm of remeshing is complicated in the case of multidimensional problems. When a shape of the solid-liquid interface becomes complex, we have to pay attention to deformation of the elements. Therefore, the moving mesh methods are frequently used for directional melting0solidification problems. The latter group~i.e., enthalpy method! introduces an artificial heat capacity containing the latent heat for phase change. This enables us to eliminate a need for the boundary condition at the solid-liquid interface. The enthalpy method requires only the single energy equation. Numerical analysis using the enthalpy method may be found in the paper by Comini et al. 8 or by Rolph and Bathe. 9 The enthalpy method is quite popular for multidimensional problems because remeshing is not required. However, there is also the drawback that isothermal phase-change phenomena cannot be modeled consistently. This is due to the inevitable assumption that phase change occurs within a certain temperature range. Application of the aforementioned two types of methods is limited by their inherent drawbacks.
A numerical analysis method using a fixed mesh can be simply applied to multidimensional problems even if it involves a complex-shaped moving interface. But calculation of heat conduction and interface tracking in a fixed mesh are complicated by the existence of a discontinuous temperature gradient at the solid-liquid interface. The extended finite element method~X-FEM!, which is a fixed-mesh-based method, introduces an enriched finite element interpolation to represent the discontinuous temperature gradient in the element. The enriched finite element interpolation consists of a standard shape function and a signed distance function. This makes it possible to track the moving solid-liquid interface accurately without remeshing. The X-FEM has the advantages of both the moving mesh method and the enthalpy method. Moës et al. 10 developed the X-FEM for arbitrary crack growth problems. Merle and Dolbow 11 and Chessa et al. 12 applied the X-FEM to melting0solidification problems. Further, Chessa and Belytschko 13 simulated a twophase flow problem involving bubble deformation successfully by using the X-FEM. These researches indicate that the X-FEM is widely applicable to moving boundary problems.
Accordingly, we have employed the X-FEM as a tool to evaluate the melting0solidification or the dissolution0 precipitation processes appearing in the nuclear fuel cycle. But the way to apply the X-FEM to the phase-change problems has not been fully elucidated. In this study, application of the X-FEM to melting0solidification problems has been discussed. Formulation of the enriched finite element interpolation and construction of the finite element equation using it are reported below. The technique of quadrature and the method to solve the problems involving liquid flows have been constructed in the present work. The numerical solutions of the example problems, i.e., a one-dimensional Stefan problem and solidification in a two-dimensional square corner, were compared to the existing solutions to verify validity of our numerical analysis method. Verification of the melting0solidification problem involving natural convection was also conducted.

II.A. Governing Equations
We consider the Stefan problem shown in Fig. 1. The domain V is divided into the solid region V s and the liquid region V l . The boundary of the domain V and the interface between the solid region and the liquid region are denoted by G and G I , respectively. The governing equations for this problem are as follows: continuity equation  Equation~4! means that the net heat flux at the solidliquid interface is translated into the latent heat for phase change per unit time~see Fig. 2!. The interface velocity can be obtained from Eq.~4!. The temperature condition at the solid-liquid interface is given by where T m is the melting temperature. Equation~5! is used as a constraint condition to Eq.~3!, as described later.

II.B. Enriched Finite Element Interpolation
In melting0solidification problems, the temperature gradient is discontinuous at the solid-liquid interface~see Fig. 2!. The basic idea of the X-FEM is to use an enriched finite element interpolation to represent the discontinuous gradient of the temperature. The formulation described in this section can be found also in the paper by Chessa et al. 12 We define the element crossed by the solid-liquid interface as a fully enriched element. The node of the fully enriched element is also defined as an enriched node. Figure 3 shows the definition of the  enriched element and the enriched node in the case of a one-or two-dimensional mesh. The enriched finite element interpolation is based on a standard interpolation: The nodal value of the temperature is defined as Figure 4a shows the temperature profile interpolated by Eq.~6! in the case of a two-node line element. To construct a discontinuous interpolation, we introduce a signed distance function f~x, t !: This definition means the shortest distance from x to the solid-liquid interface. The signed distance function varies as follows: Here, we define an enrichment function c i~x , t !: Figure 4b shows the profile of the enrichment function c 1 and c 2 . As can be seen, the enrichment function is discontinuous at the solid-liquid interface. Figure 4c shows the profile of c 1 ϩ c 2 . This summation becomes linear between the interface and the node. Adding the enrichment function to the standard interpolation function, we obtain the enriched finite element interpolation: where n e is the number of enriched nodes and a j~t ! is the additional degrees of freedom. The temperature profile interpolated by Eq.~11! is shown in Fig. 4d. The twonode line element, which does not contain any enriched nodes~i.e., regular element!, has two degrees of freedom i.e., the unknowns T 1 and T 2 !, while the fully enriched element has four degrees of freedom~i.e., the unknowns T 1 , T 2 , a 1 , and a 2 !. The additional degrees of freedom are decided so that the temperature at the solid-liquid interface becomes the melting temperature. The way to constrain the condition of the temperature at the solid-liquid interface to the finite element equation is described in Sec. II.C.2. The gradient of the temperature is written as

12!
Not all of the nodes are enriched in the element adjacent to the fully enriched element. This partially enriched element is also shown in Fig. 3. If we apply the enriched interpolation defined as Eq.~11! to the partially enriched element, the unnecessary numerical value is added to the desired profile. This problem has been discussed in the paper by Chessa et al. 14 We have used their method to avoid the incorrect interpolation.

5!
The matrix form of Eq.~3! is These coefficient matrices consist of the shape function and the enrichment function. The mass matrix in the left side is written as The other coefficient matrices are given similarly. It should be noted that the size of the global coefficient matrix changes with movement of the solid-liquid interface in the case of a multidimensional problem.

II.C.2. Enforcement of Interface Temperature Condition
The interface temperature condition, i.e., Eq.~5!, acts as a constraint condition to the finite element equation, i.e., Eq.~16!. The constraint condition has been enforced by the penalty method. 12 In the penalty method, a constrained problem is simply converted to an unconstrained one. By using Eqs.~10! and~11!, Eq.~5! is rewritten as Adding the global form of the penalty forcing term to the finite element equation, we obtain the final form: In each of the iterative calculation steps to solve Eq.~23!, the penalty number is updated to be a value larger than that in the previous step. When the penalty number becomes a large value, G QT nϩ1 Ϫ c is enforced to be small.

II.D. Gaussian Quadrature for Enriched Elements
The coefficient matrices in the finite element equation, i.e., Eq.~23!, are calculated by using the Gaussian quadrature. Since the enrichment function is discontinuous~see Fig. 4 incorrect integration. In order to achieve correct integration in the fully enriched element, the element is divided into two subelements as shown in Fig. 5~Ref. 12!.
In the present work, the method to implement the Gaussian quadrature to the subelements has been formulated. Formulation in the case of a two-node line element is described below. The method mentioned here can be simply applied to two-or three-dimensional analyses. For one example, a component of the mass matrix M, i.e., * rc p N i~x ! c j~x ! dx, is divided into two integrations: where V es and V el indicate the region of the subelements.
Here, we introduce a reference coordinate j for the enriched element and a for the subelement~see Fig. 6!. The reference coordinate ranges from Ϫ1 to 1. By using the Gaussian quadrature, the integration in the region V es is given by where P ϭ number of the quadrature points a k ϭ coordinate of the k'th quadrature point W k ϭ weight at the k'th quadrature point.
The terms x~a! and j~a! in Eq.~25! mean that the coordinates x and j are the functions of the coordinate a.
The term J~a! is the Jacobian: The integration in the region V el can be obtained similarly. In Eq.~25!, N i~x ! is replaced by N i~j ! because calculation of N i from the coordinate x is almost impossible in multidimensional analyses using an unstructured  Uchibori and Ohshima ANALYSIS OF MELTING0SOLIDIFICATION PHENOMENA mesh. The mapping from a to j and the mapping from a to x must be given to get the integration. On the coordinate a, Ϫ1 and 1 are translated into j 1 and j I , respectively. Similarly, Ϫ1 and 1 on the coordinate a are translated into x 1 and x I , respectively. From this information, the mappings are given by where N 1 sub~a ! and N I sub~a ! are the shape function shown in Fig. 6. The right sides of Eqs.~27! and~28! are known. We can calculate j~a k ! and J~a k ! from Eqs.~26!,~27!, and~28!, and then obtain the integration.

II.E. Update of Signed Distance Function
The signed distance function is updated by the algorithm of the level set method. 12,15 In this method, the interface motion is governed by the advective equation: where F is the extension velocity, described below. The Galerkin0least-squares method 16 has been applied to Eq.~29! for stabilization. The stabilized finite element equation is where df is the weight function and t e is the stabilization parameter. The formulation of the stabilization parameter has been proposed by Tezduyar. 17 To solve Eq.~30!, we first have to construct the extension velocity, which advects the signed distance function. The extension velocity has the characteristic that it be orthogonal to the signed distance function: sign~f!¹F{¹f ϭ 0 .~31! Figure 7 shows the essential boundary condition for the extension velocity. In Fig. 7, S is the intersection point between the solid-liquid interface and the perpendicular line from the enriched node P. Considering the orthogonal condition given by Eq.~31!, it is obvious that the extension velocity at the node P is equal to the interface velocity at S. Equation~31! is solved by taking the extension velocity at the enriched node as the essential boundary condition. The stabilized form of Eq.~31! is written as

II.F. Calculation of Liquid Flows
In our method, Eqs.~1! and~2! are solved by using the velocity correction method. 18 The Streamline-Upwind0 Petrov-Galerkin scheme 19 has been applied to the convective term in Eq.~2!. The effect of the buoyancy force has been considered by the Boussinesq approximation. To consider the condition that the flow velocity is 0 at the solid-liquid interface~i.e., nonslip condition!, the following essential boundary condition has been applied to the enriched nodes that exist in the solid region: Equation~33! limits the velocity to 0 at the solid-liquid interface. The u ebc approaches 0 as the solid-liquid interface comes close to the enriched node in the solid region.

II.G. Time-Stepping Procedure
The time-stepping procedure is summarized in the following steps: 1. Construct the extension velocity F from Eq.~32! under the essential boundary condition shown in Fig. 7. Equation~4! is used to give the essential boundary condition.
3. Calculate the flow velocity u nϩ1 from Eqs.~1! and~2! under the essential boundary condition given by Eq.~33!. 4. Build the global matrix in Eq.~23! and solve for T nϩ1 by using the penalty method.
5. Return to step 1 and repeat the algorithm.

III.A.1. Problem Description
Numerical analysis of a one-dimensional melting problem using the X-FEM is now presented. The temperature at one boundary is set to a value above the melting temperature, and the temperature at the other boundary is set to the melting temperature. This boundary condition keeps the temperature to the melting temperature in the whole solid region and causes movement of the solidliquid interface. This problem is the well-known onephase Stefan problem. The exact solution describing the interface motion is given by x I~t ! ϭ 2hMa l t ,~34! where x I~t ! is the interface position and a l is the thermal diffusivity of the liquid phase. The constant h satisfies the following relationship: The temperature profile in the liquid region is given by The region, which is 0.02 m in length, was divided into ten finite elements. The initial temperature profile was given by the exact solution at t ϭ 200 s~i.e., the solidliquid interface was initially at x ϭ 2.54 ϫ 10 Ϫ3 m!. The boundary conditions were T lw ϭ 108C at x ϭ 0 m and T sw ϭ 08C at x ϭ 0.02 m. The physical properties of water were used in the analysis. Figure 9 shows the calculated interface positions and the exact solution as a function of time. The numerical analysis without subdivision of the enriched element was also conducted and its result is shown for comparison. The number of the quadrature point P was taken as the parameter. As can be seen, the numerical analysis using the subdivision technique predicted the interface motion successfully. The relative error is within 1%. When the subdivision is not used, the numerical analysis could not reproduce the interface motion even in the case of P ϭ 5. Figure 10 shows the temperature profiles at the four different times. The calculated temperature profiles agree with the exact solutions excellently. From this verification, it has been demonstrated that the X-FEM could predict the interface motion and the temperature profile accurately even in a fixed mesh.

III.B.1. Problem Description
The problem of solidification in a square corner was analyzed as a part of verification of multidimensional problems. An infinitely long prism is initially filled with a fluid at the melting temperature. Figure 11 shows the square-shaped cross section of the prism. The temperature on the peripheral surface of the prism is maintained at a certain value below the melting temperature so that solidification proceeds from the surface inward. There are no changes of the physical properties, and hence there is no convection in the liquid region. As shown in Fig. 11, the verification analysis was performed for the quarter region by taking the centerline of the prism as an adiabatic boundary. The analysis region was divided into the 400 triangle elements. The physical properties of water were used in this analysis. Figure 12 shows the temperature profile and the shape of the solid-liquid interface at the different times. The solid-liquid interface is denoted by the isoline of f ϭ 0. During the initial periods of solidification, the shape of the interface is close to square. As the time progresses, the curve near the diagonal gradually flattens. The interface positions along the centerline, i.e., the adiabatic boundary, and along the diagonal line are shown against time in Figs. 13a and 13b, respectively. The approximate solutions derived by Allen and Severn, 20 Crank and Gupta, 21 Lazaridis, 22 and Rathjen and Jiji 23 are also plotted for comparison. The results of the X-FEM agree with some of the approximate solutions very well. From Fig. 13b, it can be seen that the difference between the numerical results and the solutions by Allen and Severn, Crank and Gupta, or Lazaridis is relatively large. This seems to be due to the approximations used in the previous researches.

III.C.1. Problem Description
Gau and Viskanta 24 conducted the experiment on melting of pure gallium. Their experiment was analyzed to confirm applicability of our numerical analysis method to the problem involving liquid flows. The twodimensional cavity, which is 0.04445 m in height and 0.0889 m in width~see Fig. 14!, is initially filled with pure solid gallium at its melting temperature~29.788C!. At t ϭ 0, the temperature on the left wall is raised to 388C. On the other hand, the temperature on the right wall is kept at the melting temperature. The upper and the lower walls are insulated. As time increases, melting proceeds from the left vertical wall rightward. In this problem, natural convection occurs because of the existence of the gravitational force and the density change. The cavity was divided into the 800 triangle elements. Figure 15 shows the calculated interface, the temperature profile, and the streamlines. It can be seen that the solid-liquid interface inclined rightward from the vertical axis. Since the temperature gradient near the bottom of the liquid region is relatively low, the interface at the bottom moves more slowly. Natural convection plays an important roll in determining the shape of the solidliquid interface. Figure 16 compares the calculated interface shape with the experimental results by Gau and Viskanta. The slight difference is probably due to the uncertainty in the temperature at the right wall in the experiment 24 and the anisotropic thermal conductivity of the pure gallium. 25

IV. CONCLUSIONS
In this study, the X-FEM has been applied to melting0 solidification problems. The method uses the enriched finite element interpolation to represent the discontinuous temperature gradient in the element crossed by the solid-liquid interface. Construction of the enriched finite element interpolation and the finite element equation using  it has been discussed. We have formulated the numerical integration for the enriched element. The method to solve problems involving liquid flows has also been developed in the present work. From the numerical analysis of the one-dimensional Stefan problem, it has been demonstrated that the X-FEM could predict the interface motion and the temperature profile accurately even in a fixed mesh. As a part of verification to multidimensional problems, the solidification in the two-dimensional square corner was also analyzed. The numerical results showed good agreement with some of the approximate solutions.
In the case of problems involving liquid flows, it plays an important roll in determining the shape of the solidliquid interface. In the analysis of the melting of pure gallium, the shape of the solid-liquid interface was reproduced successfully. These verifications indicate that our numerical analysis method is physically correct and is very useful as a tool to evaluate melting0solidification problems. The future work is to develop a mass transfer model to simulate multicomponent dissolution0precipitation phenomena.