Radial integral boundary element method for simulating phase change problem with mushy zone

A radial integral boundary element method (BEM) is used to simulate the phase change problem with a mushy zone in this paper. Three phases, including the solid phase, the liquid phase, and the mushy zone, are considered in the phase change problem. First, according to the continuity conditions of temperature and its gradient on the liquid-mushy interface, the mushy zone and the liquid phase in the simulation can be considered as a whole part, namely, the non-solid phase, and the change of latent heat is approximated by heat source which is dependent on temperature. Then, the precise integration BEM is used to obtain the differential equations in the solid phase zone and the non-solid phase zone, respectively. Moreover, an iterative predictor-corrector precise integration method (PIM) is needed to solve the differential equations and obtain the temperature field and the heat flux on the boundary. According to an energy balance equation and the velocity of the interface between the solid phase and the mushy zone, the front-tracking method is used to track the move of the interface. The interface between the liquid phase and the mushy zone is obtained by interpolation of the temperature field. Finally, four numerical examples are provided to assess the performance of the proposed numerical method.


Introduction
For pure media, the phase change takes place at the unique temperature, two phases including the solid phase and the liquid phase and one interface between these two phases * Citation: YAO, H. X., YAO, W. A., ZUO, C., and HU, X. F. Radial integral boundary element method for simulating phase change problem with mushy zone. to transform the arising domain integrals into boundary integrals. In the time domain, the PIM [23] is used to solve nonlinear differential equations, and then the heat fluxes on the moving boundary are obtained. Based on the energy balance equation on the moving boundary, the velocity can be obtained. The front-tracking method is used to track the move of the interface between the solid phase and the mushy zone. However, the interface between the liquid phase and the mushy zone is obtained by interpolation of the temperature field.
2 Governing equations of phase change problem with mushy zone A two-dimensional domain Ω total is considered, where the phase change with the mushy zone takes place. Contrary to the isothermal phase change problems, the interval of phase change temperature is (T ms , T ml ). Therefore, the domain Ω total is divided into three regions, namely solid domain Ω s , mushy zone Ω m , and liquid domain Ω l , by two interfaces Γ ms and Γ ml , as shown in Fig. 1. If volumetric changes are neglected and the density is constant, the heat conduction governing differential equations can be written as where x = (x 1 , x 2 ), ∇ = ∂(·) ∂x1 i + ∂(·) ∂x2 j, T denotes the temperature, t is the time, ρ is the density, and Q denotes the heat source. c represents the heat capacity, and k is the thermal conductivity. The subscripts 's', 'l', and 'm', respectively, denote the solid domain, the liquid domain, and the mushy zone. The heat capacities and thermal conductivities may be non-constant in these three domains, but in Ω total , they are continuous functions with respect to the temperature. The phase change latent heat is in Eq. (2) of the mushy zone L, and f is the solid phase fraction. It is assumed that f equals 0 on Γ ml , f equals f su on Γ ms , and f is a function with respect to temperature or distance to the interfaces Γ ms and Γ ml [25] . The boundary conditions of Ω total are given by where n denotes the unit outside normal vector, and q is the heat flux. The initial condition is given by Besides, there are interface conditions on the interfaces Γ ms and Γ ml . The temperature conditions and the energy equation on the solid-mushy interface Γ ms are given by where V n is the velocity of the moving interface. The temperature conditions and the energy equation on the liquid-mushy interface Γ ml are given by It can be found that from Eqs. (7) and (9), the temperature is continuous. Since the thermal conductivities are continuous in the whole domain, it can also be found from Eqs. (8) and (10) that the temperature gradient is continuous on Γ ml but stepped on Γ ms . However, if the latent heat change in the mushy zone is treated as the source terms, the mushy zone and the liquid phase in the simulation are treated as a whole part, namely the non-solid phase, and let Ω ns = Ω m ∪ Ω l . Then, the differential equation in the non-solid domain is written as follows: where Therefore, the phase change problem can be dealt with by separating the whole domain into two parts, namely the solid domain and the non-solid domain. Not only Eq. (1) of solid domain, Eq. (11) of non-solid domain, and the solid-mushy interface conditions (see Eqs. (7) and (8)), but also the whole domain's boundary conditions (see Eqs. (4) and (5)) and the initial condition (see Eq. (6)) need to be solved. First, the transient heat conduction problems are solved by the precise integration BEM. Then, the moving velocity of the mushy-solid interface Γ ms can be obtained by Eq. (8), and the moving position of Γ ms can be tracked by an iteration procedure of the front-tracking method. Finally, the mushy-liquid interface Γ ml can be obtained by interpolation of temperature field.

