On the Computation of the Prefactor of the Free Boundary in One Dimensional One-Phase Fractional Stefan Problems

1. Introduction

The Stefan problem of melting is a well known free boundary problem that has been deeply studied in the last century. The scope of this model is very broad, as can be seen in the review article [26] which cites and compiles more than 5000 articles on Stefan problems up to the year 2000.

When we deal with one-dimensional Stefan problems, under certain initial and boundary conditions, the motion of the free boundary (the melt front) behaves as the square of time t[1,16]. In particular, if we address the one-phase Stefan problem with a constant boundary condition, that is, if we consider the problem to find the pair of functions $\{u,s\}$ representing the temperature and the advance of the interface, such that

$$\begin{array}{ccc}\left(i\right)\hfill & \rho c{u}_t(x,t)=k{u}_{xx}(x,t)\hfill & 0<x<s\left(t\right),0<t<T,\hfill \\ \left(ii\right)\hfill & u(0,t)=U_0>0\hfill & 0<t<T,\hfill \\ \left(iii\right)\hfill & u\left(s\right(t),t)=0\hfill & 0<t<T,\hfill \\ \left(iv\right)\hfill & s\left(0\right)=0,\hfill & \\ \left(v\right)\hfill & \rho \ell \dot{s}\left(t\right)=-\underset{x\to s{\left(t\right)}^{-}}{lim}\hspace{0.17em}{u}_x(x,t)\hfill & 0<t<T,\hfill \end{array}$$

(1)

then the free boundary is given in terms of the square root of time,

$$s\left(t\right)=\xi \sqrt{t},\hspace{0.17em}t\ge 0$$

(2)

where $\xi$ is the positive solution to the transcendental equation

$${\displaystyle \frac{x}{2}}erf\left({\displaystyle \frac{x}{2}}\right)e^{x^2/4}=-{\displaystyle \frac{\rho c}{\sqrt{\pi }}}.$$

We call $\xi$ the prefactor of the free boundary because, clearly, the advance of the interface is totally characterized by this number. The relevance of computing this parameter (which is relatively easy nowadays) was exposed by Solomon in [24] by giving exact data and demonstrating the importance of similarity solutions in semi-infinite domains to approximate exact solutions.

The heat transport governing equation in the Stefan problem above, (1)$-\left(i\right)$, is the classical diffusion equation. This equation arises as a consequence of inserting the phenomenological Fourier law for heat transfer into the heat continuity equation (first law of thermodynamics). However, the Fourier law is not the only way to model heat transfer. The theory of non-Fourier heat transfer is highly diverse, as can be observed, for instance, in Zhmakin’s book [31] or the work [23]. Such treatments are suitable to model heat transfer in the presence of fractal structures, or amorphous materials like glassy polymers [2] or silica glasses [13]. The focus of this article is a particular class of non-Fourier heat transfer models, namely those that employ fractional calculus to describe anomalous diffusion. More explicitly, we consider here non-Fourier models in presence of fractional integrals or derivatives of Riemann-Liouville or Caputo type. We refer the readers to the works of Metzler and Klafter [17] and Voller [29] which provide an in-depth discussion on the modeling of anomalous diffusion processes using fractional derivatives, as well as the books [9,10,15,18,22] for the theory and applications of fractional diffusion equations.

When the presence of anomalous diffusion is considered in a free boundary problem with an imposed constant temperature at the fix boundary, an interesting feature emerges: the advance of the free boundary is proportional to a power of t that differs from $1/2$ (see e.g. [21] or [28]). That is,

$$s_{\alpha }\left(t\right)={\xi }_{\alpha }{t}^{\gamma \left(\alpha \right)},\hspace{0.17em}t\ge 0;$$

(3)

when $\gamma \left(\alpha \right)>1/2$ the process is referred to as superdiffusive and when $\gamma \left(\alpha \right)<1/2$ as subdiffusive.

We are strongly convinced that fractional models can be valuable when emphasizing the evolution of the free boundary. Thus, mirroring classical work on Stefan problems [24], our goal here is to develop algorithms and techniques to compute the prefactor ${\xi }_{\alpha }$ associated with both super- and sub- diffusive one-dimensional, one-phase Stefan melting problems. In carrying out this task we will use the similarity solutions of anomalous Stefan problems that has been previously presented in the literature. In this light, we emphasis that the key component in our work is to arrive at computations that produce explicit numerical values for the prefactors ${\xi }_{\alpha }$. Beyond their intrinsic value, such results are useful in validating the effectiveness and accuracy of numerical methods for phase-change problems in presence of anomalous diffusion, e.g., [4,7,8,11]. As a verification of our proposed computations we provide a mathematical analysis that shows that when the governing dimensionless group, the Stefan number $Ste\to 0$, our calculations for ${\xi }_{\alpha }$ converge to the known explicit values associated with quasi-stationary anomalous Stefan problems [21,28].

Our paper is laid out as follows. In Section 2, the time fractional and space fractional Stefan problems for a melting process are presented. Section 3 is devoted to the a dimensionless scaling of such problems, focusing on the role of the Stefan number in each case. Section 4 and 5 deal with the space fractional and time fractional cases, respectively, having both the same structure. The selfsimilar solutions for the general and the stationary problems are presented first. Then a method for the computation of the prefactor is presented and finally, an analytical proof of the convergence of the prefactors of the general solutions to the prefactors of the quasi-stationary solution when the Stefan number $Ste$ approaches zero, is given. Finally we present the conclusions and future work in Section 6.

2. Modeling the Space and Time Fractional Stefan Problems in a Melting Process

Suppose that $\Omega$ is a cylinder of anomalous phase change material with constant cross-section A and length $[0,l]$, where $l>>1$ is an isolated extreme. Let $u=u(x,t)$ be the temperature of the cylinder at position x and time t and let $q=q(x,t)$ be the heat flux. We assume that, initially, the bar is at its melt temperature $u=U_m$ and a melting process with a sharp moving interface $s=s\left(t\right)$ is initiated by imposing a constant temperature $U_0>U_m$ at $x=0$. This free boundary, representing the moving melt interface, is assumed to be an increasing function in time, i.e,

$$x=s\left(t\right),\mathrm{with}\mathrm{inverse}t=h\left(x\right).$$

(4)

The total energy in the model is given by the enthalpy

