Abstract
We present two simple numerical methods to find the free boundary in one-phase Stefan problem. The work is motivated by the necessity for better understanding of the solution surface (temperatures) near the free boundary. We formulate a log-transform function with the unfixed and fixed free boundary that has Lipschitz character near free boundary. We solve the quadratic equation in order to locate the free boundary in a time-recursive way. We also present several numerical results which illustrate a comparison to other methods.
1. Introduction
The free boundary in a Stefan problem can give prediction of the boundary between two opposing materials. For a famous example, the boundary of an ice in water can be found according to time. In this case, water above zero degrees Celsius and ice at zero degrees Celsius are the two materials. Stefan [1] first gave conditions necessary to find the free boundary. Stefan problems can be used for the computation of diffusion [2], preservation of human blood [3], solidification of metals [4], and estimation of the thermal state in the mantle [5]. Many researchers have studied to figure out the effective method in order to approximate the free boundary.
Methods can be largely separated into three groups: analytic methods, mixed methods of numerical and analytic methods, and another method which is numerical methods. Analytic methods are developed by researchers, such as exact solution by using a fractional differential equation of Liu and Xu [6] and an analytical solution with variable latent heat of Voller et al. [7]. Mixed methods are developed by researchers, such as the approximate analytical and a numerical solution method by using a linear perturbation method of Yigit [8] and the analytical and numerical solution method with moving boundary of Lorenzo-Trueba and Voller [9]. Numerical methods are developed by researchers, such as the numerical method with automatic differentiation of Asaithambi [10] and the finite element approximation method of Nedoma [11], Słota [12], and Hinze and Ziegenbalg [5].
Among many methods, a boundary immobilization method (BIM), suggested by Landau [13], fixes the free boundary by using a transform method. The benefit of this function is its simple domain. Whatever the shape of the free boundary is, this method changes the original domain to a rectangular domain. After the BIM method, Kutluay et al. [14] suggested another method called variable space grid (VSG) method.
The main contribution of this paper is the development of two simple numerical methods to find the free boundary in a time-recursive way. Our result is motivated by the necessity for better understanding of the solution surface near free boundary. We adopt the boundary immobilization [13] to change the unknown free boundary to a known and fixed boundary. On the other hand, we adopt the constant number of space intervals between two boundaries [14].
We exploit a log-transform function with the unfixed and fixed free boundary that has Lipschitz character which reduces the accumulating error of the solution surface near the free boundary.
We apply extrapolation near the free boundary. Thus, we can determine the free boundary by solving a quadratic equation in a time-recursive way. Our method also provides fast and accurate results for calculating the free boundary and determining solutions.
The structure of the paper is as follows. Section 2 presents the model formulation. The log-transform function with the unfixed and fixed free boundary to calculate the free boundary is presented in Section 3. Numerical results and comparative studies are presented in Section 4. Section 5 summarizes the paper.
2. Problem Formulation
In this section, we present a mathematical formula for determining solution of one-phase Stefan problem.
The temperature of liquid is denoted by , where is the liquid/solid space for and is the time to expiration for .
As seen in the previous article by Stefan [1], the temperature of liquid is considered the solution to a free boundary problem with a parabolic PDE. We suppose that the free boundary is continuously increasing with . The region where it is free boundary to hold, generally called the liquid region, is defined as , and the region where it is free boundary, generally called the solid region, is defined as . Then, and uniquely solve being subject to the boundary conditions The location of the liquid/solid interface is given by the heat balance equation known as Stefan condition Initially there is liquid region which implies the condition Referring to [14], this problem has the following exact solution for the liquid temperature distribution and the interface location , respectively: where and plug into (2). The exact solution (7) can be used to compare the numerical solutions and can also be used to initialize the numerical schemes.
A VSG method, proposed in Kutluay et al. [14], which uses a fixed node number between a fixed boundary () and a moving boundary (), was kept constant and equal to so that the moving boundary always lays on the th grid. They derive the equation and the boundary conditions with respect to as follows. First, the process finds the difference between the and and divides it by , the number of nodes. In the case below, they assign as the number of and as the number of . Therefore, so that where is the maximum of . Similarly, so that . We also define , , , and . Consider The ratio equals the ratio . Thus, Since , When , By discretizing (10), they get Plugging (12) into (13), they obtain the following: In order to get , they use
A boundary immobilization (BIM) method, proposed in Landau [13], uses a change in variables to transform the free boundary problem into a nonlinear problem on a fixed domain. The following transformation of state variable serves for such a purpose. The main point of BIM method is transformation. The transform function changes any shape of domain to a rectangular domain. when , and when . This fixes the free boundary to . In the case below, he assigns as the number of and as the number of . Therefore, so that where is the maximum of . Similarly, so that . He also define , , and Thus, When , By using central scheme, he gets Discretizing, he obtains Plugging (19) into (20) when , he gets He uses the following equation in order to trace :
3. Log-Transform Function
In this section, we present a log-transform function with the unfixed free boundary and fixed free boundary that can determine the free boundary by solving a quadratic equation in a time-recursive way. Under the assumption of the Stefan problem, the time for free boundary can be shown to be the first hitting time of a boundary, the free boundary, in the plane consisting of pairs of the space and the time to expiration. Namely, the temperature curve of Stefan problem touches the line representing the intrinsic value tangentially. With a careful examination of the solution surface near the free boundary, we find a Lipschitz surface which reduces the accumulating error of the solution surface near free boundary. To find the free boundary, we present a log-transform function with the unfixed and fixed free boundary that has Lipschitz character near the free boundary as follows: The transformed function provides that the solution surface in is a horizontal plane, and it is an inclined plain in . Namely, this function forms a gradual slope with the hyperplane corresponding to the solid region, thereby making the interpolation error more easily reducible when finding the free boundary (see Figure 1). also has a Lipschitz character with nonsingularity in . More precisely, we have . Then, is Lipschitz continuous and a natural candidate function for computation in . We also calculate the partial derivative with respect to in (10), (17).
(a)
(b)
Figure 1(a) shows that is transformed to and Figure 1(b) shows reduced interpolation error near the free boundary.
We find the log-transform function with the unfixed and fixed free boundary to decide the free boundary by the Taylor series.
3.1. Log-Transform Function with Unfixed Free Boundary (LTUF)
From (1), (3), and (4), we obtain the following relations near free boundary: From , we obtain the following relations near free boundary : Plugging (25) into (10), we obtain Using , (24), and (25), we obtain the following: We lower the error by using the optimized safety parameter for a Taylor series with moving node, a method suggested by Kim et al. [15], in order to control the error: Using , (27), and (29), we get Combining (30) with , , we have where For discretization , we introduce a two-dimensional mesh in the first quadrant of the plane. From (31) we have where When the initial values are given by (transformed value of the temperature) and (free boundary) at , we can determine (free boundary) at using (33) and find from (26). More importantly, for updating the free boundary our method does not include any iteration until the numerical solution is obtained. So, we repeat the previously mentioned process until and obtain the free boundary in a time-recursive way.
3.2. Log-Transform Function with Fixed Free Boundary (LTFF)
In the case below, we assign as the number of and as the number of . Therefore, so that , where is the maximum of . Similarly, so that . We also define , , and . The method starts with We obtain from (17) and (35) From , we obtain the following relations near free boundary (): Plugging (35) and (37) into (17), we obtain Using , (36), and (37), we obtain the following: We lower the error by using the optimized safety parameter for a Taylor series with fixed free boundary, a method suggested by Kim et al. [16], in order to control the error. We use Taylor series so that we can manage the error by employing the safety parameter : Using , (39), and (40), we get Combining (41) with , , we have where For discretization , we introduce a two-dimensional mesh in the first quadrant of the plane. From (42) we have where When the initial values are given by (transformed value of the temperature) and (free boundary) at , we can determine (free boundary) at using (44) and find from (38). More importantly, for updating the free boundary, our method does not include any iteration until the numerical solution is obtained. So, we repeat the previously mentioned process until and obtain the free boundary in a time-recursive way.
4. Numerical Result
In this section, we provide two numerical examples to illustrate our method. We also provide runtimes and computation errors used to compare with the results obtained by other numerical methods such as the VSG method developed by Kutluay et al. [14] and the BIM method developed by Landau [13]. In this section, we compare the four methods based on runtime, L2 norm, and root mean square error (RMSE). All calculations are carried out using a C++ implementation with the a 2.66 Ghz Intel 4 Core CPU with 4 GB RAM. A finite difference method with an explicit finite difference approximation is proposed for our methods. The benchmark results are obtained using the exact solution.
The parameter values used to calculate the free boundary and the values of temperature are , , , and a discrete mesh of nodes.
In Figure 2, we find a numerical optimization . Table 1 also shows the results of the free boundary and the values of temperature with the safety parameters. One can see from Table 1 that the free boundary monotonically increases as the value of increases, but the values of temperature monotonically decrease when increases. They are so gradual that they are not very susceptible to change in the value of .
In Figure 3, Table 2, Figure 4, and Table 3, we take the parameter values used in Figure 2 except for the discrete mesh, plot runtimes, and computational errors compared with the various methods.
Figure 5 and Table 4 show the RMSE and L2 norm error of the LTUF method and compare it with sufficiently inclined slope.
Figure 6 and Table 5 show the RMSE and L2 norm error of the LTFF method and compare it with sufficiently inclined slope. Especially, Figures 5 and 6 and Tables 4 and 5 show the numerical convergence of our methods. Note that the discrete meshes of , , , , and nodes are plotted in Figures 3, 4, 5, 6, 7, 8, and 9, and Tables 2, 3, 4, 5, 6, and 7.
For the sake of clarity, the pairs of (RMSE, L2 norm) of all four methods are graphed in Figure 7. Figure 7 illustrates the pairs in the solution. The figure shows that, in terms of errors (RMSE, L2 norm), LTUF approximates more accurately than VSG does, and LTFF approximates more accurately than BIM does. In Figures 8 and 9, we investigate the functional behaviors of the difference between RMSE (VSG) and RMSE (LTUF), L2 norm (VSG) and L2 norm (LTUF), RMSE (BIM) and RMSE (LTFF), and L2 norm (BIM) and L2 norm (LTFF) with respect to final time change.
Table 6 reports the values of the temperatures for the specific parameter set associated with the discrete mesh of .
In Figure 10 and Table 7, we investigate the free boundary behavior with respect to time change when , and . From (1)~(5) this problem has the following exact solution :
(a)
(b)
5. Final Remarks
We suggested LTUF and LTFF methods. There are approximated solutions that show reasonable agreement.
VSG method suggested by Kutluay et al. [14] shows that accumulating error occurred during calculation by the three-term backward difference near the free boundary, while LTUF method adopts a log-transform function to reduce the accumulating error in solution surface that causes gradual slope. Our method employing a log-transform function with the fixed free boundary solves a nonlinear problem on a fixed domain derived from a free boundary problem. Since the computation process depends on Lipschitz surface, we need to focus on the motion of the solution surface which would be simple to see the minute behavior of solution surface.
Kutluay et al. [14] use three-term backward difference scheme during calculating the free boundary with BIM method suggested by Landau [13]. However, we emphasize that there is no three-term backward difference scheme in LTFF method. Instead, we solve the quadratic equation and directly find the free boundary. In an environment with the rapidly changing free boundary, our straight forward method is a very powerful tool to understand Stefan problem.
Numerical study also shows that overall speed and accuracy comparisons have demonstrated some superiority of our methods over other methods. Further investigations on the extending of the various models could be interesting and useful. Also, more comparative and extensive studies on our method should be helpful for completing the mission of the present work. This remains as our future research work.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.