A numerical method for the solution of the two-phase fractional Lame-Clapeyron-Stefan problem

In this paper we present a numerical solution of a two-phase fractional Stefan problem with time derivative described in the Caputo sense. In the proposed algorithm, we use a special case of front-fixing method supplemented by the iterative procedure, which allows us to determine the position of the moving boundary. In the final part, we also present some examples illustrating the comparison of the analytical solution with the results received by the proposed numerical method.


Introduction
Moving boundary problems ares one of the most important area within partial differential equations. This particular kind of boundary-value problem was originally intended to describe the solid-liquid phase change process, but also refers to such phenomena as solute transport, molecular diffusion or controlled drug release [1,2,3]. Moving boundary problems, as the name implies, are characterized by having a moving boundary of the domain, which has to be determined as a part of the solution. These kinds of problems are also called Stefan problems, in connection with the early work of Slovene scientist Joseph Stefan, who investigated ice formation in the polar Arctic seas [4]. In recent years, the classical Stefan problem has been very well studied and described in many papers and monographs (compare [5,6,7,8,9] and the references therein).
The fractional Stefan problem is a natural generalization of th classical Stefan problem. The first paper devoted to this issue was published in 2004 and concerned the mathematical modeling of the controlled release of a drug from slab matrices [10]. Recently, there has been an increase in the number of scientific publications concerning the moving boundary problems modeled by the anomalous diffusion equation. Basically, published papers deal with three classes of phenomena. As mentioned earlier, the first concerned the controlled release of a drug from slab matrices [11,12]. A second class of problems relates to mathematical modeling of the thermal conductivity with phase transitions [13,14,15,16] and the last class refers to mathematical modeling a movement of the shoreline in a sedimentary ocean basin [17,18].
Most of the analytical solutions of fractional Stefan problem were obtained for the one-phase case. Liu Junyi and Xu Mingyu [10] studied the one-phase Stefan problem with fractional anomalous diffusion (Riemann-Liouville derivative with respect to time variable was used) and got an exact solution (concentration of the drug in the matrix) in terms of the Wright's function. They also showed that position of the penetration of solvent at time t moving as ∼ t α 2 , 0 < α ≤ 1. Xicheng Li et al. [19] considered one-phase fractional Stefan problem with Caputo derivative with respect to time and two types of space-fractional derivatives (Caputo and Riemann-Liouville). They got the solution in terms of the generalized Wright's function. It should be noted that in both cases deliberated by the authors, function describing the moving boundary is ∼ t α β , 0 < α ≤ 1, 1 < β ≤ 2, where α, β denotes the orders of the fractional derivatives with respect to time and spatial variable respectively. In the paper [14] Voller analytically solved a limit case fractional Stefan problem describing the melting process. For the governing equation with a Caputo derivative with respect to time of order 0 < β ≤ 1 and for the same fractional derivative with respect to space for the flux of order 0 < α ≤ 1, he showed that a melting front is described by power function t β α+1 . More analytical results were published in the papers [20,21,22].
An important result with respect to this paper is the closed analytical solution of the two-phase fractional Lamé-Clapeyron-Stefan problem in a semi-infinite region obtained by Roscani and Tarzia [23], which will be recalled in Section 3. The problem involves determination of three functions, namely, u 1 , u 2 and S fulfilling two subdiffusion equations (13), (14) and additional differential equation (interface energy balance condition) (18) governing function S. The authors showed that u 1 and u 2 can be expressed in terms of the Wright's function, moreover location of the phase-change interface is described by power function (21), where parameter p can be evaluated form transcendental equation (22). This solution will be especially useful for validation of the numerical method proposed in Section 4.
An alternative to the closed analytical solutions are those received by numerical methods. There are many techniques of numerical solving of classical Stefan problem, some of them have been generalized to the case of the fractional order. Xiaolong Gao et al. [24] generalized the boundary immobilisation technique (also known as frontfixing method [5]) to the fractional Stefan problem with a space-fractional derivative. Another variant of the front-fixing method was proposed in paper [25] and applied to the fractional Stefan problem with time-fractional derivative.
Our aim is to develop the numerical method of solving the two-phase, one-dimensional fractional Stefan problem. The proposed numerical scheme is an extension of the frontfixing method [25] to the two-phase problem. Our new approach is based on the suitable selection of the new space coordinates. The original coordinate system (x, τ ) is transformed into a two new orthogonal systems (v 1 , τ ) and (v 2 , τ ) using transformation (27) and (43), respectively. Both new spatial variables v 1 , v 2 depend on the parameter p which is unknown and chosen a priori. The proposed numerical method uses integrodifferential equations equivalent to the corresponding governing differential equations of the problem. The solutions of the integro-differential equations are obtained separately for each phase and fulfill the interface energy balance condition (the moving boundary is fixed ) only when the value of parameter p is correct. Selection of the appropriate value of parameter p is implemented by iterative algorithm on the basis of the fractional Stefan condition.
The paper is organized as follows. In the next Section, we introduce definitions of the fractional integrals and derivatives together with some of their properties. In Section 3, we formulate a mathematical model describing the melting process and recall closed analytical similarity solution for two-phase fractional Stefan problem. Section 4 is devoted to the new numerical method of solving of two-phase fractional Stefan problem, which is a extension of the method developed in the paper [25] to the twophase case. Section 5 contains the analytical and numerical results. The last Section includes a summary of the paper and conclusions.

