A Numerical Method for the Solution of the Two-Phase Fractional Lamé–Clapeyron–Stefan Problem

Abstract

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.

1. 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\left(t\right)>\sim D_{\alpha }{t}^{\alpha }$, where $<x^2\left(t\right)>$ is a mean squared displacement of the diffusing molecule in the course of time t and $D_{\alpha }$ 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 concerning the moving boundary problems modeled by the anomalous diffusion equation. Basically, published papers deal with three classes of phenomena. As mentioned above, the first concerned the controlled release of a drug from slab matrices [9,10]. A second class of problems relates to mathematical modeling of the thermal conductivity with phase transitions [11,12,13,14] and the last class refers to mathematical modeling a movement of the shoreline in a sedimentary ocean basin [15,16].

Most of the analytical solutions of fractional Stefan problem were obtained for the one-phase case, which is a special case of the two-phase problem. For this reason, we will recall the important results obtained for the one-phase problem. Liu Junyi and Xu Mingyu [8] 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^{\frac{\alpha }{2}}$, $0<\alpha \le 1$. Xicheng Li et al. [17] 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, the function describing the moving boundary is ∼$t^{\frac{\alpha }{\beta }}$, $0<\alpha \le 1$, $1<\beta \le 2$, where $\alpha$ and $\beta$ denote the orders of the fractional derivatives with respect to time and spatial variable, respectively. Voller [12] 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<\beta \le 1$ and for the same fractional derivative with respect to space for the flux of order $0<\alpha \le 1$, he showed that a melting front is described by power function $t^{\frac{\beta }{\alpha +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,\tau )$ is transformed into a two new orthogonal systems $(v_1,\tau )$ and $(v_2,\tau )$ 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.

2. 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 $\alpha \in (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+}^{\alpha }$, is given by the following formula for $Re\left(\alpha \right)>0$:

$$I_{0+}^{\alpha }f\left(t\right):=\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^t\frac{f\left(u\right)du}{{(t-u)}^{1-\alpha }},$$

(1)

where Γ is the Euler gamma function.

The left-sided Caputo derivative of order $\alpha \in (0,1)$ denoted by ${}^c{D}_{0+}^{\alpha }$ is defined via the above left-sided Riemann–Liouville integral ${}^c{D}_{0+}^{\alpha }f\left(t\right):=I_{0+}^{1-\alpha }{f}^{\prime }\left(t\right)$.

Definition 2. 

Let $Re\left(\alpha \right)\in (0,1]$. The left-sided Caputo derivative of order α is given by the formula:

$${}^c{D}_{0+}^{\alpha }f\left(t\right):=\left\{\begin{array}{cc}\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^t\frac{f^{{}^{\prime }}\left(u\right)du}{{(t-u)}^{\alpha }},\hfill & 0<\alpha <1,\hfill \\ \frac{df\left(t\right)}{dt},\hfill & \alpha =1.\hfill \end{array}\right.$$

(2)

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\left(\alpha \right)>0$ and $Re\left(\beta \right)>0$, then the equation

$$I_{0+}^{\alpha }{I}_{0+}^{\beta }f\left(t\right)=I_{0+}^{\alpha +\beta }f\left(t\right)$$

(3)

is satisfied at almost every point $t\in [0,b]$ for $f\left(t\right)\in L_p(0,b)$ where $1\le p\le \infty$. If $\alpha +\beta >1$, then the above relation holds at any point of $[0,b]$.

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\in 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:

$$I_{0+}^{\alpha }{}^c{D}_{0+}^{\alpha }f\left(t\right)=f\left(t\right)-f\left(0\right).$$

(4)

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.

Definition 3. 

Let $\gamma >-1$, $\delta \in \mathbb{C},z\in \mathbb{C}$. The two-parameter Wright function is given as the following series:

$$W(z;\gamma ,\delta ):=\sum _{k=0}^{\infty }\frac{z^k}{k!\mathsf{\Gamma }(\gamma k+\delta )}.$$

(5)

We note that for $\gamma =-\frac{1}{2}$, $\delta =1$ the above two-parameter Wright function becomes complementary error function [26]:

$$W\left(-z;-\frac{1}{2},1\right)=erfc\left(\frac{z}{2}\right).$$

(6)

For $\gamma =-\frac{1}{2}$, $\delta =-\frac{1}{2}$ the two-parameter Wright function can be expressed by the formula [27]:

$$W\left(-z;-\frac{1}{2},-\frac{1}{2}\right)=\frac{1}{\sqrt{\pi }}exp\left(-\frac{z^2}{4}\right).$$

(7)

3. 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\left(t\right)$, 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\ge 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.

Consider the following governing equations of the model

$$\begin{array}{c}\hfill c_1{\rho }_1{}^c{D}_{0+,t}^{\alpha }{U}_1(X,t)=K_1\frac{{\partial }^2{U}_1(X,t)}{\partial X^2},0<X<s\left(t\right)\le l_1,t>0,\end{array}$$

(8)

$$\begin{array}{c}\hfill c_2{\rho }_2{}^c{D}_{0+,t}^{\alpha }{U}_2(X,t)=K_2\frac{{\partial }^2{U}_2(X,t)}{\partial X^2},s\left(t\right)<X<\infty ,t>0,\end{array}$$

(9)

where $U_1$ denotes the temperature in the liquid phase, $U_2$ is the temperature in the solid phase, $s\left(t\right)$ represents the position of the moving boundary, K is the modified thermal conductivity (has SI unit [J · s${}^{-\alpha }·$ m${}^{-1}·$ K${}^{-1}$]), $\rho$ 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:

$$\begin{array}{c}\hfill U_1(0,t)=U_0,U_1(s\left(t\right),t)=U_2(s\left(t\right),t)=U_s,\underset{X\to \infty }{lim}{U}_2(X,t)=U_{\infty },t>0.\end{array}$$

(10)

At $t=0$, the liquid phase does not exist, so we use two initial conditions:

$$U_2(X,0)=U_{\infty },s\left(0\right)=0.$$

(11)

The position of the moving boundary $s\left(t\right)$ is determined by the fractional Stefan condition, which expresses the heat balance in the melting layer:

$$\begin{array}{c}\hfill L{\rho }_1{}^c{D}_{0+,t}^{\alpha }s\left(t\right)={\left.K_2\frac{\partial U_2(X,t)}{\partial X}\right|}_{X=s\left(t\right)}{\left.-K_1\frac{\partial U_1(X,t)}{\partial X}\right|}_{X=s\left(t\right)},\end{array}$$

(12)

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

$$\begin{array}{c}\hfill x=\frac{X}{l_1},\tau =t{\left(\frac{K_0}{c_0{\rho }_1{l}_1^2}\right)}^{\frac{1}{\alpha }},u_1=\frac{U_1-U_s}{U_0-U_s},u_2=\frac{U_2-U_s}{U_0-U_s},S=\frac{s}{l_1},\end{array}$$

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

$$\begin{array}{c}\hfill {}^c{D}_{0+,\tau }^{\alpha }{u}_1(x,\tau )={\kappa }_1\frac{{\partial }^2{u}_1(x,\tau )}{\partial x^2},0<x<S\left(\tau \right)\le 1,\tau >0,\end{array}$$

(13)

$$\begin{array}{c}\hfill {}^c{D}_{0+,\tau }^{\alpha }{u}_2(x,\tau )={\kappa }_2\frac{{\partial }^2{u}_2(x,\tau )}{\partial x^2},S\left(\tau \right)<x<\infty ,\tau >0,\end{array}$$

(14)

supplemented with the boundary conditions

$$\begin{array}{c}\hfill u_1(0,\tau )=1,u_1(S\left(\tau \right),\tau )=u_2(S\left(\tau \right),\tau )=0,\underset{x\to \infty }{lim}{u}_2\left(x,\tau \right)=\frac{U_{\infty }-U_s}{U_0-U_s},\tau >0,\end{array}$$

(15)

initial conditions

$$u_2(x,0)=\frac{U_{\infty }-U_s}{U_0-U_s},$$

(16)

$$S\left(0\right)=0,$$

(17)

and the fractional Stefan condition

$$\begin{array}{c}\hfill {}^c{D}_{0+,\tau }^{\alpha }S\left(\tau \right)={\left.{\lambda }_2\frac{\partial u_2(x,\tau )}{\partial x}\right|}_{x=S\left(\tau \right)}{\left.-{\lambda }_1\frac{\partial u_1(x,\tau )}{\partial x}\right|}_{x=S\left(\tau \right)},\end{array}$$

(18)

where ${\kappa }_1=(K_1/K_0)(c_0/c_1)$, ${\kappa }_2=(K_2/K_0)(c_0/c_2)$, ${\lambda }_1=\frac{(U_0-U_s)K_1{c}_0}{L{K}_0}$ and ${\lambda }_2=\frac{(U_0-U_s)K_2{c}_0}{L{K}_0}$ 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

$$u_1(x,\tau )=1-\frac{W\left(\frac{-x}{\sqrt{{\kappa }_1}{\tau }^{\alpha /2}};-\frac{\alpha }{2},1\right)-1}{W\left(\frac{-p}{\sqrt{{\kappa }_1}};-\frac{\alpha }{2},1\right)-1},$$

(19)

$$\begin{array}{c}\hfill u_2(x,\tau )=\frac{U_{\infty }-U_s}{U_0-U_s}\frac{W\left(\frac{-p}{\sqrt{{\kappa }_2}};-\frac{\alpha }{2},1\right)-W\left(\frac{-x}{\sqrt{{\kappa }_2}{\tau }^{\alpha /2}};-\frac{\alpha }{2},1\right)}{W\left(\frac{-p}{\sqrt{{\kappa }_2}};-\frac{\alpha }{2},1\right)},\end{array}$$

(20)

$$S\left(\tau \right)=p{\tau }^{\alpha /2},$$

(21)

$$\begin{array}{c}\hfill p\frac{\mathsf{\Gamma }(1+\frac{\alpha }{2})}{\mathsf{\Gamma }(1-\frac{\alpha }{2})}=\frac{{\lambda }_2}{\sqrt{{\kappa }_2}}\frac{U_{\infty }-U_s}{U_0-U_s}\frac{W\left(\frac{-p}{\sqrt{{\kappa }_2}};-\frac{\alpha }{2},1-\frac{\alpha }{2}\right)}{W\left(\frac{-p}{\sqrt{{\kappa }_2}};-\frac{\alpha }{2},1\right)}-\frac{{\lambda }_1}{\sqrt{{\kappa }_1}}\frac{W\left(\frac{-p}{\sqrt{{\kappa }_1}};-\frac{\alpha }{2},1-\frac{\alpha }{2}\right)}{W\left(\frac{-p}{\sqrt{{\kappa }_1}};-\frac{\alpha }{2},1\right)-1}.\end{array}$$

(22)

It should be noted that the above solution reduces to the fractional one-phase problem [18,19] for $u_2\to 0$.

When $\alpha \to 1$, then from (19)–(22) the following results of the classical two-phase Stefan problem are recovered:

$$u_1(x,\tau )=1-\frac{erfc\left(\frac{x}{2\sqrt{{\kappa }_1\tau }}\right)-1}{erfc\left(\frac{p}{2\sqrt{{\kappa }_1}}\right)-1},$$

(23)

$$u_2(x,\tau )=\frac{U_{\infty }-U_s}{U_0-U_s}\frac{erfc\left(\frac{p}{2\sqrt{{\kappa }_2}}\right)-erfc\left(\frac{x}{2\sqrt{{\kappa }_2\tau }}\right)}{erfc\left(\frac{p}{2\sqrt{{\kappa }_2}}\right)},$$

(24)

$$S\left(\tau \right)=p\sqrt{\tau },$$

(25)

$$\frac{p}{2}={\lambda }_2\frac{U_{\infty }-U_s}{U_0-U_s}\frac{exp\left(-\frac{p^2}{4{\kappa }_2}\right)}{\sqrt{\pi {\kappa }_2}erfc\left(\frac{p}{2\sqrt{{\kappa }_2}}\right)}-{\lambda }_1\frac{exp\left(-\frac{p^2}{4{\kappa }_1}\right)}{\sqrt{\pi {\kappa }_1}\left(erfc\left(\frac{p}{2\sqrt{{\kappa }_1}}\right)-1\right)}.$$

(26)

4. 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 $\tau =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\left(\tau \right)$, marked in blue in Figure 2. For Equation (13), it is convenient to replace the spatial variable with the following similarity variable [28]:

$$v_1=\frac{x}{S\left(\tau \right)}=\frac{x}{p{\tau }^{\alpha /2}},$$

(27)

which has an important property, namely it fixes the moving boundary at $v_1=1$ for all $\tau$. We transform the first-order derivative with respect to the spatial variable:

$$\frac{\partial u_1(x,\tau )}{\partial x}=\frac{\partial v_1}{\partial x}\frac{\partial u_1(v_1,\tau )}{\partial v_1}=\frac{1}{p{\tau }^{\alpha /2}}\frac{\partial u_1(v_1,\tau )}{\partial v_1}.$$

(28)

Subsequently, for the second-order spatial derivative, we have:

$$\frac{{\partial }^2{u}_1(x,\tau )}{\partial x^2}=\frac{1}{p^2{\tau }^{\alpha }}\frac{{\partial }^2{u}_1(v_1,\tau )}{\partial v_1^2}.$$

(29)

For the first-order partial derivative in time, we have:

$$\frac{\partial u_1(x,\tau )}{\partial \tau }=\frac{\partial u_1(v_1,\tau )}{\partial \tau }-\frac{\alpha v_1}{2\tau }\frac{\partial u_1(v_1,\tau )}{\partial v_1}.$$

(30)

The Caputo derivative of function $u_1(v_1,\tau )$ with respect to the time variable can be expressed as follows:

$$\begin{array}{cc}\hfill {}^c{D}_{0+,\tau }^{\alpha }{u}_1(v_1,\tau )& =\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\tau }{(\tau -\xi )}^{-\alpha }\frac{\partial }{\partial \xi }{u}_1(v_1,\xi )d\xi \hfill \\ & -\frac{\alpha v_1}{2\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\tau }{(\tau -\xi )}^{-\alpha }\frac{1}{\xi }\frac{\partial }{\partial v_1}{u}_1(v_1,\xi )d\xi .\hfill \end{array}$$

(31)

Using Formulas (29) and (31), we write Equation (13) in coordinate system $(v_1,\tau )$:

$$\begin{array}{cc}\hfill & \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\tau }\frac{1}{{(\tau -\xi )}^{\alpha }}\frac{\partial u_1(v_1,\xi )}{\partial \xi }d\xi -\hfill \\ \hfill & \frac{\alpha v_1}{2\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\tau }\frac{1}{{(\tau -\xi )}^{\alpha }}\frac{1}{\xi }\frac{\partial u_1(v_1,\xi )}{\partial v_1}d\xi =\frac{{\kappa }_1}{p^2{\tau }^{\alpha }}\frac{{\partial }^2{u}_1(v_1,\tau )}{\partial v_1^2}.\hfill \end{array}$$

