A Numerical Method for the Solution of the Two-Phase Fractional Lamé–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. The presented method is an extension of a front-fixing method for the one-phase problem to the two-phase case. The novelty of the method is a new discretization of the partial differential equation dedicated to the second phase, which is carried out by introducing a new spatial variable immobilizing the moving boundary. Then, the partial differential equation is transformed to an equivalent integro-differential equation, which is discretized on a homogeneous mesh of nodes with a constant spatial and time step. A new convergence criterion is also proposed in the iterative algorithm determining the location of the moving boundary. The motivation for the development of the method is that the analytical solution of the considered problem is impossible to calculate in some cases, as can be seen in the figures in the paper. Moreover, the change of the boundary conditions makes obtaining a closed analytical solution very problematic. Therefore, creating new numerical methods is very valuable. 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
Differential equations are an excellent tool for describing problems with technical applications [1]. A special class of differential equations are the so-called moving boundary problems which are 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 [2]. In recent years, the classical Stefan problem has been very well studied and described in many papers and monographs (compare [3][4][5][6][7] and references therein).
The fractional Stefan problem is a natural generalization of the classical Stefan problem and is related to the anomalous mass or heat transfer usually occurs in materials with strong heterogeneities across a range of length scales. The movement of molecules in such materials is described by the law < x 2 (t) >∼ D α t α , where < x 2 (t) > is a mean squared displacement of the diffusing molecule in the course of time t and D α is the generalized diffusion coefficient. 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 [8]. Recently, there has been an increase in the number of scientific publications flux of order 0 < α ≤ 1, he showed that a melting front is described by power function t β α+1 . More analytical results are published in [18][19][20].
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 [21], which is recalled in Section 3. The problem involves determination of three functions, namely, u 1 , u 2 and S fulfilling two subdiffusion Equations (13) and (14) and additional differential Equation (18) (interface energy balance condition) governing function S. The authors showed that u 1 and u 2 can be expressed in terms of the Wright's function, while the location of the phase-change interface is described by power function (21), where parameter p can be evaluated form transcendental Equation (22). This solution is 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. [22] generalized the boundary immobilization technique (also known as front-fixing method [3]) to the fractional Stefan problem with a space-fractional derivative. Another variant of the front-fixing method was proposed by Błasik and Klimek [23] 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 front-fixing method [23] 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 integro-differential 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 [23] to the two-phase case. Section 5 contains the analytical and numerical results. The last two sections includes a discussion and conclusions.

Preliminaries
Let us recall the basic definitions and properties from fractional calculus [24,25] which are 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 The left-sided Caputo derivative of order α is given by the formula: An important property of fractional integral operator given by Formula (1) is called a semigroup property, which allows composition of two left-sided fractional integrals. Property 1 (cf. Lemma 2.3 [24]). 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 2 (cf. Lemma 2.22, [24]). 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.