Preliminaries
Let us recall the basic definitions and properties from fractional calculus [26,27] which will be further applied to formulate and solve the two-phase fractional Stefan problem modeled by two linear equations with Caputo time derivative of order α ∈ (0, 1]. First, we define the left-sided Riemann-Liouville integral.

Definition 1
The left-sided Riemann-Liouville integral of order α, denoted as I α 0+ , is given by the following formula for Re(α) > 0: where Γ is the Euler gamma function.
The left-sided Caputo derivative of order α ∈ (0, 1) denoted by c D α 0+ is defined via the above left-sided Riemann-Liouville integral c D α 0+ f (t) := I 1−α 0+ f ′ (t). Definition 2 Let Re(α) ∈ (0, 1]. The left-sided Caputo derivative of order α is given by the formula: An important property of fractional integral operators given in Definition 1 is called a semigroup property, which allows composition of two left-sided fractional integrals. Property 3 (cf. Lemma 2.3 [26]) If Re(α) > 0, and Re(β) > 0, then the equation The composition rule of the left-sided Riemann-Liouville integral with the left-sided Caputo derivative is a consequence of the semigroup property for fractional integrals.
Property 4 (cf. Lemma 2.22, [26]) Let function f ∈ C 1 (0, b). Then, the composition rule for the left-sided Riemann-Liouville integral and the left-sided Caputo derivative is given as follows: The two-parameter Wright function determined in a complex plane plays important role in the theory of partial differential equations of fractional order. It is a generalization of the complementary error function.
The two-parameter Wright function is given as the following series: We note that for γ = − 1 2 , δ = 1 the above two-parameter Wright function becomes complementary error function [28]: For γ = − 1 2 , δ = − 1 2 two-parameter Wright function can be expressed by the formula [29]: The two-phase fractional Stefan problem is a mathematical model describing the solidification and melting process. The mathematical formulation of the problem involves an anomalus diffusion equation for the liquid and solid phases and a condition at the liquid-solid interface, called the fractional Stefan condition, which describes the position of the phase change front. At the moving boundary X = s(t), the temperature U is constant and equal to melting point U s . We are considering the melting of a semi-infinite, one dimensional slab occupying X ≥ 0, where the phase change interface moves from the heat source at a temperature U = U 0 at the boundary X = 0. A simple scheme of the model is shown in Figure 1.