(32)

We integrate the previous equation applying the left-sided Riemann–Liouville integral of order $\alpha \in (0,1)$ and we get the following equation:

$$\begin{array}{c}\hfill I_{0+,\tau }^{\alpha }{}^c{D}_{0+,\tau }^{\alpha }{u}_1(v_1,\tau )-\frac{\alpha v_1}{2}{I}_{0+,\tau }^{\alpha }{I}_{0+,\tau }^{1-\alpha }\left(\frac{\partial }{\partial v_1}\frac{u_1(v_1,\tau )}{\tau }\right)=\frac{{\kappa }_1}{p^2}{I}_{0+,\tau }^{\alpha }\left(\frac{1}{{\tau }^{\alpha }}\frac{{\partial }^2{u}_1(v_1,\tau )}{\partial v_1^2}\right)\end{array}.$$

(33)

Finally, using Properties 1 and 2 in Equation (33), we obtain an integro-differential equation in the form of:

$$\begin{array}{cc}\hfill u_1(v_1,\tau )& =u_1(v_1,0)+\frac{\alpha v_1}{2}{\int }_0^{\tau }\frac{\partial }{\partial v_1}\frac{u_1(v_1,\xi )}{\xi }d\xi +\hfill \\ & \frac{{\kappa }_1}{p^2\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^{\tau }\frac{1}{{(\tau -\xi )}^{1-\alpha }}\frac{1}{{\xi }^{\alpha }}\frac{{\partial }^2{u}_1(v_1,\xi )}{\partial v_1^2}d\xi \hfill \end{array}.$$

(34)

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:

$${\overline{u}}_1(v_1,\tau )=u_1(v_1,\tau ){\tau }^{-\alpha },$$

(35)

which leads to the integro-differential equation:

$$\begin{array}{cc}\hfill {\overline{u}}_1(v_1,\tau ){\tau }^{\alpha }& ={\overline{u}}_1(v_1,0){\tau }_0^{\alpha }+\frac{\alpha v_1}{2}{\int }_0^{\tau }\frac{\partial }{\partial v_1}\frac{{\overline{u}}_1(v_1,\xi )}{{\xi }^{1-\alpha }}d\xi +\hfill \\ & \frac{{\kappa }_1}{p^2\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^{\tau }\frac{1}{{(\tau -\xi )}^{1-\alpha }}\frac{{\partial }^2{\overline{u}}_1(v_1,\xi )}{\partial v_1^2}d\xi ,\hfill \end{array}$$

(36)

supplemented with the boundary conditions

$${\overline{u}}_1(0,\tau )={\tau }^{-\alpha },{\overline{u}}_1(1,\tau )=0,$$

and initial condition

$${\overline{u}}_1(v_1,0)=0.$$

The region marked in blue in Figure 2 in the $(x,\tau )$ plane is transformed into the unit square (in the general case to the rectangle) using the transformation (27). The rectangular mesh created in this way consists of horizontal lines spaced $\Delta \tau =1/\left(n{p}^{2/\alpha }\right)$ units apart and vertical lines spaced $\Delta v_1=1/m_1$ units apart. The points $({\left(v_1\right)}_i,{\tau }_j)=(i\Delta v_1,j\Delta \tau )$ where $i=0,\dots ,m_1$ and $j=0,\dots ,n$, of intersection of the horizontal and vertical lines are grid nodes. We denote the value of function ${\overline{u}}_1$ at point $(i\Delta v_1,j\Delta \tau )$ by ${\left({\overline{u}}_1\right)}_{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 $\tau$. Let ${\tau }_0$ be a very small positive number.

The method of discretization of the integro-differential Equation (36) is described in detail in [23] and leads to an implicit scheme:

$$\begin{array}{cc}& {\left({\overline{u}}_1\right)}_{i-1,k+1}(-r_{k+1,k+1}^1+q_{i,k+1}^1)+{\left({\overline{u}}_1\right)}_{i,k+1}({\tau }_{k+1}^{\alpha }+2r_{k+1,k+1}^1)+{\left({\overline{u}}_1\right)}_{i+1,k+1}(-r_{k+1,k+1}^1-q_{i,k+1}^1)\hfill \\ & ={\left({\overline{u}}_1\right)}_{i,0}{\tau }_0^{\alpha }+\sum _{j=0}^k{r}_{j,k+1}^1({\left({\overline{u}}_1\right)}_{i-1,j}-2{\left({\overline{u}}_1\right)}_{i,j}+{\left({\overline{u}}_1\right)}_{i+1,j})+\sum _{j=1}^k{q}_{i,j}^1({\left({\overline{u}}_1\right)}_{i+1,j}-{\left({\overline{u}}_1\right)}_{i-1,j}),\hfill \end{array}$$

(37)

where

$$r_{j,k+1}^1=\frac{c_{j,k+1}{\kappa }_1}{p^2\mathsf{\Gamma }\left(\alpha \right){(\Delta v_1)}^2},q_{i,j}^1=\frac{\alpha i{\tau }_j^{\alpha -1}\Delta \tau }{4},$$

(38)

$$c_{j,k+1}=\left\{\begin{array}{cc}\frac{{(\Delta \tau )}^{\alpha }}{\alpha (\alpha +1)}\left(k^{\alpha +1}-(k-\alpha ){(k+1)}^{\alpha }\right),\hfill & \mathrm{for}j=0,\hfill \\ \frac{{(\Delta \tau )}^{\alpha }}{\alpha (\alpha +1)}{((k-j+2)}^{\alpha +1}+{(k-j)}^{\alpha +1}\hfill & \\ -2{(k-j+1)}^{\alpha +1}),\hfill & \mathrm{for}1\le j\le k,\hfill \\ \frac{{(\Delta \tau )}^{\alpha }}{\alpha (\alpha +1)},\hfill & \mathrm{for}j=k+1.\hfill \end{array}\right.$$

(39)

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

$$\mathbf{A}{\left({\overline{\mathbf{U}}}_{\mathbf{1}}\right)}_{k+1}=\mathbf{B},$$

(40)

where ${\overline{\mathbf{U}}}_{\mathbf{1}}$ is a vector of unknown values of function ${\overline{u}}_1$ at instant ${\tau }_{k+1}$. The elements of matrices $\mathbf{A}$ and $\mathbf{B}$ are defined as follows:

$$\mathbf{A}={\left[\begin{array}{cccccccc}{a}_{k+1}^2& a_{1,k+1}^3& 0& 0& \cdots & 0& 0& 0\\ a_{2,k+1}^1& a_{k+1}^2& a_{2,k+1}^3& 0& \cdots & 0& 0& 0\\ 0& a_{3,k+1}^1& a_{k+1}^2& a_{3,k+1}^3& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& a_{i,k+1}^1& a_{k+1}^2& a_{i,k+1}^3& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& \cdots & a_{m_1-2,k+1}^1& a_{k+1}^2& a_{m_1-2,k+1}^3\\ 0& 0& 0& 0& \cdots & 0& a_{m_1-1,k+1}^1& a_{k+1}^2\end{array}\right]}_{(m_1-1)\times (m_1-1)}$$

$$\mathbf{B}={\left[\begin{array}{c}{b}_1-a_{1,k+1}^1{\left({\overline{u}}_1\right)}_{0,k+1}\\ b_2\\ b_3\\ ⋮\\ b_i\\ ⋮\\ b_{m_1-2}\\ b_{m_1-1}-a_{m_1-1,k+1}^3{\left({\overline{u}}_1\right)}_{m_1,k+1}\end{array}\right]}_{(m_1-1)\times 1}$$

where

$$\begin{array}{ccc}& & a_{i,k+1}^1:=-r_{k+1,k+1}^1+q_{i,k+1}^1,\hfill \\ & & a_{k+1}^2:={\tau }_{k+1}^{\alpha }+2r_{k+1,k+1}^1,\hfill \\ & & a_{i,k+1}^3:=-r_{k+1,k+1}^1-q_{i,k+1}^1,\hfill \\ & & b_i:={\left({\overline{u}}_1\right)}_{i,0}{\tau }_0^{\alpha }+\sum _{j=0}^k{r}_{j,k+1}^1({\left({\overline{u}}_1\right)}_{i-1,j}-2{\left({\overline{u}}_1\right)}_{i,j}+{\left({\overline{u}}_1\right)}_{i+1,j}),\hfill \\ & & +\sum _{j=1}^k{q}_{i,j}^1({\left({\overline{u}}_1\right)}_{i+1,j}-{\left({\overline{u}}_1\right)}_{i-1,j}).\hfill \end{array}$$

We recover the values of function $u_1$ in coordinate system $(x,\tau )$ by applying formulas

$${\left(u_1\right)}_{i,j}={\left({\overline{u}}_1\right)}_{i,j}{\tau }_j^{\alpha },$$

(41)

$$x_{i,j}={\left(v_1\right)}_{i}p{\tau }_j^{\alpha /2}.$$

(42)

The region marked in red in Figure 2 in the $(x,\tau )$ plane is transformed into the unit square using the following transformation:

$$v_2=\frac{x-p{\tau }^{\alpha /2}}{l_2/l_1-p{\tau }^{\alpha /2}},$$

(43)

which fixes the moving boundary at $v_2=0$ for all $\tau$. 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:

$$u_2\left(l_2/l_1,\tau \right)=\frac{U_{\infty }-U_s}{U_0-U_s},$$

(44)

where $l_2$ is a sufficiently large positive real number.

We transform the first-order derivative with respect to the spatial variable:

$$\frac{\partial u_2(x,\tau )}{\partial x}=\frac{\partial v_2}{\partial x}\frac{\partial u_2(v_2,\tau )}{\partial v_2}=\frac{1}{l_2/l_1-p{\tau }^{\alpha /2}}\frac{\partial u_2(v_2,\tau )}{\partial v_2}.$$

(45)

Consequently, for the second-order spatial derivative, we have:

$$\frac{{\partial }^2{u}_2(x,\tau )}{\partial x^2}=\frac{1}{{(l_2/l_1-p{\tau }^{\alpha /2})}^2}\frac{{\partial }^2{u}_2(v_2,\tau )}{\partial v_2^2}.$$

(46)

For the first-order partial derivative with respect to the time variable, we have:

$$\frac{\partial u_2(x,\tau )}{\partial \tau }=\frac{\partial u_2(v_2,\tau )}{\partial \tau }-\frac{\alpha p(v_2-1)}{2\tau (p-l_2/l_1{\tau }^{-\alpha /2})}\frac{\partial u_2(v_2,\tau )}{\partial v_2}.$$

(47)

The Caputo derivative of function $u_1(v_1,\tau )$ with respect to the time variable can be expressed as follows:

$$\begin{array}{cc}\hfill & {}^c{D}_{0+,\tau }^{\alpha }{u}_2(v_2,\tau )=\frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\tau }{(\tau -\xi )}^{-\alpha }\frac{\partial u_2(v_2,\xi )}{\partial \xi }d\xi \hfill \\ & -\frac{\alpha p(v_2-1)}{2\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\tau }\frac{{(\tau -\xi )}^{-\alpha }}{\xi (p-l_2/l_1{\xi }^{-\alpha /2})}\frac{\partial u_2(v_2,\xi )}{\partial v_2}d\xi .\hfill \end{array}$$

(48)

Using Formulas (46) and (48), we can write Equation (14) in the new coordinate system $(v_2,\tau )$:

$$\begin{array}{cc}\hfill & \frac{1}{\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\tau }\frac{1}{{(\tau -\xi )}^{\alpha }}\frac{\partial u_2(v_2,\xi )}{\partial \xi }d\xi -\frac{\alpha p(v_2-1)}{2\mathsf{\Gamma }(1-\alpha )}{\int }_0^{\tau }\frac{{(\tau -\xi )}^{-\alpha }}{\xi (p-l_2/l_1{\xi }^{-\alpha /2})}\frac{\partial u_2(v_2,\xi )}{\partial v_2}d\xi =\hfill \\ & \frac{{\kappa }_2}{{(l_2/l_1-p{\tau }^{\alpha /2})}^2}\frac{{\partial }^2{u}_2(v_2,\tau )}{\partial v_2^2}.\hfill \end{array}$$

(49)

We integrate both sites of Equation (49) applying the left-sided Riemann–Liouville integral of order $\alpha \in (0,1)$:

$$\begin{array}{cc}\hfill & I_{0+,\tau }^{\alpha }{}^c{D}_{0+,\tau }^{\alpha }{u}_2(v_2,\tau )-\frac{\alpha p(v_2-1)}{2}{I}_{0+,\tau }^{\alpha }{I}_{0+,\tau }^{1-\alpha }\left(\frac{\partial }{\partial v_2}\frac{u_2(v_2,\tau )}{\tau (p-l_2/l_1{\tau }^{-\alpha /2})}\right)=\hfill \\ & I_{0+,\tau }^{\alpha }\left(\frac{{\kappa }_2}{{(l_2/l_1-p{\tau }^{\alpha /2})}^2}\frac{{\partial }^2{u}_2(v_2,\tau )}{\partial v_2^2}\right).\hfill \end{array}$$

(50)

Finally, using Properties 1 and 2 in Equation (50), we obtain an integro-differential equation in the form of:

$$\begin{array}{cc}\hfill u_2(v_2,\tau )& =u_2(v_2,0)+\frac{\alpha p(v_2-1)}{2}{\int }_0^{\tau }\frac{\partial }{\partial v_2}\frac{u_2(v_2,\xi )}{\xi (p-l_2/l_1{\xi }^{-\alpha /2})}d\xi +\hfill \\ & \frac{{\kappa }_2}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^{\tau }\frac{{(\tau -\xi )}^{\alpha -1}}{{(l_2/l_1-p{\xi }^{\alpha /2})}^2}\frac{{\partial }^2{u}_2(v_2,\xi )}{\partial v_2^2}d\xi .\hfill \end{array}$$

(51)

We introduce the second auxiliary function

$${\overline{u}}_2(v_2,\tau )=\frac{u_2(v_2,\tau )}{{(l_2/l_1-p{\tau }^{\alpha /2})}^2},$$

(52)

which leads to the integro-differential equation:

$$\begin{array}{cc}\hfill & {\overline{u}}_2(v_2,\tau ){(l_2/l_1-p{\tau }^{\alpha /2})}^2={\overline{u}}_2(v_1,0){(l_2/l_1-p{\tau }_0^{\alpha /2})}^2+\hfill \\ & \frac{\alpha p(v_2-1)}{2}{\int }_0^{\tau }\frac{\partial }{\partial v_2}\left(p{\xi }^{\alpha -1}-\frac{l_2}{l_1}{\xi }^{\alpha /2-1}\right){\overline{u}}_2(v_2,\xi )d\xi +\hfill \\ & \frac{{\kappa }_2}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^{\tau }\frac{1}{{(\tau -\xi )}^{1-\alpha }}\frac{{\partial }^2{\overline{u}}_2(v_2,\xi )}{\partial v_2^2}d\xi ,\hfill \end{array}$$

(53)

supplemented with the boundary conditions

$${\overline{u}}_2(0,\tau )=0,{\overline{u}}_2(1,\tau )=\frac{1}{{(l_2/l_1-p{\tau }^{\alpha /2})}^2}\frac{U_{\infty }-U_s}{U_0-U_s},$$

(54)

and initial condition

$${\overline{u}}_2(v_2,0)=\frac{1}{{(l_2/l_1-p{\tau }_0^{\alpha /2})}^2}\frac{U_{\infty }-U_s}{U_0-U_s}.$$

(55)

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 $\Delta \tau =1/\left(n{p}^{2/\alpha }\right)$ units apart and vertical lines spaced $\Delta v_2=1/m_2$ units apart. The points $({\left(v_2\right)}_i,{\tau }_j)=(i\Delta v_2,j\Delta \tau )$ where $i=0,\dots ,m_2$ and $j=0,\dots ,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:

$$\begin{array}{cc}\hfill & {\left({\overline{u}}_2\right)}_{i-1,k+1}(-r_{k+1,k+1}^2+q_{i,k+1}^2)+{\left({\overline{u}}_2\right)}_{i,k+1}({(l_2/l_1-p{\tau }_{k+1}^{\alpha /2})}^2\hfill \\ & +2r_{k+1,k+1}^2)+{\left({\overline{u}}_2\right)}_{i+1,k+1}(-r_{k+1,k+1}^2-q_{i,k+1}^2)=\hfill \\ & ={\left({\overline{u}}_1\right)}_{i,0}{(l_2/l_1-p{\tau }_0^{\alpha /2})}^2+\sum _{j=0}^k{r}_{j,k+1}^2({\left({\overline{u}}_2\right)}_{i-1,j}-2{\left({\overline{u}}_2\right)}_{i,j}+{\left({\overline{u}}_2\right)}_{i+1,j})+\hfill \\ & \sum _{j=1}^k{q}_{i,j}^2({\left({\overline{u}}_2\right)}_{i+1,j}-{\left({\overline{u}}_2\right)}_{i-1,j}),\hfill \end{array}$$

(56)

where

$$r_{j,k+1}^2=\frac{c_{j,k+1}{\kappa }_2}{\mathsf{\Gamma }\left(\alpha \right){(\Delta v_2)}^2},q_{i,j}^2=\frac{\alpha p(i\Delta v_2-1)(p{\tau }_j^{\alpha -1}-l_2/l_1{\tau }_j^{\alpha /2-1})\Delta \tau }{4\Delta v_2}.$$

(57)

The obtained numerical scheme can be written in the matrix form

$$\mathbf{D}{\left({\overline{\mathbf{U}}}_{\mathbf{2}}\right)}_{k+1}=\mathbf{E},$$

(58)

where ${\overline{\mathbf{U}}}_{\mathbf{2}}$ is a vector of unknown values of function ${\overline{u}}_2$ at instant ${\tau }_{k+1}$.

$$\mathbf{D}={\left[\begin{array}{cccccccc}{d}_{k+1}^2& d_{1,k+1}^3& 0& 0& \cdots & 0& 0& 0\\ d_{2,k+1}^1& d_{k+1}^2& d_{2,k+1}^3& 0& \cdots & 0& 0& 0\\ 0& d_{3,k+1}^1& d_{k+1}^2& d_{3,k+1}^3& \cdots & 0& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& d_{i,k+1}^1& d_{k+1}^2& d_{i,k+1}^3& 0& 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮& ⋮\\ 0& 0& 0& 0& \cdots & d_{m_2-2,k+1}^1& d_{k+1}^2& d_{m_2-2,k+1}^3\\ 0& 0& 0& 0& \cdots & 0& d_{m_2-1,k+1}^1& d_{k+1}^2\end{array}\right]}_{(m_2-1)\times (m_2-1)}$$

$$\mathbf{E}={\left[\begin{array}{c}{e}_1-d_{1,k+1}^1{\left({\overline{u}}_2\right)}_{0,k+1}\\ e_2\\ e_3\\ ⋮\\ e_i\\ ⋮\\ e_{m_2-2}\\ e_{m_2-1}-d_{m_2-1,k+1}^3{\left({\overline{u}}_2\right)}_{m_2,k+1}\end{array}\right]}_{(m_2-1)\times 1}$$

The elements of matrices $\mathbf{D}$ and $\mathbf{E}$ are defined as follows:

$$\begin{array}{ccc}& & d_{i,k+1}^1:=-r_{k+1,k+1}^2+q_{i,k+1}^2,\hfill \\ & & d_{k+1}^2:={(l_2/l_1-p{\tau }_{k+1}^{\alpha /2})}^2+2r_{k+1,k+1}^2,\hfill \\ & & d_{i,k+1}^3:=-r_{k+1,k+1}^2-q_{i,k+1}^2,\hfill \\ & & e_i:={\left({\overline{u}}_2\right)}_{i,0}{(l_2/l_1-p{\tau }_0^{\alpha /2})}^2,+\sum _{j=0}^k{r}_{j,k+1}^2({\left({\overline{u}}_2\right)}_{i-1,j}-2{\left({\overline{u}}_2\right)}_{i,j}+{\left({\overline{u}}_2\right)}_{i+1,j})\hfill \\ & & +\sum _{j=1}^k{q}_{i,j}^2({\left({\overline{u}}_2\right)}_{i+1,j}-{\left({\overline{u}}_2\right)}_{i-1,j}).\hfill \end{array}$$

Applying the following formulas:

$${\left(u_2\right)}_{i,j}={\left({\overline{u}}_2\right)}_{i,j}{(l_2/l_1-p{\tau }_j^{\alpha /2})}^2,$$

(59)

$$x_{i,j}={\left(v_2\right)}_i(l_2/l_1-p{\tau }_j^{\alpha /2})+p{\tau }_j^{\alpha /2},$$

(60)

we recover the values of function $u_2$ in coordinate system $(x,\tau )$.

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 $\alpha \in (0,1)$. Next, using Property 2, we obtain:

$$\begin{array}{c}\hfill S\left(\tau \right)-S\left(0\right)={\left.{\lambda }_2{I}_{0+,\tau }^{\alpha }\frac{\partial u_2(x,\tau )}{\partial x}\right|}_{x=S\left(\tau \right)}{\left.-{\lambda }_1{I}_{0+,\tau }^{\alpha }\frac{\partial u_1(x,\tau )}{\partial x}\right|}_{x=S\left(\tau \right)}.\end{array}$$

(61)

Applying the difference quotient approximation for the first order derivative with respect to space variable, weights (39) and condition (17), we get:

$$S\left({\tau }_n\right)=\frac{{\lambda }_2}{\mathsf{\Gamma }\left(\alpha \right)}\sum _{j=0}^n{c}_{j,n}\frac{{\left(u_2\right)}_{1,j}-{\left(u_2\right)}_{0,j}}{x_{1,j}-x_{0,j}}-\frac{{\lambda }_1}{\mathsf{\Gamma }\left(\alpha \right)}\sum _{j=0}^n{c}_{j,n}\frac{{\left(u_1\right)}_{m_1,j}-{\left(u_1\right)}_{m_1-1,j}}{x_{m_1,j}-x_{m_1-1,j}},j=0,\dots ,n.$$

(62)

Let us note that the value of function S at the final time instant ${\tau }_n$ for the correct value of parameter p should be equal to 1. From this relation, the following criterion of convergence results in:

$$|1-S({\tau }_n,p\left)\right|<ϵ,$$

(63)

where $ϵ>0$ is a certain arbitrarily small real number. We denote as $S({\tau }_n,p)$ the value of $S\left({\tau }_n\right)$ 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].

• Choose interval $[p_a,p_b]$ for parameter p, define the value of parameter $ϵ$
• Calculate numerically the solution of Equations (36) and (53) for $p_a$ and $p_b$
• If one of the conditions
  $|1-S({\tau }_n,p_a\left)\right|<ϵ$ or $|1-S({\tau }_n,p_b\left)\right|<ϵ$
  is fulfilled, then end the algorithm, otherwise go to the next step
• If the condition
  $(1-S({\tau }_n,p_a))·(1-S({\tau }_n,p_b))<0$
  is fulfilled, then go to the next step, otherwise go to Step 1
• Calculate $p_c:=\frac{p_a+p_b}{2}$
• Calculate numerically the solution of Equations (36) and (53) for $p_c$
• If $|1-S({\tau }_n,p_c\left)\right|<ϵ$
  is fulfilled, then end the algorithm, otherwise go to the next step
• If $(1-S({\tau }_n,p_a))·(1-S({\tau }_n,p_c))>0$,
  then substitute $p_a:=p_c$ and go to Step 5, or
  if $(1-S({\tau }_n,p_b))·(1-S({\tau }_n,p_c))>0$,
  then substitute $p_b:=p_c$ and go to Step 5.

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.

5. 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 $\alpha \in \{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. Value of parameter p from exact solution.

				$\mathit{\alpha }$
${\mathit{\lambda }}_{\mathbf{1}}$	${\mathit{\lambda }}_{\mathbf{2}}$	${\mathit{\kappa }}_{\mathbf{1}}$	${\mathit{\kappa }}_{\mathbf{2}}$	0.25	0.5	0.75	1
1	1	1	1	0.6834	0.7472	0.8299	0.9397
1	2	1	1	0.5496	0.6013	0.668	0.7555
1	1	2	1	0.7218	0.7868	0.8697	0.9783

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 $\alpha$. This means that for $\alpha$ 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 $\alpha$. Comparing the first and second rows of Table 1, we notice that p is a decreasing function of ${\lambda }_2$. This observation results directly from Equation (18). For a growing heat a flux in solid zone (${\lambda }_2$ is increasing), heat balance in the melting layer is preserved only when the value of the Caputo derivative of function S decreases. For $\alpha =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 $\alpha$, coefficient p is an increasing function of parameter ${\kappa }_1$.

The results given in

Table 2. Value of parameter p from numerical solution.

				$\mathit{\alpha }$
${\mathit{\lambda }}_{\mathbf{1}}$	${\mathit{\lambda }}_{\mathbf{2}}$	${\mathit{\kappa }}_{\mathbf{1}}$	${\mathit{\kappa }}_{\mathbf{2}}$	0.25	0.5	0.75	1
1	1	1	1	0.7053	0.7358	0.8041	0.9311
1	2	1	1	0.5691	0.5935	0.6497	0.7512
1	1	2	1	0.739	0.7801	0.8514	0.9674

were obtained by applying the criterion of convergence formulated in (63) and the iterative algorithm discussed in the previous section.

Table 3. Value of variable $\tau$ fulfilling equation $S\left(\tau \right)=1$ obtained from numerical solution.

				$\mathit{\alpha }$
${\mathit{\lambda }}_{\mathbf{1}}$	${\mathit{\lambda }}_{\mathbf{2}}$	${\mathit{\kappa }}_{\mathbf{1}}$	${\mathit{\kappa }}_{\mathbf{2}}$	0.25	0.5	0.75	1
1	1	1	1	16.329	3.411	1.789	1.153
1	2	1	1	90.885	8.06	3.158	1.772
1	1	2	1	11.241	2.7	1.536	1.069

contains the values of the variable $\tau$ for which the liquid phase propagates to the region bounded by $S\left(\tau \right)=1$. The collected data lead to the following conclusion: changing parameter ${\lambda }_2$ has a greater effect on process dynamics than the changing parameter ${\kappa }_1$ for the fixed value of $\alpha$.

In

Table 4. Average absolute error generated by front-fixing method.

				$\mathit{\alpha }$
${\mathit{\lambda }}_{\mathbf{1}}$	${\mathit{\lambda }}_{\mathbf{2}}$	${\mathit{\kappa }}_{\mathbf{1}}$	${\mathit{\kappa }}_{\mathbf{2}}$	0.25	0.5	0.75	1
1	1	1	1	0.0109	0.0082	0.0076	0.0012
1	2	1	1	0.0108	0.0079	0.0064	0.0009
1	1	2	1	0.0106	0.0081	0.0059	0.0015

, we present the average absolute error generated by the proposed method in mesh nodes fulfilling condition $\frac{x_i}{{\tau }_j^{\alpha /2}}<10$. The results collected in Table 4 lead to the following conclusion: the average error is a decreasing function of the order $\alpha$.

In Figure 3, Figure 4, Figure 5 and Figure 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\in [0,10]$. In two cases, for $\alpha =1$ and $\alpha =0.75$, we adopted $x\in [0,2]$ and $x\in [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 $\frac{x}{{\tau }^{\alpha /2}}$.

In Figure 7, Figure 8, Figure 9 and Figure 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.

Figure 11, Figure 12, Figure 13 and Figure 14 lead to an important observation. In all considered cases, we note the highest absolute error for small values of variable $\tau$ and in the area close to the moving boundary $x=1$. The calculations were performed in the mesh nodes fulfilling condition $\frac{x_i}{{\tau }_j^{\alpha /2}}<10$.

6. 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(\Delta \tau )$) and the central differential quotient (error $O(\Delta v_1^2)$ and $O(\Delta v_2^2)$). The second integral term of Equations (36) and (53) is approximated by the fractional trapezoidal rule (error $O(\Delta {\tau }^{1+\alpha })$ [29]) and the difference quotient for the second derivative (error $O(\Delta v_1^2)$ and $O(\Delta 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 $\alpha$ and the accuracy of the method.

Figure 11, Figure 12, Figure 13 and Figure 14 show the distribution of the absolute errors for the four cases presented in the paper. Their analysis leads to the following observations:

• The transformation (27) has a singularity at $\tau =0$. Applying the numerical approximation consisting in assuming a very small positive number as ${\tau }_0$ (the calculations assumed ${\tau }_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.
• 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 Figure 11, Figure 12, Figure 13 and Figure 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.

7. 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.