$$H(x,t)=\left\{\begin{array}{cc}c\hspace{0.17em}u(x,t)+\ell \hfill & \mathrm{if}\hspace{0.17em}0<x<s\left(t\right),t>0,\hfill \\ c{U}_m\hfill & x>s\left(t\right),t>0,\hfill \end{array}\right.$$

(5)

where the thermophysical parameters at the liquid phase are the specific heat c and latent heat ℓ, which is the energy used in the phase change. With this definition we can use the basic thermodynamic principles to recover two balance statements:

$\left(\mathrm{i}\right)$The balance of energy in the liquid melt$0<x<s$,

$$\rho {\displaystyle \frac{\partial }{\partial t}}H+{\displaystyle \frac{\partial }{\partial x}}q=0,$$

(6)

where $\rho$ is the material density and $\left(\mathrm{ii}\right)$ the of balance energy at the interface where the phase change occurs,

$$〚\rho H{〛}_l^s\mathbf{v}=〚q{〛}_l^s.$$

(7)

This is recognized as the classical Rankine–Hugonoit condition, the double brackets representing jumps in the enthalpy and the heat flux respectively and $\mathbf{v}$ the velocity of the free boundary.

It is well known that equation (1)$-\left(i\right)$ is obtained after replacing the classical Fourier law in the balance equation (6). This law states that the local heat flux q at a point x is proportional to the gradient of temperature

$$q(x,t)=-k{u}_x(x,t),$$

(8)

where we have assumed that the conductivity of the material k is constant. By contrast, here we consider a non-local flux law defined by the relation

$$q_{\alpha }^{nl}(x,t)=-{\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_0^x\hspace{0.17em}{k}_{\alpha }{u}_x(p,t){(x-p)}^{-\alpha }\mathrm{dp}.$$

(9)

This equality states that the flux at every time t and position x is a generalized sum of all the local fluxes at every position between the left extreme of the slab ($x=0$) and the current position, imposing the condition that the local fluxes “closer” to the current position have more weighting than the local fluxes “further” away. Note the constant $k_{\alpha }$ is defined by

$$k_{\alpha }:=k{\nu }_{\alpha }$$

(10)

where ${\nu }_{\alpha }$ is a parameter that has been added to preserve the physical dimensions.

At this point let us recall the classical definitions of fractional calculus that will be used in this article.

For every $f\in L^1(a,b)$ we define the fractional integral of Riemann-Liouville of order $\alpha \in (0,1)$ by

$${}_{a}I^{\alpha }f\left(z\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}{\int }_a^{z}f\left(w\right){(z-w)}^{\alpha -1}\mathrm{d}w,\mathrm{a}.\mathrm{e}.\mathrm{in}(a,b).$$

Now, for $f\in AC[a,b]$ we define the fractional derivatives of Riemann-Liouville and Caputo of order $\alpha$, respectively, by

$${}_a{}^{RL}D^{\alpha }f\left(z\right)={\displaystyle \frac{d}{dz}}{}_{a}I^{1-\alpha }f\left(z\right)={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\displaystyle \frac{d}{dz}}{\int }_a^{z}f\left(w\right){(z-w)}^{-\alpha }\mathrm{d}w$$

and

$${}_a{}^{C}D^{\alpha }f\left(z\right)={}_{a}I^{1-\alpha }{f}^{\prime }\left(z\right)={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_a^z{f}^{\prime }\left(w\right){(z-w)}^{-\alpha }\mathrm{d}w,$$

a.e. in $(a,b)$. The possibility to define these functions at $z=a$ must be analyzed in every case.

With these fractional derivative definitions in hand we can provide non-local definitions of the flux. In particular we note that (9) can be expressed in terms of the Caputo derivative

$$q_{\alpha }^{nl}(x,t)=-k_{\alpha }{}_0{}^{C}D_x^{\alpha }u(x,t),$$

(11)

referred to as the Caputo flux. In this way, replacing (11) in the balance equations (6) and (7) we obtain the governing equations of the space fractional Stefan problems considered in this article

$$\rho c{\displaystyle \frac{\partial u}{\partial t}}(x,t)-k_{\alpha }{\displaystyle \frac{\partial }{\partial x}}{}_0{}^{C}D_x^{\alpha }u(x,t)=0$$

(12)

and

$$\rho \ell \dot{s}\left(t\right)=k_{\alpha }\underset{x\to s{\left(t\right)}^{-}}{lim}{}_0{}^{C}D_x^{\alpha }u(x,t).$$

(13)

Here we have used (5) to compute

$$〚\rho H{〛}_l^s=\rho 〚H{〛}_l^s=\rho \ell .$$

Also, we have assumed that the non-local flux at the solid phase is given in terms of the local fluxes between the current position and the right extreme $x=l$, that is by

$$q_{\alpha }^{nl}(x,t)=-{\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_x^l\hspace{0.17em}{k}_{\alpha }{u}_x(p,t){(x-p)}^{-\alpha }\mathrm{dp}={}_x{}^{C}D_b^{\alpha }u(x,t),$$

which, according to (5), yields that $q_{\alpha }^{nl}(x,t)=0$ for every $x>s\left(t\right)$. Then

$$〚q{〛}_l^s=\hspace{0.17em}\underset{x\to s{\left(t\right)}^{-}}{lim}{q}_{\alpha }^{nl}(x,t)-\underset{x\to s{\left(t\right)}^{+}}{lim}{q}_{\alpha }^{nl}(x,t)=-k_{\alpha }\underset{x\to s{\left(t\right)}^{-}}{lim}{}_0{}^{C}D_x^{\alpha }u(x,t).$$

Thus, on defining the region $Q_s:=\{(x,t):0<x<s\left(t\right),0<t\}$ and providing appropriate initial and boundary conditions we obtain the Space Fractional Stefan Problem: Find the pair of functions $u:Q_s\to \mathbb{R}$ and $s:{\mathbb{R}}_0^{+}\to \mathbb{R}$ with enough regularity such that

$$\begin{array}{ccc}\left(i\right)\hfill & \rho c{u}_t(x,t)=k_{\alpha }{\displaystyle \frac{\partial }{\partial x}}{}_0{}^{C}D_x^{\alpha }u(x,t)\hfill & 0<x<s\left(t\right),0<t,\hfill \\ \left(ii\right)\hfill & u(0,t)=U_0>U_m\hfill & 0<t,\hfill \\ \left(iii\right)\hfill & u(s\left(t\right),t)=U_m\hfill & 0<t,\hfill \\ \left(iv\right)\hfill & s\left(0\right)=0,\hfill & \\ \left(v\right)\hfill & \rho \ell \dot{s}\left(t\right)=-k_{\alpha }\underset{x\to s{\left(t\right)}^{-}}{lim}{}_0{}^{C}D_x^{\alpha }u(x,t)\hfill & 0<t.\hfill \end{array}$$

(14)

We move now to defining a time fractional model, achieved through introducing the concept of memory enthalpy in the balance equations.

$$H^{\alpha }(x,t):={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_0^t{\eta }_{\alpha }[H(x,\tau )-H_0\left(x\right)]{(t-\tau )}^{-\alpha }\mathrm{d}\tau ={\eta }_{\alpha }{}_{0}I^{1-\alpha }[H(x,·)-H_0\left(x\right)]\left(t\right)$$

(15)

where $H_0\left(x\right)\equiv c{U}_m$ and ${\eta }_{\alpha }$ is a parameter that has been added to preserve dimensional consistency. It is important to stress that expression (15) corresponds to a continuous function that converges pointwise1 to $H-H_0$ when $\alpha \to 1^{+}$, and thus it can be interpreted as a regularization of the piecewise function $H-H_0$ which presents a finite jump at the interface. The idea of considering a regularized enthalpy is a common idea in the literature [3,27,30], the novelty of the approach here is to achieve this through a fractional integral as opposed to the more standard approach of allowing the phase change to occur smoothly over a small temperature region–a mushy zone.

On joining (4), (5) and (15) we have that

$$H^{\alpha }(x,t)=\left\{\begin{array}{cc}{\eta }_{\alpha }{}_{h\left(x\right)}I^{1-\alpha }[H(x,·)-H_0\left(x\right)]\left(t\right),\hfill & 0<x<s\left(t\right),\hfill \\ 0\hfill & x\ge s\left(t\right).\hfill \end{array}\right.$$

(16)

Then, by replacing (16) and (8) in the balance equation (6) we obtain the governing equation in the liquid region

$$\rho {\eta }_{\alpha }{\displaystyle \frac{\partial }{\partial t}}\left[{}_{h\left(x\right)}I_t^{1-\alpha }(H(x,t)-H_0\left(x\right))\right]=k{\displaystyle \frac{{\partial }^2}{{\partial }^{2}x}}u(x,t)$$

(17)

or equivalently

$$\rho {\eta }_{\alpha }{}_{h\left(x\right)}{}^{RL}D_t^{\alpha }[cu(x,t)+\ell -c{U}_m]=k{\displaystyle \frac{{\partial }^2}{{\partial }^{2}x}}u(x,t)$$

(18)

For the condition at the interface, we use the continuity of the memory enthalpy to deduce that

$$〚\rho {\eta }_{\alpha }{H}^{\alpha }{〛}_l^s=\rho {\eta }_{\alpha }\left[{\displaystyle \underset{x\to s{\left(t\right)}^{-}}{lim}{H}^{\alpha }(x,t)-\underset{x\to s{\left(t\right)}^{+}}{lim}{H}^{\alpha }(x,t)}\right]=0,$$

(19)

Then the classical Stefan condition at the interface (1)$-\left(v\right)$ is “lost” in the memory enthalpy model and it can be replaced by

$$u_x(s{\left(t\right)}^{-},t)=0.$$

(20)

Remark 1. 

Condition (20) was previously obtained for the memory flux Stefan problem in [14] and later in [21]. It is a natural condition for the memory enthalpy problem, which comes from the continuity of the “memory enthalpy” at the interface.

Following from the above analysis, we can present the Time Fractional Stefan Problem: Find the pair of functions $u:Q_s\to \mathbb{R}$ and $s:{\mathbb{R}}_0^{+}\to \mathbb{R}$ with enough regularity such that

$$\begin{array}{ccc}\left(i\right)\hfill & \rho {\eta }_{\alpha }{}_{h\left(x\right)}{}^{RL}D_t^{\alpha }\left[cu(x,·)+\ell -c{U}_m\right]\left(t\right)=k{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}(x,t),\hfill & 0<x<s\left(t\right),\hspace{0.17em}0<t,\hspace{0.17em}\hspace{0.17em}\hfill \\ \left(ii\right)\hfill & u(0,t)=U_0,\hfill & 0<t,\hfill \\ \left(iii\right)\hfill & \hspace{0.17em}u(s\left(t\right),t)=U_m,\hfill & x\ge s\left(t\right),\hfill \\ \left(iv\right)\hfill & s\left(0\right)=0,\hfill & \\ \left(v\right)\hfill & u_x(s{\left(t\right)}^{-},t)=0,\hfill & t>0.\hfill \end{array}$$

(21)

3. The Dimensionless Problems

In this subsection, we aim to rewrite problems (14) and (21) in a dimensionless form, allowing us to recover the quasi-stationary problems associated with each case in the limit as the dimensionless Stefan number approaches zero. Note that the dimensionless form of the time fractional problem has previously been derived in [19][Prop. 8], the addition of the spacial fractional case is new.

To start our derivations of the dimensionless problem formulations we identify the dimensions of the problem properties in Table 1 and provide a Proposition 1 to determine the dimensions of fractional derivative operators.

Proposition 1. 

For every $\alpha \in (0,1)$ it holds that: 

(1)

$\left[{}^{C}D_x^{\alpha }f\right]={\displaystyle \frac{\left[f\right]}{{\left[\mathrm{L}\right]}^{\alpha }}}$.

(2)

$\left[{\displaystyle \frac{\partial }{\partial x}}{}^{C}D_x^{\alpha }f\right]={\displaystyle \frac{\left[f\right]}{{\left[\mathrm{L}\right]}^{1+\alpha }}}$.

(3)

$\left[{}_a{}^{RL}D_t^{\alpha }f\right]={\displaystyle \frac{\left[f\right]}{{\left[\mathrm{T}\right]}^{\alpha }}}$.

Proof. 

We only give the proof of 1, which is a simple compute.

$${[}^C{D}_x^{\alpha }f]=\left[{\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_0^x{\displaystyle \frac{f^{\prime }}{{(x-s)}^{\alpha }}}ds\right]=\left[f^{\prime }\right]{\left[{(x-·)}^{\alpha }\right]}^{-1}\left[ds\right]=\left[f\right]{\left[{(x-·)}^{\alpha }\right]}^{-1}={\displaystyle \frac{\left[f\right]}{{\left[\mathrm{L}\right]}^{\alpha }}}.$$

□

Remark 2. 

At this point it is important to stress that the classical change of variable for the time, which is generally defined by $x={\displaystyle \frac{y}{x_0}}$ and $\tau ={\displaystyle \frac{{\lambda }^2}{x_0^2}}t$, is not used in this article because we are specially interested in obtaining a non-dimensional model that permits to obtain the quasi-stationary case when the dimensionless Stefan number (see Table 1) approaches zero.

The importance given to the quasi-steady-state case lies in the fact that we exactly know the value of the prefactor of the free boundary for any given fractional order $0<\alpha \le 1$.

3.1. The Dimensionless Spacial Fractional Stefan Problem

To simplify, without loss of generality, we rewrite problem (14) with the melt temperature $U_m=0$.

$$\begin{array}{ccc}\left(i\right)\hfill & \rho c{u}_t(x,t)=k{\nu }_{\alpha }{\displaystyle \frac{\partial }{\partial x}}{}_0{}^{C}D_x^{\alpha }u(x,t),\hfill & 0<x<s\left(t\right),0<t,\hfill \\ \left(ii\right)\hfill & u(0,t)=U_0,\hfill & 0<t,\hfill \\ \left(iii\right)\hfill & u\left(s\right(t),t)=0,\hfill & 0<t,\hfill \\ \left(iv\right)\hfill & s\left(0\right)=0,\hfill \\ \left(v\right)\hfill & \rho \ell s^{\prime }\left(t\right)=-k{\nu }_{\alpha }{}_0{}^{C}D_x^{\alpha }u(s{\left(t\right)}^{-},t),\hfill & 0<x<s\left(t\right).\hfill \end{array}$$

(22)

First, let us clarify the units of measure of the parameter ${\nu }_{\alpha }$ which was added to give physical dimension consistency. The appropriate dimensions follow directly from the governing equation (22)-(i). By applying Proposition 1 and the dimensions noted in Table 1 to this equation we arrive at the following dimensional balance

$${\displaystyle \frac{\left[\mathrm{M}\right]}{\left[{\mathrm{L}}^3\right]}}·{\displaystyle \frac{\left[{\mathrm{L}}^2\right]}{\left[{\mathrm{MT}}^2\right]}}·{\displaystyle \frac{\left[\mathrm{K}\right]}{\left[\mathrm{T}\right]}}={\displaystyle \frac{\left[\mathrm{ML}\right]}{\left[{\mathrm{KT}}^3\right]}}·\left[{\nu }_{\alpha }\right]·{\displaystyle \frac{\left[\mathrm{K}\right]}{\left[{\mathrm{L}}^{1+\alpha }\right]}},$$

(23)

which implies that

$$\left[{\nu }_{\alpha }\right]=\left[{\mathrm{L}}^{\alpha -1}\right]\mathrm{and}\underset{\alpha \to 1^{-}}{lim}{\nu }_{\alpha }=[-].$$

(24)

Change of variables for the Space Fractional case. If we let $x_0$ be a characteristic position we can define the the following dimensionless space and time variables in terms of the quantities in Table 1

$$y={\displaystyle \frac{x}{x_0}},\tau =Ste\hspace{0.17em}{\displaystyle \frac{{\lambda }^2}{x_0^2}}t.$$

(25)

In addition, on using $U_0$ as the characteristic temperature, we can construct dimensionless dependent variables for temperature and front position as

$$\tilde{u}(y,\tau )={\displaystyle \frac{1}{U_0}}u(x,t)\mathrm{and}\tilde{s}\left(\tau \right)={\displaystyle \frac{s\left(t\right)}{x_0}}.$$

(26)

Then, by making the substitution $\xi =\nu x_0$ in the fractional integral, we have

$$\begin{array}{cc}\hfill {}_0{}^{C}D_x^{\alpha }u(x,t)& ={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_0^x{\displaystyle \frac{u_x(\xi ,t)}{{(x-\xi )}^{\alpha }}}d\xi ={\displaystyle \frac{U_0}{\mathrm{\Gamma }(1-\alpha )}}{\int }_0^x{\displaystyle \frac{{\tilde{u}}_y(y\left(\xi \right),\tau )}{x_0{(x-\xi )}^{\alpha }}}d\xi \hfill \\ \hfill & ={\displaystyle \frac{U_0}{\mathrm{\Gamma }(1-\alpha )}}{\int }_0^y{\displaystyle \frac{{\tilde{u}}_y(\nu ,\tau )}{x_0^{\alpha }{(y-\nu )}^{\alpha }}}d\nu ={\displaystyle \frac{U_0}{x_0^{\alpha }}}{}_0{}^{C}D_y^{\alpha }\tilde{u}(y,\tau ).\hfill \end{array}$$

(27)

$${\displaystyle \frac{\partial }{\partial x}}{}_0{}^{C}D_x^{\alpha }u(x,t)={\displaystyle \frac{U_0}{x_0^{\alpha }}}{\displaystyle \frac{\partial }{\partial y}}{}_0{}^{C}D_y^{\alpha }\tilde{u}(y,\tau )·{\displaystyle \frac{\partial y}{\partial x}}={\displaystyle \frac{U_0}{x_0^{1+\alpha }}}{\displaystyle \frac{\partial }{\partial y}}{}_0{}^{C}D_y^{\alpha }\tilde{u}(y,\tau ).$$

(28)

And

$$u_t(x,t)={\tilde{u}}_{\tau }(y,\tau ){\displaystyle \frac{\partial \tau }{\partial t}}=Ste\hspace{0.17em}{U}_0{\displaystyle \frac{{\lambda }^2}{x_0^2}}{\tilde{u}}_{\tau }(y,\tau ).$$

(29)

From (22)$-\left(i\right)$, (28) and (29) we obtain the dimensionless governing equation for the liquid phase

$$Ste\hspace{0.17em}{x}_0^{\alpha -1}{\nu }_{\alpha }^{-1}{\tilde{u}}_{\tau }(y,\tau )={\displaystyle \frac{\partial }{\partial y}}{}_0{}^{C}D^{\alpha }\tilde{u}(y,\tau ).$$

(30)

From (22)$-\left(v\right)$, (26) and (27) we have

$${\tilde{s}}^{\prime }\left(\tau \right)=-x_0^{1-\alpha }{\nu }_{\alpha }{}_0{}^{C}D_y^{\alpha }\tilde{u}(\tilde{s}{\left(\tau \right)}^{-},\tau ).$$

(31)

To retain the dimensionless nature of the equations we set the admissible parameter as

$${\nu }_{\alpha }:=x_0^{\alpha -1},$$

(32)

which has dimensions consistent with (24).

Finally we can rewrite (30) and (31) as

$$Ste\hspace{0.17em}{\tilde{u}}_{\tau }(y,\tau )={\displaystyle \frac{\partial }{\partial y}}{}_0{}^{C}D_y^{\alpha }\tilde{u}(y,\tau ).$$

(33)

$${\tilde{s}}^{\prime }\left(\tau \right)=-{}_0{}^{C}D_y^{\alpha }\tilde{u}(\tilde{s}{\left(\tau \right)}^{-},\tau ).$$

(34)

Adding appropriate initial and boundary conditions we arrive at the desired dimensionless space fractional Stefan problem: Find the pair of functions $\tilde{u}:Q_{\tilde{s}}\to \mathbb{R}$ and $\tilde{s}:{\mathbb{R}}_0^{+}\to \mathbb{R}$ with enough regularity such that

$$\begin{array}{ccc}\left(i\right)\hfill & Ste\hspace{0.17em}{\tilde{u}}_{\tau }(y,\tau )={\displaystyle \frac{\partial }{\partial y}}{}_0{}^{C}D^{\alpha }\tilde{u}(y,\tau ),\hfill & 0<y<\tilde{s}\left(\tau \right),0<\tau ,\hfill \\ \left(ii\right)\hfill & \tilde{u}(0,\tau )=1,\hfill & 0<\tau ,\hfill \\ \left(iii\right)\hfill & \tilde{u}(\tilde{s}\left(\tau \right),\tau )=0,\hfill & 0<\tau ,\hfill \\ \left(iv\right)\hfill & \tilde{s}\left(0\right)=0,\hfill \\ \left(v\right)\hfill & {\tilde{s}}^{\prime }\left(\tau \right)=-{}_0{}^{C}D^{\alpha }\tilde{u}(\tilde{s}{\left(\tau \right)}^{-},\tau ),\hfill & 0<\tau .\hfill \end{array}$$

(35)

3.2. The Dimensionless Time Fractional Stefan Problem

Let us work now with problem (21) for the case $U_m=0$. That is, consider the problem to find the pair of functions $u:Q_s\to \mathbb{R}$ and $s:{\mathbb{R}}_0^{+}\to \mathbb{R}$ with enough regularity such that

$$\begin{array}{ccc}\left(i\right)\hfill & \rho {\eta }_{\alpha }{}_{h\left(x\right)}{}^{RL}D_t^{\alpha }\left[cu(x,·)+\ell \right]\left(t\right)=k{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}(x,t),\hfill & 0<x<s\left(t\right),\hspace{0.17em}0<t<T,\hspace{0.17em}\hspace{0.17em}\hfill \\ \left(ii\right)\hfill & u(0,t)=U_0,\hfill & 0<t\le T,\hfill \\ \left(iii\right)\hfill & u\left(s\right(t),t)=0,\hfill & 0<t\le T,\hfill \\ \left(iv\right)\hfill & s\left(0\right)=0,\hfill & \\ \left(v\right)\hfill & u_x(s{\left(t\right)}^{-},t)=0,\hfill & 0<t\le T.\hfill \end{array}$$

(36)

As before, we start by clarifying the role of the parameter ${\eta }_{\alpha }$. But first, it is important to introduce an equivalent problem to (36) that will be used in this subsection. More precisely, according to [21, Prop. 5], problem (36) is equivalent to the problem to find the pair of functions $u:Q_s\to \mathbb{R}$ and $s:{\mathbb{R}}_0^{+}\to \mathbb{R}$ such that

$$\begin{array}{ccc}\left(i\right)\hfill & {\eta }_{\alpha }\rho c{}_{h\left(x\right)}{}^{C}D_t^{\alpha }u(x,t)+{\displaystyle \frac{{\eta }_{\alpha }\rho \ell }{\mathrm{\Gamma }(1-\alpha ){(t-h\left(x\right))}^{\alpha }}}=k{\displaystyle \frac{{\partial }^{2}u}{\partial x^2}}(x,t),\hfill & 0<x<s\left(t\right),\hspace{0.17em}0<t<T,\hspace{0.17em}\hspace{0.17em}\hfill \\ \left(ii\right)\hfill & u(0,t)=U_0,\hfill & 0<t\le T,\hfill \\ \left(iii\right)\hfill & u\left(s\right(t),t)=0,\hfill & 0<t\le T,\hfill \\ \left(iv\right)\hfill & s\left(0\right)=0,\hfill & \\ \left(v\right)\hfill & u_x(s{\left(t\right)}^{-},t)=0,\hfill & 0<t\le T.\hfill \end{array}$$

(37)

Note that we can rewrite equation (37)$-\left(i\right)$ by using the diffusion coefficient ${\lambda }^2$ and defining ${\mu }_{\alpha }:={\left({\eta }_{\alpha }\right)}^{-1}$ in the following way.

$${}_{h\left(x\right)}{}^{C}D_t^{\alpha }u(x,t)+{\displaystyle \frac{\ell }{c\mathrm{\Gamma }(1-\alpha ){(t-h\left(x\right))}^{\alpha }}}={\lambda }^2{\mu }_{\alpha }{u}_{xx}(x,t).$$

(38)

From this equation, using Proposition 1, Table 1, we obtain the dimensional balance

$${\displaystyle \frac{\left[\mathrm{K}\right]}{\left[{\mathrm{T}}^{\alpha }\right]}}={\displaystyle \frac{\left[{\mathrm{L}}^2\right]}{\left[\mathrm{T}\right]}}\left[{\mu }_{\alpha }\right]{\displaystyle \frac{\left[\mathrm{K}\right]}{\left[{\mathrm{L}}^2\right]}}$$

(39)

From which it follows that

$$\left[{\mu }_{\alpha }\right]=\left[T^{1-\alpha }\right]\mathrm{and}\underset{\alpha \to 1^{-}}{lim}{\mu }_{\alpha }=1.$$

(40)

Remark 3. 

Consistent with Remark 2, we consider the non-classical change of variable (25) in order to achieve the quasi-stationary problem when the Stefan number approaches zero.

Change of variables for the Time Fractional case. Let $x_0$ be a characteristic position and consider the change of variables (25). We define as before,

$$\tilde{u}(y,\tau )={\displaystyle \frac{1}{U_0}}u(x,t),\tilde{s}\left(\tau \right)={\displaystyle \frac{s\left(t\right)}{x_0}},\tilde{h}\left(y\right)=Ste\hspace{0.17em}{\displaystyle \frac{{\lambda }^2}{x_0^2}}h\left(x\right).$$

(41)

where $x,y,t$ and $\tau$ are related by (25). Applying the substitution $r={\displaystyle \frac{x_0^2}{Ste\hspace{0.17em}{\lambda }^2}}v$ in the following integral, we deduce that

$$\begin{array}{cc}\hfill {}_{h\left(x\right)}{}^{C}D_t^{\alpha }u(x,t)& ={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{h\left(x\right)}^t{\displaystyle \frac{u_r(x,r)}{{(t-r)}^{\alpha }}}dr={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{h\left(x\right)}^t{\displaystyle \frac{U_0{\tilde{u}}_{\tau }(y,\tau \left(r\right))}{{(t-r)}^{\alpha }}}Ste\hspace{0.17em}{\displaystyle \frac{{\lambda }^2}{x_0^2}}dr\hfill \\ \hfill & ={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\tilde{h}\left(y\right)}^{\tau }{\displaystyle \frac{U_{0}Ste\hspace{0.17em}\frac{{\lambda }^2}{x_0^2}{\tilde{u}}_{\tau }(y,v)}{{\left(t-\frac{x_0^2}{Ste\hspace{0.17em}{\lambda }^2}v\right)}^{\alpha }}}{\displaystyle \frac{x_0^2}{Ste\hspace{0.17em}{\lambda }^2}}dv={\displaystyle \frac{{\lambda }^{2\alpha }}{x_0^{2\alpha }}}St{e}^{\alpha }{U0}_{\tilde{h}\left(y\right)}^C{D}_{\tau }^{\alpha }\tilde{u}(y,\tau ).\hfill \end{array}$$

(42)

Also,

$$t=h\left(x\right)\iff {\displaystyle \frac{x_0^2}{Ste\hspace{0.17em}{\lambda }^2}}\tau =h\left(x\right)\iff \tau =\tilde{h}\left(y\right),$$

$$t=h\left(x\right)\iff s\left(t\right)=x\iff s\left(t\right)=x_{0}y\iff \tilde{s}\left(\tau \right)=y.$$

Then

$${\displaystyle \frac{\ell }{c\mathrm{\Gamma }(1-\alpha ){(t-h\left(x\right))}^{\alpha }}}={\displaystyle \frac{\ell }{c\mathrm{\Gamma }(1-\alpha ){\left(\frac{x_0^2}{Ste\hspace{0.17em}{\lambda }^2}\tau -\frac{x_0^2}{Ste\hspace{0.17em}{\lambda }^2}\tilde{h}\left(y\right)\right)}^{\alpha }}}={\displaystyle \frac{{\lambda }^{2\alpha }}{x_0^{2\alpha }}}{\displaystyle \frac{St{e}^{\alpha }\ell }{c\mathrm{\Gamma }(1-\alpha ){\left(\tau -\tilde{h}\left(y\right)\right)}^{\alpha }}}.$$

(43)

From the other side,

$${\lambda }^2{\mu }_{\alpha }{u}_{xx}(x,t)={\lambda }^2{\mu }_{\alpha }{\displaystyle \frac{U_0}{x_0^2}}{\tilde{u}}_{yy}(y,\tau ).$$

(44)

Substituting (42), (43) and (44) in (38) gives

$${\displaystyle \frac{{\lambda }^{2\alpha }}{x_0^{2\alpha }}}St{e}^{\alpha }{U}_0{}_{\tilde{h}\left(y\right)}{}^{C}D_{\tau }^{\alpha }\tilde{u}(y,\tau )+{\displaystyle \frac{{\lambda }^{2\alpha }}{x_0^{2\alpha }}}{\displaystyle \frac{St{e}^{\alpha }\ell }{c\mathrm{\Gamma }(1-\alpha ){\left(\tau -\tilde{h}\left(y\right)\right)}^{\alpha }}}={\lambda }^2{\mu }_{\alpha }{\displaystyle \frac{U_0}{x_0^2}}{\tilde{u}}_{yy}(y,\tau ).$$

(45)

Note that ${\mu }_{\alpha }:={\left({\displaystyle \frac{x_0^2}{Ste\hspace{0.17em}{\lambda }^2}}\right)}^{1-\alpha }$ is an admissible parameter verifying (40). Then by replacing this parameter into (45) and addressing with initial and boundary conditions, the following dimensionless problem associated to the Time Fractional Stefan Problem is obtained: Find the pair of functions $\tilde{u}:Q_{\tilde{s}}\to \mathbb{R}$ and $\tilde{s}:{\mathbb{R}}_0^{+}\to \mathbb{R}$ enough regular such that

$$\begin{array}{ccc}\left(i\right)\hfill & {Ste}_{\tilde{h}\left(y\right)}^C{D}_{\tau }^{\alpha }\tilde{u}(y,\tau )+{\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha ){\left(\tau -\tilde{h}\left(y\right)\right)}^{\alpha }}}={\tilde{u}}_{yy}(y,\tau ),\hfill & 0<y<\tilde{s}\left(\tau \right),\hspace{0.17em}0<\tau ,\hfill \\ \left(ii\right)\hfill & \tilde{u}(0,\tau )=1,\hfill & 0<\tau ,\hfill \\ \left(iii\right)\hfill & \tilde{u}(\tilde{s}\left(\tau \right),\tau )=0,\hfill & 0<\tau ,\hfill \\ \left(iv\right)\hfill & \tilde{s}\left(0\right)=0,\hfill & \\ \left(v\right)\hfill & {\tilde{u}}_y(\tilde{s}{\left(\tau \right)}^{-},\tau )=0,\hfill & 0<\tau .\hfill \end{array}$$

(46)

4. Computing the Prefactor for the Space Fractional Case

4.1. Closed Solutions

Let us start with the quasi-stationary case. This problem is usually attained by making the specific heat $c\to 0$. The resulting problem is trivial in (14) or (21) but we focus now in problem (35). Thus, by assuming that the latent heat ℓ and the imposed temperature $U_0$ are finite positive numbers, it is immediately clear from the Stefan number definition in Table 1 that

$$c\to 0\iff Ste\to 0.$$

(47)

Then, the limit problem obtained from (35) by making $c\to 0$ is the Space Fractional quasi-Stationary Problem: Find the pair of functions $\tilde{u}:Q_{\tilde{s}}\to \mathbb{R}$ and $\tilde{s}:{\mathbb{R}}_0^{+}\to \mathbb{R}$ with enough regularity such that

$$\begin{array}{ccc}\left(i\right)\hfill & {\displaystyle \frac{\partial }{\partial y}}{}_0{}^{C}D^{\alpha }\tilde{u}(y,\tau )=0,\hfill & 0<y<\tilde{s}\left(\tau \right),0<\tau ,\hfill \\ \left(ii\right)\hfill & \tilde{u}(0,\tau )=1,\hfill & 0<\tau ,\hfill \\ \left(iii\right)\hfill & \tilde{u}(\tilde{s}\left(\tau \right),\tau )=0,\hfill & 0<\tau ,\hfill \\ \left(iv\right)\hfill & \tilde{s}\left(0\right)=0,\hfill \\ \left(v\right)\hfill & {\tilde{s}}^{\prime }\left(\tau \right)=-{}_0{}^{C}D^{\alpha }\tilde{u}(\tilde{s}{\left(\tau \right)}^{-},\tau ),\hfill & 0<\tau .\hfill \end{array}$$

(48)

Remark 4. 

The change of variable performed in (25)-(26) and (41) remains unaffected when considering the limit $c\to 0$. Indeed,

$$Ste{\lambda }^2={\displaystyle \frac{U_{0}k}{\rho \ell }}$$

is independent of c.

Problem (48) was solved in [28] and its close solution is given by the pair

$${\tilde{u}}_{\alpha ,0}(y,\tau )=1-{\displaystyle \frac{y^{\alpha }}{{\left[\mathrm{\Gamma }(2+\alpha )\right]}^{\frac{\alpha }{1+\alpha }}{\tau }^{\frac{\alpha }{1+\alpha }}}},{\tilde{s}}_{\alpha ,0}\left(\tau \right)={\left[\mathrm{\Gamma }(2+\alpha )\right]}^{{\displaystyle \frac{1}{1+\alpha }}}{\tau }^{{\displaystyle \frac{1}{1+\alpha }}}.$$

(49)

We now present an exact solution to problem (35), originally introduced in [20], which has been adapted here for the dimensionless case with the Stefan number as the main parameter. To that end, let us introduce the special functions involved

Definition 1. 

Let $\alpha >0,m>0$, and l such that $\alpha (jm+l)\ne -1,-2,-3,\cdots \hspace{0.17em}(j=0,1,2,\cdots )$. The three-parametric Mittag-Leffler function $E_{\alpha ,m,l}\left(z\right)$ is defined by

$$E_{\alpha ,m,l}\left(z\right)=\sum _{n=0}^{\infty }{c}_n{z}^n,\hspace{0.17em}\mathrm{with}\hspace{0.17em}{c}_0=1,\hspace{0.17em}\hspace{0.17em}{c}_n=\prod _{j=0}^{n-1}{\displaystyle \frac{\mathrm{\Gamma }\left(\alpha \right(jm+l)+1)}{\mathrm{\Gamma }\left(\alpha \right(jm+l+1)+1)}},\hspace{0.17em}\hspace{0.17em}(n=1,2,3,\cdots ).$$

(50)

Remark 5. 

A particular case is $E_{1,1,0}\left(z\right)=e^z$ and we recover the classical Mittag-Leffler function for $m=1$ and $l=0$$E_{\alpha ,1,0}\left(z\right)=E_{\alpha }\left(z\right)$. Also, a two parametric Mittag–Leffler function is recovered for the case $E_{\alpha ,1,l}\left(z\right)=\mathrm{\Gamma }(\alpha l+1)E_{\alpha ,\alpha l+1}\left(z\right)$ and the special case of our interest is $E_{1,2,1}\left(-{\displaystyle \frac{z^2}{2}}\right)=e^{-{\left({\displaystyle \frac{z}{2}}\right)}^2}.$ For the role of the three-parametric Mittag-Leffler function we refer the reader to the original works of Kilbas and Siago

We will consider the parameters $l=1,m=1+{\displaystyle \frac{1}{\alpha }}$. For this case we have

$$E_{\alpha ,1+{\displaystyle \frac{1}{\alpha }},1}\left(z\right)=\sum _{n=0}^{\infty }{c}_n{z}^n,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{with}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{c}_0=1,\hspace{0.17em}\hspace{0.17em}{c}_n=\prod _{j=0}^{n-1}{\displaystyle \frac{\mathrm{\Gamma }\left(\alpha \right(j(1+1/\alpha )+1)+1)}{\mathrm{\Gamma }\left(\alpha \right(j(1+1/\alpha )+2)+1)}},\hspace{0.17em}\hspace{0.17em}n\in \mathbb{N}$$

(51)

and the coefficients $c_n$ verify the next recursive form

$$c_0=1,\hspace{0.17em}\hspace{0.17em}{c}_n=c_{n-1}{\displaystyle \frac{\mathrm{\Gamma }\left(\right(n-1\left)\right(\alpha +1)+\alpha +1)}{\mathrm{\Gamma }\left(\right(n-1\left)\right(\alpha +1)+2\alpha +1)}}=c_{n-1}{\displaystyle \frac{\mathrm{\Gamma }\left(n\right(1+\alpha \left)\right)}{\mathrm{\Gamma }\left(n\right(1+\alpha )+\alpha )}},\hspace{0.17em}\hspace{0.17em}n\in \mathbb{N}.$$

(52)

We know from [12,Th.1] that the three-parametric Mittag-Leffler functions (50) already defined are entire functions for every $\alpha >0$.

Define the function $\sigma :{\mathbb{R}}_0^{+}\to \mathbb{R}$ by

$${\sigma }_{\alpha }\left(w\right):=w^{\alpha -1}{E}_{\alpha ,1+{\displaystyle \frac{1}{\alpha }},1}\left(-{\displaystyle \frac{w^{1+\alpha }}{1+\alpha }}\right)=\sum _{n=0}^{\infty }{c}_n{(-1)}^n{\displaystyle \frac{w^{(n+1)(1+\alpha )-2}}{{(1+\alpha )}^n}}.$$

(53)

Now, by mimicking the steps in [20, Section5], it is straightforward that the pair $\{{\tilde{u}}_{\alpha },{\tilde{s}}_{\alpha }\}$ given by

$${\tilde{u}}_{\alpha ,Ste}(y,\tau )=1-{\displaystyle \frac{1}{{\displaystyle {\int }_0^{{\xi }_{\alpha ,Ste}}{\sigma }_{\alpha }\left(w\right)\mathrm{d}w}}}{\int }_0^{y/{(\tau /Ste)}^{1/(1+\alpha )}}{\sigma }_{\alpha }\left(w\right)\mathrm{d}w$$

(54)

and

$${\tilde{s}}_{\alpha ,Ste}\left(\tau \right)={\xi }_{\alpha ,Ste}{(\tau /Ste)}^{{\displaystyle \frac{1}{1+\alpha }}},\tau \hspace{0.17em}\in \hspace{0.17em}(0,T),$$

(55)

where ${\xi }_{\alpha ,Ste}\in \hspace{0.17em}{\mathbb{R}}^{+}$ is the unique solution to the equation

$${\displaystyle Ste\left[{\displaystyle \mathrm{\Gamma }\left(\alpha \right)(1+\alpha )-{\int }_0^{y}w{\sigma }_{\alpha }\left(w\right)\mathrm{d}w}\right]=y{\int }_0^y{\sigma }_{\alpha }\left(w\right)\mathrm{d}w,}$$

(56)

is the solution to problem (35).

Remark 6. 

Let us highlight that the expression for theprefactor of the free boundarythat we are trying to approximate is given by

$${\theta }_{\alpha ,Ste}={\displaystyle \frac{{\xi }_{\alpha ,Ste}}{St{e}^{\frac{1}{1+\alpha }}}}.$$

(57)

The results in the next proposition can be founded in [20, Prop.8 and Section5].

Proposition 2. 

The function ${\sigma }_{\alpha }:{\mathbb{R}}^{+}\to \mathbb{R}$ given in (53) verifies the following properties: 

1.

${\sigma }_{\alpha }$ is a non-negative function such that ${\sigma }_{\alpha }\neg \equiv 0$.

2.

The following limits hold:

$$\underset{y\searrow 0}{lim}{\int }_0^{y}w{\sigma }_{\alpha }\left(w\right)dw=0^{+},\underset{y\searrow 0}{lim}{\int }_0^y{\sigma }_{\alpha }\left(w\right)dw=0^{+}.$$

(58)

3.

The function $H:{\mathbb{R}}_0^{+}\to \mathbb{R}$ given by ${\displaystyle H\left(y\right)=\frac{\mathrm{\Gamma }\left(\alpha \right)(1+\alpha )-{\int }_0^{y}w{\sigma }_{\alpha }\left(w\right)dw}{{\int }_0^y{\sigma }_{\alpha }\left(w\right)dw}}$ is a decreasing function.

4.2. Computing ${\theta }_{\alpha ,Ste}$

Using series expansion of $E_{\alpha ,m,l}\left(z\right)$ in the definition of ${\sigma }_{\alpha }\left(\omega \right)$ we can rewrite the equation (56) as

$$\sum _{n=0}^{\infty }{\displaystyle \frac{c_n{(-1)}^n}{{(1+\alpha )}^n}}\left[{\displaystyle \frac{Ste}{(1+\alpha )(n+1)}}+{\displaystyle \frac{1}{(1+\alpha )(n+1)-1}}\right]x^{(n+1)(1+\alpha )}-Ste·\mathrm{\Gamma }\left(\alpha \right)(1+\alpha )=0.$$

(59)

The solution to the last equation (59) can be interpreted as the root of the following function in the left hand side. Thus, a simple method like the bisection method can be applied to find the solution ${\xi }_{\alpha }$ to (59), after making an efficient computation of the 3-parametric Mittag-Leffler functions (in [5], a simple code is presented to see how it has been computed).

4.3. Analysis of the Convergence to the Quasi-Stationary Case

As it was stated in Section 3, the quasi-stationary problem is obtained by making the sensible heat $c\to 0$ in (35). However taking into account (47) and being the Stefan number the visible parameter in our equations, the convergence to the quasi-stationary case will be analyzed by making $Ste\to 0$. It is also worth noting that the convergence analysis is not affected by the change of variable, according to Remark 4.

In view of the results obtained from the previous algorithm for ${\theta }_{\alpha ,Ste}$ for a fix $\alpha$ (see Table 2), the next theorem is presented.

Theorem 1. 

The prefactor of the free boundary (57) corresponding to the Space Fractional Stefan Problem (35) converges to the prefactor of the quasi-stationary Space Fractional Stefan problem (48) when the Stefan number approaches zero. In other words,

$$\underset{Ste\to 0}{lim}{\theta }_{\alpha ,Ste}=\mathrm{\Gamma }{(2+\alpha )}^{1/(1+\alpha )}.$$

(60)

Proof. 

For every $Ste>0$ let ${\xi }_{\alpha ,Ste}$ the unique positive solution to (56). Equivalently, ${\xi }_{\alpha ,Ste}$ verifies the equation

$$Ste\hspace{0.17em}H\left({\xi }_{\alpha ,Ste}\right)={\xi }_{\alpha ,Ste}$$

where H is the function defined in Proposition 2. Let $St{e}_1<St{e}_2$. If we suppose that ${\xi }_{\alpha ,St{e}_2}\le {\xi }_{\alpha ,St{e}_1}$, from Proposition 2 item 3 we deduce that $H\left({\xi }_{\alpha ,St{e}_2}\right)\ge H\left({\xi }_{\alpha ,St{e}_1}\right)$. Then, ${\xi }_{\alpha ,St{e}_2}=St{e}_{2}H\left({\xi }_{\alpha ,St{e}_2}\right)>St{e}_{1}H\left({\xi }_{\alpha ,St{e}_1}\right)={\xi }_{\alpha ,St{e}_1}$, which is a contradiction. Then the family ${\left\{{\xi }_{\alpha ,Ste}\right\}}_{Ste}$ is decreasing respect on the parameter $Ste$ and being the positiveness of the elements we deduce that ${\displaystyle \underset{Ste\to 0}{lim}{\xi }_{\alpha ,Ste}=\beta \ge 0.}$ Now, making $Ste\to 0$ in (56) yields that $0\le \beta {\int }_0^{\beta }{\sigma }_{\alpha }\left(w\right)\mathrm{d}w=0,$ and from Proposition 2 we have

$$\underset{Ste\to 0}{lim}{\xi }_{\alpha ,Ste}=0.$$

(61)

Now, using (56) again

$${\theta }_{\alpha ,Ste}^{1+\alpha }={\left({\displaystyle \frac{{\xi }_{\alpha ,Ste}}{St{e}^{1/(1+\alpha )}}}\right)}^{1+\alpha }={\displaystyle \frac{{\xi }_{\alpha ,Ste}^{\alpha }}{{\int }_0^{{\xi }_{\alpha ,Ste}}{\sigma }_{\alpha }\left(w\right)dw}}·\left({\displaystyle \mathrm{\Gamma }\left(\alpha \right)(1+\alpha )-{\int }_0^{{\xi }_{\alpha ,Ste}}w{\sigma }_{\alpha }\left(w\right)dw}\right).$$

(62)

From (61) and Proposition 2,

$${\displaystyle \underset{Ste\to 0}{lim}\mathrm{\Gamma }\left(\alpha \right)(1+\alpha )-{\int }_0^{{\xi }_{\alpha ,Ste}}w{\sigma }_{\alpha }\left(w\right)dw=\mathrm{\Gamma }\left(\alpha \right)(1+\alpha ).}$$

(63)

By using the series expansion and (61) we compute the limit in the first factor in (62)

$$\underset{Ste\to 0}{lim}{\displaystyle \frac{{\int }_0^{{\xi }_{\alpha ,Ste}}{\sigma }_{\alpha }\left(w\right)dw}{{\xi }_{\alpha ,Ste}^{\alpha }}}=\underset{Ste\to 0}{lim}\left[{\displaystyle \frac{1}{\alpha }}+\sum _{n=1}^{\infty }{\displaystyle \frac{c_n{(-1)}^n}{{(1+\alpha )}^n}}{\displaystyle \frac{{\xi }_{\alpha ,Ste}^{n(1+\alpha )}}{(n+1)(1+\alpha )-1}}\right]={\displaystyle \frac{1}{\alpha }}.$$

(64)

Finally, we make $Ste\to 0$ in (62) and use (63), (64) and the property of the Gamma function saying that $z\mathrm{\Gamma }\left(z\right)=\mathrm{\Gamma }(z+1)$ to obtain

$$\underset{Ste\to 0}{lim}{\theta }_{\alpha ,Ste}=\underset{Ste\to 0}{lim}{\displaystyle \frac{{\xi }_{\alpha ,Ste}}{St{e}^{1/(1+\alpha )}}}={\left(\alpha \mathrm{\Gamma }\left(\alpha \right)(1+\alpha )\right)}^{1/(1+\alpha )}=\mathrm{\Gamma }{(2+\alpha )}^{1/(1+\alpha )}.$$

(65)

□

Corollary 1. 

The solution (54)-(55) to the Space Fractional Stefan Problem (35) converges to the solution (49) to the quasi-stationary Space Fractional Stefan problem (48) when the Stefan number approaches zero.

Proof. 

Going back to (55), we use (65) to claim that

$$\underset{Ste\to 0}{lim}{\tilde{s}}_{\alpha ,Ste}\left(\tau \right)=\underset{Ste\to 0}{lim}{\displaystyle \frac{{\xi }_{\alpha ,Ste}}{St{e}^{\frac{1}{1+\alpha }}}}{\tau }^{{\displaystyle \frac{1}{1+\alpha }}}=\mathrm{\Gamma }{(2+\alpha )}^{{\displaystyle \frac{1}{1+\alpha }}}{\tau }^{{\displaystyle \frac{1}{1+\alpha }}}={\tilde{s}}_{\alpha ,0}\left(\tau \right).$$

(66)

For the limit in (54) we use the series expansion approach again

$$\begin{array}{cc}\hfill {\tilde{u}}_{\alpha ,Ste}(y,\tau )& =1-{\displaystyle \frac{{\int }_0^{y/{(\tau /Ste)}^{1/(1+\alpha )}}{\sigma }_{\alpha }\left(w\right)dw}{{\int }_0^{{\xi }_{Ste}}{\sigma }_{\alpha }\left(w\right)dw}}\hfill \\ \hfill & =1-{\displaystyle \frac{{\displaystyle \sum _{n=0}^{\infty }}\frac{c_n{(-1)}^n}{{(1+\alpha )}^n}\frac{y^{(n+1)(1+\alpha )-1}}{(n+1)(1+\alpha )-1}{\left(\frac{Ste}{\tau }\right)}^{\frac{(n+1)(1+\alpha )-1}{1+\alpha }}}{{\displaystyle \sum _{n=0}^{\infty }}\frac{c_n{(-1)}^n}{{(1+\alpha )}^n}\frac{{\xi }_{Ste}^{(n+1)(1+\alpha )-1}}{(n+1)(1+\alpha )-1}}}\hfill \\ \hfill & =1-{\displaystyle \frac{y^{\alpha }}{{\tau }^{\frac{\alpha }{1+\alpha }}}}{\left({\displaystyle \frac{St{e}^{\frac{1}{1+\alpha }}}{{\xi }_{Ste}}}\right)}^{\alpha }{\displaystyle \frac{{\displaystyle \sum _{n=0}^{\infty }}\frac{c_n{(-1)}^n}{{(1+\alpha )}^n}\frac{y^{n(1+\alpha )}}{(n+1)(1+\alpha )-1}{\left(\frac{Ste}{\tau }\right)}^{\frac{n(1+\alpha )}{1+\alpha }}}{{\displaystyle \sum _{n=0}^{\infty }}\frac{c_n{(-1)}^n}{{(1+\alpha )}^n}\frac{{\xi }_{Ste}^{n(1+\alpha )}}{(n+1)(1+\alpha )-1}}}\hfill \end{array}$$

(67)

By applying (61), (65) and the uniform convergence of the series we conclude that

$$\underset{Ste\to 0}{lim}{\tilde{u}}_{\alpha ,Ste}(y,\tau )=1-{\displaystyle \frac{y^{\alpha }}{{\tau }^{\frac{\alpha }{1+\alpha }}}}{\left({\displaystyle \frac{1}{\mathrm{\Gamma }{(2+\alpha )}^{1/(1+\alpha )}}}\right)}^{\alpha }={\tilde{u}}_{\alpha ,0}(y,\tau ).$$

(68)

□

5. Computing the Parameter for the Time Fractional Case

5.1. Close Solutions

Let us start by introducing the special functions involved in the representation of the solution for the time fractional case. For more details we refer the reader to [25]

Definition 2. 

For every $y\in [0,1]$, the upper incomplete Beta function of parameters $\alpha ,\beta >0$ is given by

$$\overline{B}(y;\alpha ,\beta )={\int }_y^1{z}^{\alpha -1}{(1-z)}^{\beta -1}\mathrm{d}z.$$

(69)

Also, the regularized upper incomplete Beta functions is defined by

$$\overline{I}(y;\alpha ,\beta )={\displaystyle \frac{\overline{B}(y;\alpha ,\beta )}{B(\alpha ,\beta )}}.$$

(70)

Here, $B(\alpha ,\beta )$ is the classical Beta function defined by $B(\alpha ,\beta )={\int }_0^1{z}^{\alpha -1}{(1-z)}^{\beta -1}\mathrm{d}z$.

Remark 7. 

Note that the former functions are not symmetric respect to the parameters α and β, in contrasts to the classical Beta function which verifies that

$$B(\alpha ,\beta )={\displaystyle \frac{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }\left(\beta \right)}{\mathrm{\Gamma }(\alpha +\beta )}}.$$

(71)

Proposition 3. 

For all $0<\alpha <1$, it holds that

$${\int }_0^1\overline{I}\left(z^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dz={\displaystyle \frac{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }\left(1-\frac{\alpha }{2}\right)}{\mathrm{\Gamma }\left(\frac{\alpha }{2}\right)}}.$$

(72)

Proof. 

Applying Fubini’s theorem, we have that

$$\begin{array}{cc}\hfill {\int }_0^1\overline{I}\left(z^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dz& ={\displaystyle \frac{1}{B\left(\frac{\alpha }{2},1-\alpha \right)}}{\int }_0^1{\int }_{z^{2/\alpha }}^1{w}^{{\displaystyle \frac{\alpha }{2}}-1}{(1-w)}^{-\alpha }dwdz\hfill \\ \hfill & ={\displaystyle \frac{1}{B\left(\frac{\alpha }{2},1-\alpha \right)}}{\int }_0^1{\int }_0^{w^{\alpha /2}}{w}^{{\displaystyle \frac{\alpha }{2}}-1}{(1-w)}^{-\alpha }dzdw\hfill \end{array}$$

(73)

and the thesis holds by using the Beta function definition and (71). □

We start with the quasi-stationary problem. That is, we consider the limit problem obtained from (46) by making $c\to 0$ (or equivalently $Ste\to 0$), which is the Time Fractional Quasi-Stationary Problem: Find the pair of functions $\tilde{u}:Q_{\tilde{s}}\to \mathbb{R}$ and $\tilde{s}:{\mathbb{R}}_0^{+}\to \mathbb{R}$ with enough regularity such that

$$\begin{array}{ccc}\left(i\right)\hfill & {\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha ){(\tau -\tilde{h}\left(y\right))}^{\alpha }}}={\displaystyle \frac{{\partial }^2}{\partial y^2}}\tilde{u}(y,\tau ),\hfill & 0<y<\tilde{s}\left(\tau \right),\hspace{0.17em}0<\tau ,\hspace{0.17em}\hspace{0.17em}\hfill \\ \left(ii\right)\hfill & \tilde{s}\left(0\right)=0,\hfill & \\ \left(iii\right)\hfill & \tilde{u}(0,\tau )=1,\hfill & 0<\tau ,\hfill \\ \left(iv\right)\hfill & \tilde{u}(\tilde{s}\left(\tau \right),\tau )=0,\hfill & 0<\tau ,\hfill \\ \left(v\right)\hfill & {\displaystyle \underset{y\to \tilde{s}{\left(\tau \right)}^{-}}{lim}{\tilde{u}}_y(y,\tau )=0,}\hfill & 0<\tau ,\hfill \end{array}$$

(74)

where $\tilde{h}:{\mathbb{R}}_0\to \mathbb{R}$ is the function given by

$$\tilde{h}\left(y\right)={\tilde{s}}^{-1}\left(y\right),y\ge 0.$$

(75)

Problem (74) was already addressed in [21], where the following solution was obtained:

$$\begin{array}{c}{\displaystyle {\tilde{u}}_{\alpha ,0}(y,\tau )=\frac{2}{\mathrm{\Gamma }(\alpha +1)}\frac{\mathrm{\Gamma }(1+\alpha /2)}{\mathrm{\Gamma }(1-\alpha /2)}{\int }_{\frac{y}{\sqrt{\frac{2}{\mathrm{\Gamma }(\alpha +1)}}{\tau }^{\alpha /2}}}^1\overline{I}(z^{2/\alpha };1-\alpha ,\alpha /2)\mathrm{d}z,0\le y\le \tilde{s}\left(\tau \right),\hspace{0.17em}0<\tau <T,}\hfill \\ {\tilde{s}}_{\alpha ,0}\left(\tau \right)=\sqrt{{\displaystyle \frac{2}{\mathrm{\Gamma }(1+\alpha )}}}{\tau }^{\alpha /2}.\hfill \end{array}$$

(76)

Regarding the non-stationary problem (46), we must refer first to the work of Kubica and Ryszewska [14] where a similarity solution was obtained. However, unlike the solution here, this solution was not presented in a compact form in terms of calculable transcendental functions. To more clearly see this, we present a selfsimilar solution following the lines in [14]. To begin, let us recall the formulation of Problem (46).

$$\begin{array}{ccc}\left(i\right)\hfill & {Ste}_{\tilde{h}\left(y\right)}^C{D}_{\tau }^{\alpha }\tilde{u}(y,\tau )+{\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha ){\left(\tau -\tilde{h}\left(y\right)\right)}^{\alpha }}}={\tilde{u}}_{yy}(y,\tau ),\hfill & 0<y<\tilde{s}\left(\tau \right),\hspace{0.17em}0<\tau ,\hfill \\ \left(ii\right)\hfill & \tilde{u}(0,\tau )=1,\hfill & 0<\tau ,\hfill \\ \left(iii\right)\hfill & \tilde{u}(\tilde{s}\left(\tau \right),\tau )=0,\hfill & 0<\tau ,\hfill \\ \left(iv\right)\hfill & \tilde{s}\left(0\right)=0,\hfill & \\ \left(v\right)\hfill & {\tilde{u}}_y(\tilde{s}{\left(\tau \right)}^{-},\tau )=0,\hfill & 0<\tau .\hfill \end{array}$$

(77)

Let us define the one variable function in terms of a self-similar variable (which is different than the one considered in [14]).

$$F\left(\mu \right):=\tilde{u}(y,\tau ),\mu ={\displaystyle \frac{y}{a{\tau }^{\alpha /2}}}.$$

(78)

And let the free boundary be given by

$$\tilde{s}\left(\tau \right)=a{\tau }^{\alpha /2}(\mathrm{and}\hspace{0.17em}\hspace{0.17em}\tilde{h}\left(y\right)={(y/a)}^{2/\alpha }).$$

(79)

Remark 8. 

We have made a simplification in the notation by denoting $a:={\theta }_{\alpha ,Ste}$, in the aim to simplify the reading of the next computes. And the correct notation will be recovered at the end of this section.

We have

$$\begin{array}{cc}\hfill {\tilde{u}}_{\tau }(y,\tau )& =F^{\prime }\left(\mu \right)\left(-{\displaystyle \frac{\alpha }{2}}\right){\displaystyle \frac{y}{a{\tau }^{\alpha /2+1}}}\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill {\tilde{u}}_y(y,\tau )& =F^{\prime }\left(\mu \right){\displaystyle \frac{1}{a{\tau }^{\alpha /2}}}\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill {\tilde{u}}_{yy}(y,\tau )& =F^{\prime \prime }\left(\mu \right){\displaystyle \frac{1}{a^2{\tau }^{\alpha }}}\hfill \end{array}$$

(82)

From the boundary condition (77)$-\left(iv\right)$ and making the substitution $\eta =p\tau$, it holds that

$$\begin{array}{cc}\hfill {}_{\tilde{h}\left(y\right)}{}^{C}D_{\tau }^{\alpha }\tilde{u}(y,\tau )& ={}_{\tilde{h}\left(y\right)}I_{\tau }^{1-\alpha }{\tilde{u}}_{\tau }(y,\tau )\hfill \\ \hfill & ={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\tilde{h}\left(y\right)}^{\tau }{F}^{\prime }\left({\displaystyle \frac{y}{a{\eta }^{\alpha /2}}}\right)\left(-{\displaystyle \frac{\alpha }{2}}\right){\displaystyle \frac{y}{a{\eta }^{\alpha /2+1}}}{(\tau -\eta )}^{-\alpha }\mathrm{d}\eta \hfill \\ \hfill & =-{\displaystyle \frac{\alpha /2}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{{\left({\displaystyle \frac{y}{a{\tau }^{\alpha /2}}}\right)}^{2/\alpha }}^1{F}^{\prime }\left({\displaystyle \frac{y}{a{\tau }^{\alpha /2}}}{\displaystyle \frac{1}{p^{\alpha /2}}}\right){\displaystyle \frac{y}{a{\left(\tau p\right)}^{\alpha /2+1}}}{t}^{-\alpha }{(1-p)}^{-\alpha }\tau \mathrm{d}p\hfill \\ \hfill & =-{\tau }^{-\alpha }{\displaystyle \frac{\alpha /2}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{{\mu }^{2/\alpha }}^1{F}^{\prime }\left(\mu p^{-\alpha /2}\right)\mu p^{-\alpha /2-1}{(1-p)}^{-\alpha }\mathrm{d}p.\hfill \end{array}$$

Also

$$\begin{array}{c}\hfill {\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha ){(\tau -\tilde{h}\left(y\right))}^{\alpha }}}={\displaystyle \frac{{\tau }^{-\alpha }}{\mathrm{\Gamma }(1-\alpha )}}{\displaystyle \frac{1}{{(1-{\left(\frac{y}{a{\tau }^{\alpha /2}}\right)}^{2/\alpha })}^{\alpha }}}={\displaystyle \frac{{\tau }^{-\alpha }}{\mathrm{\Gamma }(1-\alpha )}}{\displaystyle \frac{1}{{(1-{\mu }^{2/\alpha })}^{\alpha }}}.\end{array}$$

Then u is a solution to the fractional PDE (77)$-\left(i\right)$ if and only if F is a solution to the fractional ODE

$$-a^{-2}{\tau }^{-\alpha }{F}^{\prime \prime }\left(\mu \right)={\displaystyle \frac{Ste·\alpha /2}{\mathrm{\Gamma }(1-\alpha )}}{\tau }^{-\alpha }\mu {\int }_{{\mu }^{2/\alpha }}^1{F}^{\prime }\left(\mu p^{-\alpha /2}\right)p^{-\alpha /2-1}{(1-p)}^{-\alpha }\mathrm{d}p-{\displaystyle \frac{{\tau }^{-\alpha }}{\mathrm{\Gamma }(1-\alpha ){(1-{\mu }^{2/\alpha })}^{\alpha }}}$$

(83)

or equivalently

$$-F^{\prime \prime }\left(\mu \right)={\displaystyle \frac{Ste·a^2}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\mu }^1{F}^{\prime }\left(w\right){\left(1-{\left({\displaystyle \frac{\mu }{w}}\right)}^{2/\alpha }\right)}^{-\alpha }\mathrm{d}w-{\displaystyle \frac{a^2}{\mathrm{\Gamma }(1-\alpha ){(1-{\mu }^{2/\alpha })}^{\alpha }}}.$$

(84)

Note that condition (77)$-\left(iv\right)$ implies that $F^{\prime }\left(1\right)=0$. Then, integrating (84) between $\mu$ and 1 yields

$$F^{\prime }\left(\mu \right)={\displaystyle \frac{Ste·a^2}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\mu }^1{\int }_z^1{F}^{\prime }\left(w\right){\left(1-{\left({\displaystyle \frac{z}{w}}\right)}^{2/\alpha }\right)}^{-\alpha }\mathrm{d}w\mathrm{d}z-{\displaystyle \frac{a^2}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\mu }^1{(1-z^{2/\alpha })}^{-\alpha }\mathrm{d}z.$$

(85)

Let us work with the integrals in the r.h.s of (85) starting from the second one.

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\mu }^1{(1-z^{2/\alpha })}^{-\alpha }\mathrm{d}z& ={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{{\mu }^{2/\alpha }}^1{\displaystyle \frac{\alpha }{2}}{y}^{\alpha /2-1}{(1-y)}^{-\alpha }\mathrm{d}y\hfill \\ \hfill & ={\displaystyle \frac{\mathrm{\Gamma }(1+\alpha /2)}{\mathrm{\Gamma }(1-\alpha /2)}}\overline{I}({\mu }^{2/\alpha };\alpha /2,1-\alpha ).\hfill \end{array}$$

(86)

In the first one we apply Fubini and proceed as in (86) to get

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\mu }^1{\int }_z^1{F}^{\prime }\left(w\right){\left(1-{\left({\displaystyle \frac{z}{w}}\right)}^{2/\alpha }\right)}^{-\alpha }\mathrm{d}w\mathrm{d}z& ={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\mu }^1{F}^{\prime }\left(w\right){\int }_{\mu }^w{\left(1-{\left({\displaystyle \frac{z}{w}}\right)}^{2/\alpha }\right)}^{-\alpha }\mathrm{d}z\mathrm{d}w\hfill \\ \hfill & ={\displaystyle \frac{1}{\mathrm{\Gamma }(1-\alpha )}}{\int }_{\mu }^1{F}^{\prime }\left(w\right){\int }_{(\mu /w)}^1{\left(1-y^{2/\alpha }\right)}^{-\alpha }w\mathrm{d}y\mathrm{d}w\hfill \\ \hfill & ={\displaystyle \frac{\mathrm{\Gamma }(1+\alpha /2)}{\mathrm{\Gamma }(1-\alpha /2)}}{\int }_{\mu }^1{F}^{\prime }\left(w\right)w\overline{I}\left({\left({\displaystyle \frac{\mu }{w}}\right)}^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)\mathrm{d}w.\hfill \end{array}$$

(87)

By replacing (86) and (87) in (85) we obtain

$$F^{\prime }\left(\mu \right)={\displaystyle \frac{a^2\mathrm{\Gamma }(1+\alpha /2)}{\mathrm{\Gamma }(1-\alpha /2)}}\left[Ste {\int }_{\mu }^1{F}^{\prime }\left(w\right)w\overline{I}\left({\left({\displaystyle \frac{\mu }{w}}\right)}^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)\mathrm{d}w-\overline{I}({\mu }^{2/\alpha };\alpha /2,1-\alpha )\right].$$

(88)

Now, proceeding as in [14] we define the function

$$\begin{array}{c}\hfill G\left(\mu \right)=a^2{\displaystyle \frac{\mathrm{\Gamma }(1+\alpha /2)}{\mathrm{\Gamma }(1-\alpha /2)}}\overline{I}({\mu }^{2/\alpha };\alpha /2,1-\alpha )\end{array}$$

(89)

and the operator $L:C^1[0,1]\to C^1[0,1]$ such that

$$\begin{array}{c}\hfill \left(Lj\right)\left(\mu \right)=Ste a^2{\displaystyle \frac{\mathrm{\Gamma }(1+\alpha /2)}{\mathrm{\Gamma }(1-\alpha /2)}} {\int }_{\mu }^{1}j\left(w\right)w\overline{I}\left({\left({\displaystyle \frac{\mu }{w}}\right)}^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)\mathrm{d}w,\forall \hspace{0.17em}j\in C^1[0,1].\end{array}$$

(90)

Then $F^{\prime }$ verifies

$$F^{\prime }\left(\mu \right)=\left(L{F}^{\prime }\right)\left(\mu \right)-G\left(\mu \right),\forall \mu \in [0,1].$$

(91)

It is important to recall that we have written the operator L and function G in terms of the incomplete Beta function, but the work and estimations done in Section 5.3 and 5.4 of [14] are valid for these expressions. The recursive argument is also valid, consisting of applying the operator L to both sides of equation (91) and use eq. (91) to recover a new equation for $F^{\prime }$$n-$times, yielding that

$${\displaystyle F^{\prime }\left(\mu \right)=\left(L^n{F}^{\prime }\right)\left(\mu \right)-\sum _{k=0}^n\left(L^{k}G\right)\left(\mu \right),\forall \hspace{0.17em}n\in \mathbb{N}.}$$

(92)

In [14,Secc.5.2] it was proved that taking the limit when $n\to \infty$ in (92) yields that

$${\displaystyle F^{\prime }\left(\mu \right)=-\sum _{n=0}^{\infty }\left(L^{n}G\right)\left(\mu \right).}$$

(93)

And integrating between $\mu$ and 1 and using condition (46)$-\left(iii\right)$ we get an expression for F which is

$${\displaystyle F\left(\mu \right)={\int }_{\mu }^1\sum _{n=0}^{\infty }\left(L^{n}G\right)\left(z\right)\mathrm{d}z.}$$

(94)

The uniform convergence of the series in (93) as well as the regularity of function F were analyzed in the mentioned paper. We look now for the parameter a, by using the boundary condition (77)$-\left(ii\right)$. In other words, we look for a positive value a such that

$${\displaystyle 1={\int }_0^1\sum _{n=0}^{\infty }\left(L^{n}G\right)\left(z\right)\mathrm{d}z,}$$

(95)

where from (90) and denoting $C_{\alpha }:={\displaystyle \frac{\mathrm{\Gamma }(1+\alpha /2)}{\mathrm{\Gamma }(1-\alpha /2)}}$, the operators $L^n$ are given by

$$L^{n}G\left(\mu \right)=Ste\hspace{0.17em}{a}^2{C}_{\alpha }\hspace{0.17em}{\int }_{\mu }^1\left[L^{n-1}G\right]\left(w\right)w\overline{I}\left({\left({\displaystyle \frac{\mu }{w}}\right)}^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)\mathrm{d}w,$$

for $n\in \mathbb{N}$ and we define $L^0$ as $L^{0}G\left(\mu \right)=G\left(\mu \right)$.

In order to obtain an expression in terms of the powers of the parameter $Ste·a^2{C}_{\alpha }$, we redefine the operators $L^n$ as follows. First we define

$${\tilde{L}}^{0}j\left(\mu \right)=\overline{I}({\mu }^{2/\alpha };\alpha /2,1-\alpha ),{\tilde{L}}^{n}j\left(\mu \right)={\int }_{\mu }^1\left[{\tilde{L}}^{(n-1)}j\right]\left(w\right)w\overline{I}\left({\left({\displaystyle \frac{\mu }{w}}\right)}^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)\mathrm{d}w$$

(96)

and then we replace it in (95), multiply both sides by $Ste$, and we deduce that the prefactor a verifies that

$${\displaystyle Ste=Ste·a^2{C}_{\alpha }\sum _{n=0}^{\infty }{b}_n{(Ste\hspace{0.17em}{a}^2{C}_{\alpha })}^n,b_n={\int }_0^1{\tilde{L}}^{n}G\left(z\right)\mathrm{d}z.}$$

(97)

Now, define the function

$${\displaystyle \mathcal{H}\left(y\right)=\sum _{n=0}^{\infty }{b}_n{y}^n,0\le y\le 1.}$$

(98)

We know from [14][Secc- 5.2] that the series in (98) is absolutely convergent and thus the function $\mathcal{H}$ is a well defined continuous function. Moreover, $\mathcal{H}$ is an increasing and positive function. Thus, the equation

$$y\mathcal{H}\left(y\right)=Ste$$

(99)

admits a unique positive solution, named $y_{Ste}$ and we can recover a asking $y_{Ste}$ that

$$y_{Ste}=Ste·a^2{C}_{\alpha }.$$

More precisely, recalling the definition of $C_{\alpha }$, we state that the prefactor a is given by

$$a=\sqrt{{\displaystyle \frac{\mathrm{\Gamma }(1-\alpha /2)}{\mathrm{\Gamma }(1+\alpha /2)}}{\displaystyle \frac{y_{Ste}}{Ste}}}.$$

(100)

From the previous analysis we can state that the solution to problem (77) is given by the pair

$$\begin{array}{cc}\hfill {\tilde{u}}_{\alpha ,Ste}(y,\tau )& {\displaystyle ={\int }_{y/a{\tau }^{\alpha /2}}^1\sum _{n=0}^{\infty }\left(L^{n}G\right)\left(z\right)\mathrm{d}z,}\hfill \end{array}$$

(101)

$$\begin{array}{cc}\hfill {\tilde{s}}_{\alpha ,Ste}\left(\tau \right)& ={\theta }_{\alpha ,Ste}{\tau }^{\alpha /2}.\hfill \end{array}$$

(102)

where ${\theta }_{\alpha ,Ste}:=a$ is defined by (100) and $y_{Ste}$ is the unique solution to equation (99).

5.2. Computing ${\theta }_{\alpha ,Ste}$

In this section, we outline the elements needed to compute the prefactor a. i.e. the coefficients $b_n$. After this computation, the value of a is achieved using the nonlinear equation (99) in a similar way than before (see [5] for details on the computation).

For that purpose, let us focus on the coefficients $b_n$ defined in (97).

Proposition 4. 

For every $n\in \mathbb{N}$ the coefficients $b_n$ defined in (97) verify that

$$b_n=\prod _{i=1}^{n+1}\left({\int }_0^1{t}^{2(i-1)}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)\mathrm{d}t\right).$$

(103)

Proof. 

We prove first that, for every $k=0,...,n-1$

$$\begin{array}{cc}\hfill & {\int }_0^1{\tilde{L}}^{n-k}G\left(z\right)z^{2k}dz={\int }_0^1{\tilde{L}}^{n-(k+1)}G\left(z\right)z^{2(k+1)}dz·{\int }_0^1{t}^{2k}\hspace{0.17em}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dt.\hfill \end{array}$$

(104)

In fact, applying (96), Fubini’s Theorem and making the substitution $z=tw$ in the penultimate equality, we get

$$\begin{array}{cc}\hfill {\int }_0^1{\tilde{L}}^{n-k}G\left(z\right)z^{2k}dz& ={\int }_0^1\tilde{L}\left({\tilde{L}}^{n-k-1}G\right)\left(z\right)z^{2k}dz\hfill \\ \hfill & ={\int }_0^1{\int }_z^1{\tilde{L}}^{n-k-1}G\left(w\right)w\overline{I}\left({\left({\displaystyle \frac{z}{w}}\right)}^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)z^{2k}dwdz\hfill \\ \hfill & ={\int }_0^1{\int }_0^w{\tilde{L}}^{n-k-1}G\left(w\right)w\overline{I}\left({\left({\displaystyle \frac{z}{w}}\right)}^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)z^{2k}dz\hspace{0.17em}dw\hfill \\ \hfill & ={\int }_0^1{\tilde{L}}^{n-k-1}G\left(w\right)w{\int }_0^w\overline{I}\left({\left({\displaystyle \frac{z}{w}}\right)}^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)z^{2k}dzdw\hfill \\ \hfill & ={\int }_0^1{\tilde{L}}^{n-k-1}G\left(w\right)w{\int }_0^1\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)t^{2k}{w}^{2k}wdtdw\hfill \\ \hfill & ={\int }_0^1{\tilde{L}}^{n-(k+1)}G\left(w\right)w^{2(k+1)}dw{\int }_0^1{t}^{2k}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dt.\hfill \end{array}$$

(105)

Repeating the former argument $n-$times and using (96) we obtain

$$\begin{array}{cc}\hfill b_n& ={\int }_0^1{\tilde{L}}^{n}G\left(z\right)dz\hfill \\ \hfill & ={\int }_0^1{\tilde{L}}^{n-1}G\left(w\right)w^{2}dw{\int }_0^1{t}^{2.0}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dt\hfill \\ \hfill & ={\int }_0^1{\tilde{L}}^{n-2}G\left(w\right)w^{2.2}dw{\int }_0^1{t}^{2.1}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dt·{\int }_0^1{t}^{2.0}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dt\hfill \\ \hfill & =...\hfill \\ \hfill & ={\int }_0^1{\tilde{L}}^{0}G\left(w\right)w^{2n}dw\prod _{i=1}^n\left({\int }_0^1{t}^{2(i-1)}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dt\right)\hfill \\ \hfill & ={\int }_0^1\overline{I}\left(w^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)w^{2n}dw·\prod _{i=1}^n\left({\int }_0^1{t}^{2(i-1)}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dt\right)\hfill \\ \hfill & =\prod _{i=1}^{n+1}\left({\int }_0^1{t}^{2(i-1)}\overline{I}\left(t^{2/\alpha },{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dt\right),\hfill \end{array}$$

(106)

which completes the proof.

□

We also present the following property needed for next section.

Proposition 5. 

For all $n\in {\mathbb{N}}_0$, it holds that $b_n>0$.

Proof. 

It immediately follows from the positiveness of

$${\tilde{L}}^{0}G\left(z\right)=\overline{I}\left(z^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)>0,\forall z\ne 1,$$

and formula (103). □

5.3. Analysis of the Convergence to the Quasi-Stationary Case

Following a similar approach presented in Section 4.3, and recalling (47), the analysis of convergence will be state by making $Ste\to 0$.

Theorem 2. 

The prefactor of the free boundary in equation () associated with the Time Fractional Stefan Problem (77) converges to the prefactor of the quasi-stationary Time Fractional Stefan Problem (74) as the Stefan number tends to zero. That is,

$$\underset{Ste\to 0}{lim}{\theta }_{\alpha ,Ste}=\sqrt{{\displaystyle \frac{2}{\mathrm{\Gamma }(1+\alpha )}}}.$$

(107)

Proof. 

Let $y_{Ste}$ the unique positive solution to $y\mathcal{H}\left(y\right)=Ste$.

By Proposition 5, we deduce that $y_{Ste}\mathcal{H}\left(y_{Ste}\right)>0$, and $\mathcal{H}\left(y_{Ste}\right)\ge b_0>0$. Then, from this comment and equality

$$\underset{Ste\to 0}{lim}{y}_{Ste}\mathcal{H}\left(y_{Ste}\right)=\underset{Ste\to 0}{lim}Ste=0,$$

it holds that $\underset{Ste\to 0}{lim}{y}_{Ste}=0$.

We can compute now the following limit

$$\begin{array}{cc}\hfill \underset{Ste\to 0}{lim}{\displaystyle \frac{y_{Ste}}{Ste}}& =\underset{Ste\to 0}{lim}{\displaystyle \frac{y_{Ste}}{y_{Ste}\mathcal{H}\left(y_{Ste}\right)}}=\underset{Ste\to 0}{lim}{\displaystyle \frac{1}{{\displaystyle \sum _{n=0}^{\infty }{b}_n{y}_{Ste}^n}}}={\displaystyle \frac{1}{b_0}}\hfill \\ \hfill & ={\displaystyle \frac{1}{{\int }_0^1\overline{I}\left(z^{2/\alpha };\frac{\alpha }{2},1-\alpha \right)dz}}={\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\alpha }{2}\right)}{\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }\left(1-\frac{\alpha }{2}\right)}}.\hfill \end{array}$$