Formulation of the Problem
The two-phase fractional Stefan problem is a mathematical model describing the solidification and melting process. The mathematical formulation of the problem involves an anomalous 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 is the temperature in the solid phase, s(t) represents the position of the moving boundary, K is the modified thermal conductivity (has SI unit [J· s −α · m −1 · K −1 ]), ρ is the density, c is the specific heat and subscripts 1 and 2 indicate liquid and solid phase, respectively. Equations (8) and (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) and (9) in a non-dimensional form supplemented with the boundary conditions initial conditions and the fractional Stefan condition where are the standard values of the respective variables.
According to results obtained by Roscani and Tarzia [21], 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 [18,19] for u 2 → 0 .
When α → 1, then from (19)-(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 [23] 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 Equation (13), it is convenient to replace the spatial variable with the following similarity variable [28]: which has an important property, namely it 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: For the first-order partial derivative in time, we have: The 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 1 and 2 in 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ū 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) The method of discretization of the integro-differential Equation (36) is described in detail in [23] and leads to an implicit scheme: where Weights (39) are the result of application of the trapezoidal rule [29] to the last integral term in Formula (36). The 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, we have: For the first-order partial derivative with respect to the time variable, we have: The Caputo derivative of function u 1 (v 1 , τ) with respect to the time variable can be expressed as follows: 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 1 and 2 in 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) and (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 2, we obtain: Applying the difference quotient approximation for the first order derivative with respect to space variable, weights (39) and condition (17), we get: 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 [23].

1.
Choose interval [p a , p b ] for parameter p, define the value of parameter 2.
Calculate numerically the solution of Equations (36) and (53) for p a and p b 3.
If one of the conditions |1 − S(τ n , p a )| < or |1 − S(τ n , p b )| < is fulfilled, then end the algorithm, otherwise go to the next step 4.
Calculate p c := p a +p b 2 6.
Calculate numerically the solution of Equations (36) and (53) for p c 7.
If |1 − S(τ n , p c )| < is fulfilled, then end the algorithm, otherwise go to the next step 8.
A small value of causes that the algorithm need many iterations to achieve the desired accuracy. In this case, we can define an additional stop criterion. For example, we can study the variance of the value of p from the last few iterations-a value close to zero will show that the assumed accuracy has been achieved.

Numerical Examples
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 and n = 400. Moreover, the calculations assumed = 0.001, p a = p exact − 0.05 and p b = p exact + 0.05, where p exact was determined from the exact solution. The calculations usually required about eight iterations. Table 1 shows the values of coefficient p obtained from transcendental Equation (22) depending on the order of the Caputo derivative and three sets of physical parameters.  The results collected in Table 1 indicate 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 rows of 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 the value of the Caputo derivative of function S decreases. For α = 1, the 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 .
The 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 lead 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 Table 4, we present the average absolute error generated by the proposed method in mesh nodes fulfilling condition x i τ α/2 j < 10. The results collected in Table 4 lead to the following conclusion: the average error is a decreasing function of the order α. Table 4. Average absolute error generated by front-fixing method. In Figures 3-6, 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, for α = 1 and α = 0.75, we adopted x ∈ [0, 2] and x ∈ [0, 6.5], respectively, because the analytical solution (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 Figures 7-10, we present graphs of the function S which describe the position of the phase change front. The black and red linea represent the analytical and numerical solutions, respectively.

Discussion
The proposed numerical method is complicated due to the non-local differential operators used in the equations. The numerical scheme presented in the paper allows to solve three partial differential Equations (13), (14) and (18) supplemented with a set of boundary and initial conditions. These equations are related to each other through the function S, describing the position of the moving boundary, which reflects the presence of the parameter p in the numerical scheme. The analysis of the convergence of the presented method would require the simultaneous consideration of the implicit (multi-step) scheme (40) and (58) and Formula (62), which is a very complicated task. Therefore, let us look at the proposed method not as a whole, but at its smaller fragments that make up this whole.
At this point, it should be noted that the discretization of the integro-differential Equations (36) and (53) is performed using convergent approximations and numerical methods known from the literature. Let us point out that the first integral term of Equations (36) and (53) is approximated by the rectangular rule (error O(∆τ)) and the central differential quotient (error O(∆v 2 1 ) and O(∆v 2 2 )). The second integral term of Equations (36) and (53) is approximated by the fractional trapezoidal rule (error O(∆τ 1+α ) [29]) and the difference quotient for the second derivative (error O(∆v 2 1 ) and O(∆v 2 2 )). The same methods are used in Formula (62). The above observations are in agreement with the results in Table 4, which shows the relationship between the order of the Caputo derivative α and the accuracy of the method. Figures 11-14 show the distribution of the absolute errors for the four cases presented in the paper. Their analysis leads to the following observations: 1.
The transformation (27) has a singularity at τ = 0. Applying the numerical approximation consisting in assuming a very small positive number as τ 0 (the calculations assumed τ 0 = 10 −30 ), we make a certain error at the beginning of the calculations. This initial error is transferred to subsequent time layers of the mesh, by the preceding time layers. Its highest value was observed for the initial time layers. We can significantly improve the quality of the solution using the numerical integration algorithm supported by an artificial neural network, as shown in [30]. This approach is very time-consuming and is planned at a later stage of the research.

2.
Formula (62) allows us to determine the position of the moving boundary for the nodes of the mesh shown on the left side of Figure 2. It is a non-uniform grid with a constant time step and a variable spatial step. The unequal spatial step causes that the flux of the liquid and solid phases are approximated with different accuracy, which results in some inaccuracy in determining the position of the moving boundary, which is visible in Figures 11-14. At a further stage of the research, it is planned to determine the optimal relationships among m 1 , m 2 and n. An artificial neural network will be used to implement this task.

Conclusions
In this paper, we develop 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 to 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.
As part of further research, it is planned to test the stability and convergence of the proposed numerical method.
Funding: This research received no external funding.