Liquid Solid
Interface U(x,t) X Figure 1: Semi-infinite slab melting from x = 0 due to the high temperature U 0 .
Consider the following governing equations of the model where U 1 denotes the temperature in the liquid phase, U 2 the temperature in the solid phase, s(t) represents the position of the moving boundary, K the modified thermal conductivity (has SI unit [J · s −α · m −1 · K −1 ]), ρ the density, c the specific heat and subscripts 1 and 2 indicate liquid and solid phase, respectively. Equations (8,9) should be supplemented by Dirichlet conditions: At t = 0 the liquid phase does not exist, so we use two initial conditions: The position of the moving boundary s(t) is determined by the fractional Stefan condition, which expresses the heat balance in the melting layer: where L is the latent heat. Let us note that the above mathematical model depends on eight parameters. To simplify of the studied problem, we introduce the following dimensionless variables which reduces the number of free parameters to five and makes it possible study their impact on the solution. We can rewrite governing equations (8,9) in a non-dimensional form supplemented with the boundary conditions and the fractional Stefan condition where , are standard values of the respective variables.
According to results obtained by Roscani and Tarzia [23] the closed analytical solution of the two-phase fractional Stefan problem (13)-(18) is given by the functions It should be noted that the above solution reduces to the fractional one-phase problem [20,21] for u 2 → 0 . When α → 1, then from (19)(20)(21)(22) the following results of the classical two-phase Stefan problem are recovered:

Numerical solution
The numerical method proposed in the paper is an extension of the technique developed in [25] to the two-phase problem. Applying the finite difference method to the fractional Stefan problem, we encounter a number of difficulties related to the discretization of the fractional derivative with respect to the time variable. As we know, the Caputo derivative is a non-local operator, which is defined on an interval. Suppose the positions, that the moving boundary will reach at different times are known and are at uniform spaced intervals in the same time layer, see at the left side of Figure 2.
For such a grid, the discretization of the Caputo derivative at some node with respect to time is very difficult because it requires values of function u 1 (or u 2 ) for all previous times (where the spatial variable is fixed) in points, which do not overlap with the mesh nodes. Let us also notice one important fact regarding discretization of equation (13) leading to some ambiguity. The Caputo derivative requires for each point of the liquid phase their history from time τ = 0, but they do not exist before they reach the melting point. It seems that the only way to solve a mesh problem is a suitable choice of new space coordinates for equations (13) and (14). Let us first consider the region occupied by the first phase bounded by x = 0 and x = S(τ ), marked in blue in Figure 2. For the equation (13), it is convenient to replace the spatial variable with the following similarity variable [30]: which has an important property, namely fixes the moving boundary at v 1 = 1 for all τ . We transform the first order derivative with respect to the spatial variable: Subsequently, for the second order spatial derivative we have: Respectively, for the first order partial derivative in time: Caputo derivative of function u 1 (v 1 , τ ) with respect to the time variable can be expressed as follows: Using formulas (29) and (31) we write equation (13) in coordinate system (v 1 , τ ): We integrate the previous equation applying the left-sided Riemann-Liouville integral of order α ∈ (0, 1) and we get the following equation: Finally, using Properties 3 and 4 to equation (33), we obtain an integro-differential equation in the form of: The kernel of the second integral on the right-hand side of formula (34) causes some difficulties in deriving a numerical scheme. For this reason we are introducing an auxiliary function:ū which leads to the integro-differential equation: supplemented with the boundary conditions and initial conditionū 1 (v 1 , 0) = 0. The region marked in blue in Figure 2 in the (x, τ ) plane is transformed into the unit square (in the general case to the rectangle) using the transformation (27). Rectangular mesh created in this way consists of horizontal lines spaced ∆τ = 1/(np 2/α ) units apart and vertical lines spaced ∆v 1 = 1/m 1 units apart. The points ((v 1 ) i , τ j ) = (i∆v 1 , j∆τ ) where i = 0, ..., m 1 and j = 0, ..., n, of intersection of the horizontal and vertical lines are grid nodes. We denote the value of functionū 1 at point (i∆v 1 , j∆τ ) by (ū 1 ) i,j . At this stage, the value of parameter p is not known and chosen a priori. Construction of similarity variable (27) requires an additional assumption for variable τ . Let τ 0 be a very small positive number.
The method of discretization of the integro-differential equation (36) is described in detail in the paper [25] and leads to an implicit scheme: where Weights (39) are the result of application of the trapezoidal rule [31] to the last integral term in formula (36). Obtained numerical scheme can be written in the matrix form whereŪ 1 is a vector of unknown values of functionū 1 at instant τ k+1 .
The elements of matrices A and B are defined as follows: We recover the values of function u 1 in coordinate system (x, τ ) by applying formulas The region marked in red in Figure 2 in the (x, τ ) plane is transformed into the unit square using the following transformation: which fixes the moving boundary at v 2 = 0 for all τ . Let us note, that the mesh for the semi-infinite region contains an infinite number of nodes, which makes it impossible to perform any numerical calculations. Therefore, from a practical point of view the Dirichlet boundary condition in infinity (15) can be replaced by the following condition: where l 2 is a sufficiently large positive real number. We transform the first order derivative with respect to the spatial variable: Consequently, for the second-order spatial derivative Respectively, for the first order partial derivative with respect to the time variable: Caputo derivative of function u 1 (v 1 , τ ) with respect to the time variable can be expressed as follows: (48) Using formulas (46) and (48) we can write equation (14) in the new coordinate system (v 2 , τ ): We integrate both sites of equation (49) applying the left-sided Riemann-Liouville integral of order α ∈ (0, 1): Finally, using Properties 3 and 4 to equation (50), we obtain an integro-differential equation in the form of: We introduce the second auxiliary function which leads to the integro-differential equation: supplemented with the boundary conditions and initial conditionū The integro-differential equation (53) is discretized in analogy to equation (36). For the second phase we also operate on a rectangular uniform grid consist of horizontal lines spaced ∆τ = 1/(np 2/α ) units apart and vertical lines spaced ∆v 2 = 1/m 2 units apart. The points ((v 2 ) i , τ j ) = (i∆v 2 , j∆τ ) where i = 0, ..., m 2 and j = 0, ..., n, of intersection of the horizontal and vertical lines are grid nodes. At this stage of the calculation, the value (for both phases the same) of parameter p is not known and chosen a priori.
The implicit numerical scheme for second phase is in the form of: where The obtained numerical scheme can be written in the matrix form whereŪ 2 is a vector of unknown values of functionū 2 at instant τ k+1 .
The elements of matrices D and E are defined as follows: Applying the following formulas: we recover the values of function u 2 in coordinate system (x, τ ). The presented numerical scheme allows us solve the governing equations (13,14), but only when the correct value of the parameter p is known. We can determine it using the fractional Stefan condition. First, we integrate both sides of equation (18) applying the left-sided Riemann-Liouville integral of order α ∈ (0, 1). Next, using Property 4 we obtain: Applying the difference quotient approximation for the first order derivative with respect to space variable, weights (39) and condition (17), we get: (62) Let us note that the value of function S at the final time instant τ n for the correct value of parameter p should be equal to 1. From this relation the following criterion of convergence results in: where ǫ > 0 is a certain arbitrarily small real number. We denote as S(τ n , p) the value of S(τ n ) for a fixed value of parameter p. Below we present an iterative algorithm for determining the value of parameter p, based on the bisection method [25]. . if |1 − S(τ n , p c )| < ǫ is fulfilled, then end the algorithm, otherwise go to the next step 8. if (1 − S(τ n , p a )) · (1 − S(τ n , p c )) > 0 then substitute p a := p c and go to the step 5 or if (1 − S(τ n , p b )) · (1 − S(τ n , p c )) > 0 then substitute p b := p c and go to the step 5.
The discussed algorithm was used to determine the value of the parameter p, whose values are presented in Table 2.