(108)

Returning to (100), and using (108) and the Gamma function property $\mathrm{\Gamma }(z+1)=z\mathrm{\Gamma }\left(z\right)$ several times we get

$$\underset{Ste\to 0}{lim}{\theta }_{\alpha ,Ste}=\underset{Ste\to 0}{lim}\sqrt{{\displaystyle \frac{\mathrm{\Gamma }\left(1-\frac{\alpha }{2}\right)}{\mathrm{\Gamma }\left(1+\frac{\alpha }{2}\right)}}{\displaystyle \frac{y_{Ste}}{Ste}}}=\sqrt{{\displaystyle \frac{2}{\mathrm{\Gamma }(1+\alpha )}}}.$$

(109)

□

Corollary 2. 

The solution (101)-() to the Time Fractional Stefan Problem (77) converges to the solution (76) to the quasi-stationary Time Fractional Stefan problem (74) when the Stefan number approaches zero.

Proof. 

First, note that, being $\overline{I}\left(z;{\displaystyle \frac{\alpha }{2}},1-\alpha \right)\le 1$ for all $0\le z\le 1$, it follows that

$$0\le L^{0}G\left(z\right)={\theta }_{\alpha ,Ste}^2{C}_{\alpha }\overline{I}\left(z^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)\le {\theta }_{\alpha ,Ste}^2{C}_{\alpha },$$