Implementation of precise integration BEM
The transient heat conduction problems in the solid domain and the non-solid domain are solved by the precise integration BEM, respectively. The governing equation in the solid domain or non-solid domain can be rewritten in a uniform equation as follows: According to the whole domain's boundary conditions (see Eqs. (4) and (5)) and the solid-mushy interface condition (see Eq. (7)), the boundary conditions for the solid domain or the non-solid domain can be rewritten as follows: The Green function G(x, y) is adopted as the weight function to derive the boundary-domain integral equation. For two-dimensional problems, it is expressed as Both sides of Eq. (13) are multiplied by G. Then, the resulting equation is integrated over the domain Ω. Applying the integration by parts and using Gauss' divergence theorem, one can get where c(y) = 1 when y ∈ Ω, and c(y) = ϕ(y)/(2π) when y ∈ Γ, in which ϕ denotes the interior angle in radians. Equation (17) can be rewritten as [26] c where In Eq. (18), the domain integral, where the source function Q(x, t, T ) is involved, can be directly transformed by the RIM [20] as follows: where In order to transform the third domain integral in Eq. (18), the following approximation needs to be done: where is the radial basis function (RBF), and α i and α i are coefficients.
According to Ref. [27], the transformation is done as follows: where the jth-component of the row vector D y is where the meaning of Ψ ij can be found in Ref. [27], and The same approximation procedure is done for the first domain integral in Eq. (18), The domain integral can be transformed as follows: where the jth-component of the row vector V y is where Substituting Eqs. (23) and (29) into Eq. (18), one can get the integral equation with only boundary integrals as follows: The boundary Γ is discretized into N e boundary elements with N b1 boundary nodes on Γ 1 and N b2 boundary nodes on Γ 2 , which ensures that the total number of boundary nodes is N b = N b1 + N b2 , and N i internal nodes are distributed into the domain Ω. The total number of nodes is N = N b + N i . Then, the temperature vector can be written as After applying Eq. (34) on both boundary and internal nodes, the following differential equations are obtained: where C y = diag(c(y 1 ), c(y 2 ), · · · , c(y N b )), G, H, V , f , and C are respectively generated by discretization of the five boundary integrals on the right-hand side of Eq. (34), and the subscripts 'b' and 'i' denote boundary nodes and internal nodes, respectively.
According to the boundary conditions, for N b1 boundary nodes whose T b1 and˙ T b1 are known, q b1 is to be solved, and for N b2 boundary nodes whose q b2 is known, T b2 and˙ T b2 are to be solved. According to the N b1 equations of Eq. (36), q b1 can be eliminated [27] . Then, the differential equations, where the unknown vector X is involved, can be obtained,

