Numerical Solution for Two-Sided Stefan Problem

In this paper, we consider a two-phase Stefan problem in one-dimensional space for parabolic heat equation with non-homogenous Dirichlet boundary condition. This problem contains a free boundary depending on time. Therefore, the shape of the problem is changing with time. To overcome this issue, we use a simple transformation to convert the free-boundary problem to a fixed-boundary problem. However, this transformation yields a complex and nonlinear parabolic equation. The resulting equation is solved by the finite difference method with CrankNicolson scheme which is unconditionally stable and second-order of accuracy in space and time. The numerical results show an excellent accuracy and stable solutions for two test examples.


1.Introduction
The classical Stefan problem is the name given to an initial -boundary value problem, which involves both fixed and moving boundaries. Stefan problems model has been involved in many realworld engineering situations in which there is a freezing or melting that are causing a boundary to change in time, as in melts solidification, water and food freezing, crystal growth, etc. Stefan direct problems are boundary value problems for parabolic heat equations in regions with unknown and moving boundaries, which requires determining the temperature [1]. In a previous article [2], the numerical solution for the free boundary problems has been investigated through the finite difference method (FDM) and the minimax approach. Whilst, the numerical solution for the free problem in fluid mechanics such as gas-liquid interface problem was considered in another ISSN: 0067-2904

Adel and Hussein
Iraqi Journal of Science, 2020, Vol. 61, No. 2, pp: 444-452 445 study [3]. Also, a one-side Stefan problem for parabolic equation was solved numerically using shifted Chebyshev operational matrix [4]. Other authors [5,6] were concerned with the numerical solution of one-phase Stefan problem in addition to the reconstruction of some unknowns which hold physical meaning such as the thermal conductivity, heat capacity, and fluid velocity. It is worth to mention that Stefan problems are basically restricted to heat transfer problems. However, they can be applied to full parabolic equations.
In this work, we consider the numerical solution for a two -sided Stefan problem where the free boundary depends on time only. This problem is solved by FDM with Crank-Nicolson scheme. This paper is organized as follows. In Section 2 the mathematical formulation of the two-sided Stefan problem is given. The numerical solution for the problem under consideration using FDM is discussed in Section 3. Whilst, the stability analysis for the proposed method is considered in Section 4. In Section 5, the results and discussion are presented. Finally, conclusions are highlighted in Section 6.

2.Mathematical formulation
Consider the one-dimensional two-phase Stefan problem [7], ( ) and the non-homogenous Dirichlet boundary conditions are are given and u is the solution of problem, i.e. temperature. In order to solve this problem, we change the variables  [7] and no numerical solution was obtained. After performing the transformation, the complicated problem turns to a fixed domain problem. However a nonlinear equation is obtained.
The unique solvability of the direct problem is guaranteed by the continuity of the coefficients , and , as previously reported [8]. Now, the problem in fixed domain has the following form which will be solved numerically in the next section using a finite-difference scheme.

Solution of the direct problem
In this section, consider the direct initial boundary value problem (IBVP) ( ) ( ) where the functions ( ) ( ) ( ) ( ) and ( ) are known and the solution ( ) is to be computed. We employ the Crank-Nicolson finite-difference scheme which is unconditionally stable and second order accurate in time and space [9].

Discretization of the problem
The discrete form of the problem ( ) ( ) is as follows: We divide the domain ( ) ( ) into M and N subintervals of equal lengths and , where the uniform space and time increment are , , respectively. We denote the solution at the node point ( ) In order to apply the CN-scheme for equation (4), we simply assume the right hand side as follows: After simple arrangement we obtain the following difference equation

Stability Analysis
It is worth to mention that the numerical scheme still has the stability properties, which are  [9] for more details. This proves that the Crank-Nicolson scheme is unconditionally stable.

5.Numerical Results and discussion 5.1 Example 1
In this example, we consider the case when the coefficients in equation ( ) are set of polynomials of first order in and . Moreover, we assume the free boundary as a linear function in time, as illustrated in the following quantities.
We test our numerical scheme for multivalues of and which starts as { } For simplicity, we take . Figures and Table show the numerical solution for the direct problem ( ) ( ) with various mesh sizes. In comparison with the exact solution, the absolute error between exact and numerical solutions is included. From these figures, it can be easily seen that an excellent agreement is obtained for all the selected mesh sizes. A very low magnitude of error is obtained, which is of the order ( ) ( ) This result is expected since the scheme is unconditionally stable and a second order accurate for space and time.

Example 2
In this example we consider a nonlinear case for the coefficients, and a linear for the free boundary .  Table 2 show a very good agreement with exact solution has been obtained. However, the magnitude of error is ( ) ( ) which is higher than that in the previous example, and this is expected since nonlinear coefficients are used. However, the results are still acceptable and free from occilations and unstable behavior.

Conclusions
The problem of free (Stefan) boundary in the full parabolic heat equation with non-homogenous Dirichlet boundary conditions has been numerically investigated. The transformed problem which is in a fixed domain has been solved using FDM with Crank-Nicolson scheme which is unconditionally stable and second order accurate in space and time. The proposed method is efficient and produces a highly accurate solution. Two test examples are introduced in order to test and explain the accuracy and stability of the used scheme. The results are accurate and satisfactory.