By induction, we deduce that

$$0\le L^{n}G\left(z\right)\le St{e}^n{\theta }_{\alpha ,Ste}^{2(n+1)}{C}_{\alpha }^{n+1}.$$

Then,

$$0\le {\tilde{u}}_{Ste}(y,\tau )={\int }_{{\displaystyle \frac{y}{{\theta }_{\alpha ,Ste}{\tau }^{\alpha /2}}}}^1\sum _{n=0}^{\infty }\left(L^{n}G\right)\left(z\right)dz\le {\int }_{{\displaystyle \frac{y}{{\theta }_{\alpha ,Ste}{\tau }^{\alpha /2}}}}^1\sum _{n=0}^{\infty }St{e}^n{\theta }_{\alpha ,Ste}^{2(n+1)}{C}_{\alpha }^{n+1}dz,$$

and for $Ste<1$ the series in the right-hand-side converges absolutely. Hence, using that $L^{n}G\left(z\right)\to 0$ for all $n\in \mathbb{N}$ when $Ste\to 0$, we obtain that

$$\begin{array}{cc}\hfill \underset{Ste\to 0}{lim}{\tilde{u}}_{Ste}(y,\tau )& =\underset{Ste\to 0}{lim}{\int }_{{\displaystyle \frac{y}{a{\tau }^{\alpha /2}}}}^1\sum _{n=0}^{\infty }\left(L^{n}G\right)\left(z\right)dz\hfill \\ \hfill & ={\displaystyle \frac{\mathrm{\Gamma }\left(1+\frac{\alpha }{2}\right)}{\mathrm{\Gamma }\left(1-\frac{\alpha }{2}\right)}}{\displaystyle \frac{2}{\mathrm{\Gamma }\left(1+\alpha \right)}}{\int }_{{\displaystyle \frac{y}{\sqrt{\frac{2}{\mathrm{\Gamma }\left(1+\alpha \right)}}{\tau }^{\alpha /2}}}}^1\overline{I}\left(z^{2/\alpha };{\displaystyle \frac{\alpha }{2}},1-\alpha \right)dz,\hfill \end{array}$$