Numerical examples
In order to validate the results obtained by the proposed numerical method with the analytical solution one, twelve computer simulations were performed. We assumed four values of order of the Caputo derivative α ∈ {0.25, 0.5, 0.75, 1} and three sets of physical parameter values. Due to the fact that the proposed method is iterative (involves multiple solving of governing equations for different values of parameter p), the following mesh parameters were adopted: l 2 /l 1 = 10, m 1 = 100, m 2 = 500, n = 400. Table 1 shows the values of coefficient p obtained from transcendental equation (22) depending on order of the Caputo derivative and three sets of physical parameters. The   Table 1 indicates, that for a fixed set of physical parameters, p is an increasing function of order α. This means that for α values less than one the melting process is slower than in the case of the classical Stefan problem.
Let us now analyze the impact of the physical parameters on the solution for the fixed value of α. Comparing the first and second row of the Table 1, we notice that p is a decreasing function of λ 2 . This observation results directly from equation (18). For a growing heat a flux in solid zone (λ 2 is increasing), heat balance in the melting layer is preserved only when value of the Caputo derivative of function S decreases. For α = 1 fractional derivative of function S can be interpreted as a velocity of movement of the melting front. The opposite behavior of the melting process is observed by analyzing the first and third row of Table 1 for the fixed value of α -coefficient p is an increasing function of parameter κ 1 .
Results given in Table 2 were obtained by applying the criterion of convergence formulated in (63) and the iterative algorithm discussed in the previous section. Table 3 contains the values of the variable τ for which the liquid phase propagates to the region bounded by S(τ ) = 1. The collected data leads to the following conclusion -changing parameter λ 2 has a greater effect on process dynamics, than the changing parameter κ 1 for the fixed value of α.
In Figure 3,5,7 and 9 we present dimensionless temperature profiles. On each graph we show the solutions obtained numerically (blue, yellow and green line) and the corresponding analytical solutions (black lines) plotted for x ∈ [0, 10]. In two cases    (20) is difficult to calculate for large values of the spatial variable x. This is related to the definition of the Wright function, see formula (5). Therefore, due to time constraints of calculations, we used 500 terms of the series. However, this is not sufficient for large values of the similarity variable x τ α/2 . In Figure 4

Conclusions
In this paper, we have developed a numerical method for solving the two-phase fractional Stefan problem, which can be used as an effective tool for determining a solution alternative to the closed analytical one. The applied special case of front-fixing technique allows us use the mesh with constant spatial and time step, which provides discretization of the Caputo derivative with respect to time. Let us point out that the proposed method has some limitations directly derived from the construction of the new spatial variables requiring the analytical explicit form of the function describing the moving boundary depending on the unknown parameter p. Based on the numerical results presented in Section 5, we can draw the following conclusions: 1. The error generated by the numerical solution of the governing subdiffusion equations is a decreasing function of the order of the Caputo derivative. In the special case for α = 1 ,we observe very good matching of the numerical solution to the analytical one.
2. The relationship between the accuracy of the numerical solution of function S and the order of the Caputo derivative is not as obvious as it was listed in the previous point. This is confirmed by the graphs in the Figures 6 and 8. It should be emphasized that the accuracy of the numerical solution of the function S is the result of many factors, which explains formula (62). First, the numerical solutions of the functions u 1 and u 2 are burdened by different errors. Second, the approximation of the heat flux calculated in the liquid phase is usually more accurate than that calculated for the solid phase, because for the example considered in the paper, there is the following relationship x m 1 ,j − x m 1 −1,j < x 1,j − x 0,j (see formula (62) and Figure 2). In summary, the accuracy of the numerical solution of the function S is not only dependent on the order of Caputo derivative α, but also depends on mesh parameters and the model's physical parameters.