Procedures of predictor-corrector PIM
In order to use the PIM to solve Eq. (38), in the time interval [t k , t k+1 ], Eq. (38) can be rewritten as where A 0 = A(T (t k )), and K 0 = K(T (t k )). Then, one can obtainẊ where ∆t = t k+1 − t k . The calculation of Eq. (41) is carried out by the predictor-corrector iteration technique [28] as follows.
Step 1 Let n iter = 0, and do the following calculation: where X e = EX(t k ), E = exp(B∆t), K = K(T (X e )), and r i0 = r i (t k , T k ).
Step 2 Let n iter := n iter + 1, and assume that r i is linear with respect to t, which is expressed as where 0 < τ ∆t, θ is the relaxing factor, and 0 < θ 1, which is taken as 0.5 in the calculation processes of the numerical examples of this paper.
Do the following calculation: where Step 3 Check the convergence. If X k+1,niter − X k+1,niter−1 2 < ε (where · 2 denotes a vector's 2-norm, and ε is a prespecified acceptable error), the convergence is achieved, X k+1,niter is the solution at t k+1 , and then go to the next time step. Otherwise, go to Step 2.
After X(t k+1 ) is obtained, the temperature field and temperature gradient at t k+1 are all known. Then, the moving interface needs to be tracked.

Numerical steps of front-tracking method
To determine the position of the interface, the moving velocity and the direction are needed. The moving velocity can be calculated by Eq. (8) once the heat flux on the solid-mushy interface is obtained. The moving direction can be determined by a length weighted unite normal vector as the following equation [6,29] : where the subscript j indicates the jth-node on the moving interface, the subscripts i − 1 and i indicate the two adjacent boundary elements of the jth-node on the moving interface, l i−1 and l i denote the lengths of these adjacent boundary elements, and · 2 denotes a vector's 2-norm. The position of the moving interface can be obtained by an iterative algorithm. The front-tracking procedures in one time step [t k , t k+1 ] are described in detail by the following steps.
Step 1 Specify a tiny distance ∆s beforehand for the purpose of controlling the maximum moving distance of Γ ms .
Step 2 At the initial time t k , the heat flux q(t k ) at all the boundaries is obtained from the results of the last time step, boundary conditions, or initial conditions. Then, the normal velocity V j (t k ) of the jth-node on the moving interface Γ ms at t k can be calculated by Eq. (8), and V (t k ) = {V j (t k )}.
Step 3 The time step ∆t k+1 is determined by the following equation: where · ∞ denotes a vector's infinity norm, and t k+1 = t k + ∆t k+1 .
Step 4 Estimate V p (t k+1 ) = 0.7V (t k ), where V p (t k+1 ) denotes the predicted velocity of the moving interface Γ ms at t k+1 .
Step 5 Calculate the predicted position of the jth-node on Γ ms at t k+1 by and update the geometry at t k+1 .
Step 6 Apply the radial integration BEM in the solid domain and non-solid domain, and then use the procedures described in the last section to calculate the temperature at t k+1 by the PIM.
Step 7 Calculate the unknown heat flux at t k+1 in the solid domain and the non-solid domain, respectively, and then calculate the normal velocity V j (t k+1 ) of the jth-node on the moving interface Γ ms at t k+1 by Eq. (10), and V (t k+1 ) = {V j (t k+1 )}.
Step 8 Check the convergence.
where ε is a prespecified acceptable error, the convergence is achieved. Otherwise, let V p (t k+1 ) = V (t k+1 ), and then go to Step 5.
Step 9 The position of Γ ms obtained in Step 5 and the temperatures calculated in Step 6 are assumed to be the results at t k+1 . Obtain the mushy-liquid interface by interpolation of temperature field. Do the re-meshing. Add or delete the boundary elements, boundary nodes, and internal nodes in both solid and non-solid domains if needed. Let t k := t k+1 , and go to Step 1. Execute the iterative algorithm for the next time step until the ending time is reached.

Numerical examples
In the numerical processes, the heat capacities and thermal conductivities are assumed to be constant in the solid domain and the liquid domain, but linear with respect to temperature in the mushy zone [30] , which are given by where f is assumed to be linear with respect to temperature as follows: Four numerical examples for the phase change problems with the mushy zone are presented to assess the performance of the proposed numerical algorithm. Analytical solutions are used as reference for Example 1, while the solutions by a combination of fixed-domain method and FEM are used as reference for Example 2, Example 3, and Example 4.
Example 1 The freezing problem in a cylindrical symmetry region with an extended freezing temperature range is studied. The considered region is liquid initially with the constant temperature T 0 (T 0 > T ml ). When t > 0 s, a heat sink of strength Q is located at r = 0. According to Ref. [25], the analytical solutions of the solid-mushy interface and liquid-mushy interface are given by The analytical solutions of the temperature field are given by where E i (x) is the first-order exponential integral function, α s = k s /(ρc s ), α m = k m /(ρc m ), and α * ms = (α m k m (T ml − T ms ))/(ρLf su α m + k m (T ml − T ms )). This problem can be approximated by an annulus region 0.1 < r < 1 as shown in Fig. 2. Consider the phase change material-1 (PCM-1) as shown in Table 1. The initial condition is specified by The boundary conditions are specified by The thermal-physical properties, the initial condition, and the boundary conditions determine λ = 0.343 331, and η = 0.758 593.
The relative errors of the interface positions of the present method and the analytical solutions are shown in Figs. 3 and 4 when ∆s is different, and the total number of nodes is 320. The maximum relative error of the solid-mushy interface position is 4.29%, while the maximum relative error of the liquid-mushy interface position is 2.79%. The relative errors of the solid-mushy interface positions of the present method and the analytical solutions are shown in Fig. 5, when the total number of nodes is different, and ∆s = 4 × 10 −3 . The maximum relative error of the solid-mushy interface position is 3.45%. If the total number of nodes is 320, the maximum relative error of the solid-mushy interface position is 0.99%. It can be found that the results of the present method agree well with the analytical solutions, and both the decrease in the step ∆s and the increase in the total number of BEM nodes would increase the accuracy. The present method is coded by MATLAB and computed by Intel Core i5-4460 CPU, and the CPU time is listed in Table 2. Example 2 A solidification problem of a 1 × 1 domain is studied as shown in Fig. 6. The considered domain is full of PCM with the constant initial temperature T 0 = 3 • C. When t > 0 s, a temperature boundary condition T w = 0 • C is applied at x 1 = 0 m and x 2 = 0 m. Consider the properties of the PCM-2 in Table 1. The locations of the interface at various time are shown in Fig. 7, where 'S-BEM' and 'L-BEM' represent the BEM for solid-mushy interfaces and liquid-mushy interfaces, respectively.
It can be found that there is good agreement between the interface position results of the present method and the FEM solutions.  Example 3 The solidification problem of a phase change heat storage system as shown in Fig. 8 is studied. The considered region is full of PCM with the constant initial temperature T 0 = 3 • C. When t > 0 s, a temperature boundary condition T w = 0 • C is applied at both the interior and exterior walls. Consider the properties of the PCM-3 in Table 1.
The locations of the interface at various time are shown in Figs. 9-12. It can be found that there is good agreement between the interface position results of the present method and the FEM solutions. The topology change of the solid domain can be seen from Fig. 10. From Figs. 11 and 12, the topology changes of the liquid domain can also be found. The proposed method can solve phase change problems with topology changes for both the solid domain and the liquid domain accurately.
Example 4 The solidification problem of a phase change heat storage system as shown in Fig. 13 is studied. The area between the interior and exterior copper walls with fins is occupied by the PCM, where the phase change takes place. The initial temperature of the considered domain is T 0 = 40 • C. When t > 0 s, a temperature boundary condition T w = 15 • C is applied at both the interior and exterior walls. Consider the properties of the PCM-4 in Table 1. In the copper domain, the thermal conductivity is 400 W/(m · K), the specific heat is 380 J/(kg · K), and the density is 8 920 kg/m 3 . There is not any phase change in the copper domain, but the heat conduction in this domain needs to be considered and numerically analyzed. The BEM is applied in the copper domain and the adjoining PCM domain, respectively, to get two sets of differential equations. At the interface between the copper domain and the PCM domain, the temperature is continuous, and the heat flux is equilibrated. Therefore, according to this interface condition, these two sets of equations can be assembled into one. Then, the PIM results of temperature and heat flux can be obtained for both the copper domain and PCM domain.
Because of the symmetry, only 1/8 part of the storage device is considered in the numerical calculation. The locations of the interface at various time are shown in Fig. 14. It can be found that there is good agreement between the interface position results of the present method and the FEM solutions. The topology change of the solid domain can be seen from Figs. 14(b) and 14(c). The topology change of the liquid domain can be seen from Fig. 14(b). The proposed method can solve phase change problems with topology changes and multi-medium heat conduction accurately.

Conclusions
This article proposes a precise integration BEM to solve the phase change problem with the mushy zone. The changes of latent heat are approximated by heat source, and the liquid phase and the mushy zone are considered as a whole part. The front-tracking method is used to track the move of the interface between the solid phase and the mushy zone, and the interface between the liquid phase and the mushy zone is obtained by interpolation of the temperature field.
Compared with the reference solutions, the results of the numerical examples show that the proposed method can determine the phase change interface positions correctly. It can be found that both the step ∆s and the number of BEM nodes would influence the computing accuracy dramatically. The proposed scheme is very convenient to be applied in the moving-mesh of the grid, and can properly deal with topology changes including separation merging and vanishment. It can also deal with multi-medium problems by carrying out the BEM in different domains. It is reliable and accurate to apply the proposed numerical method for solving phase change problems with mushy zones.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International
License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.