(110)

and the convergence ${\tilde{u}}_{\alpha ,Ste}\to {\tilde{u}}_{\alpha ,0}$ is proven.

Going back to definition (100), applying (108) and properties of Gamma function we obtain that

$$\underset{Ste\to 0}{lim}{\theta }_{\alpha ,Ste}{\tau }^{\alpha /2}=\underset{Ste\to 0}{lim}\sqrt{{\displaystyle \frac{\mathrm{\Gamma }\left(1-\frac{\alpha }{2}\right)}{\mathrm{\Gamma }\left(1+\frac{\alpha }{2}\right)}}{\displaystyle \frac{y_{Ste}}{Ste}}}{\tau }^{\alpha /2}=\sqrt{{\displaystyle \frac{2}{\mathrm{\Gamma }(1+\alpha )}}}{\tau }^{\alpha /2}.$$

(111)

□

6. Conclusions

We have presented selfsimilar solutions for one-dimensional time and space fractional Stefan melting problems; problems directly derived from thermodynamically consistent balance laws. It is well known that, when a constant temperature boundary ($x=0$) is applied and the initial position of the sharp interface is the origin ($s\left(0\right)=0$), these solution exhibit sub- and super- diffusion behaviors respectively, i.e., the melt front advance is determined as the product of a prefactor and time to a power n, different from the diffusion value of ${\displaystyle \frac{1}{2}}$. In both the fractional time and space cases, the time exponent is given in terms of known functions of the orders of the space and time derivatives in the problem formulations. The main contribution of the current work has been to, for the first time, make explicit computations for the values of the prefactors. We expect that our analysis will be useful when proposing fractional models for anomalous diffusion. Besides, the exact solutions (prefactor values) computed here will be a strong tool for testing numerical methods for fractional free boundary problems.