1 Introduction

The phase-change problems (Stefan problems) involve one or more moving boundaries that separate the different phases of the material. These problems arise in many natural and manufacturing phenomena. The applicability of these problems and the presence of the moving boundaries make it interesting from industrial as well as mathematical point of views. Moreover, the presence of moving boundary is also a key reason for these problems to be a non-linear even in its simplest form. In the classical Stefan problems (Gupta, 2017), the velocity of phase change material has been assumed as zero, and the thermal coefficients have been taken as constants. But it is not always appropriate with the many materials. Hence, the variable thermal coefficients have been attracted many scientists and engineers in the field of phase-change processes (Oliver and Sunderland, 1987; Rogers and Broadbridge, 1988; Ramos et al., 1994; Tritscher and Broadbridge, 1994; Broadbridge and Pincombe, 1996). Mondal et al. (2015) also assumed temperature-dependent thermal conductivity in the study of thermal radiation on an unsteady MHD axisymmetric stagnation point flow over a shrinking sheet. The temperature-dependent thermal coefficients in the one-dimensional phase change problem are considered by Briozzo et al. (2007), and they discussed the exact solution to the problem. Many other authors (Briozzo and Natale, 2015; Briozzo and Natale, 2017) also took variable thermal coefficients in their study and presented either exact or approximate solutions or both. Khader (2016) also considered temperature-dependent thermal conductivity in the problem of flow of Newtonian fluid over an impermeable stretching sheet. Ceretani et al. (2018) assumed thermal conductivity which linearly varies with temperature and Robin boundary condition in a phase change problem, and presented the similarity solution of the problem. Recently, Kumar et al. (2019) presented a Stefan problem involving thermal conductivity as a function of time and temperature. They discussed similarity solution for a limit case and approximate solution for the general case. Another Stefan problem containing temperature-dependent specific heat and thermal conductivity is mentioned by Kumar et al. (2018).

The occurrence of the phase change when the material is moving itself during the process is not much studied in the literature (Fila and Souplet, 2001; Lombardi and Tarzia, 2001). However, this type of physical situation may arise in many phase change processes. Recently, Turkyilmazoglu (2018) discussed the problems concerning melting and solidification processes that include moving phase change material (PCM). He has presented some analytical solutions to the problem by taking constant thermal coefficients. Singh et al. (2018) also discussed a freezing problem including convective boundary condition, moving phase change material and variable thermal coefficients.

Inspired by all these works, we consider a phase change problem related to melting process in which the phase change material moves with a speed $u$ in the positive direction of $x$ -axis which depends on time. Simultaneously, the variable thermal conductivity $k(T)$ and specific heat $c(T)$ are assumed in the problem. The constant melting temperature $T_m$ is assumed as initial temperature of the material. The mathematical model governing the process is given below:

Due to the complexity associated with the phase change problems, the establishment of analytical solutions always draws the attention of investigators. Some existing exact solutions of phase change problems can also be seen in the (Voller et al., 2004; Voller and Falcini, 2013; Zhou and Li-jiang, 2015). In this article, the similarity solution to the problem (mentioned in Eqs. (1)–(5)) is discussed for $\alpha =\beta$ and the uniqueness of this solution is also deliberated. Beside analytical method, we also present a spectral approach with the aid of shifted Legendre polynomials and collocation technique to the problem for all $\alpha$ and $\beta$ .

Due to the exponential rate of convergence, spectral methods have been used by many researchers to solve differential equations of various orders, and few of them are Canuto et al., 1988; Gottlieb and Hesthaven, 2001; Doha and Abd-Elhameed, 2006; Guo and Yan, 2009; Doha et al., 2012; Atabakzadeh et al., 2013; Hosseini et al., 2013; Agbaje et al., 2018. The applications of spectral relaxation method in fluid flow can be seen in Haroun et al., 2015a; Haroun et al., 2015b; Haroun et al., 2015c; Oyelakin et al., 2016; Haroun et al., 2016. Some other applications of spectral method, viz. spectral quasi linearization method, multi-domain quasilinearization method and multi-domain collocation method are reported by Mondal et al., 2016; Ahamed et al., 2016; Mahapatra et al., 2012; Almakki et al., 2018; Goqo et al., 2018; Noreldin et al., 2018; Mondal et al., 2019. Ahmadian et al. (2013) discussed the operational matrix based on shifted Legendre polynomials to solve the fuzzy differential equations of fractional order. Khader and Babatin (2014) used Legendre spectral collocation method to solve SIRC model and influenza A. Bhrawy and Zaky (2015) proposed shifted Jacobi collocation technique based on Jacobi operational matrix for Caputo fractional derivatives and solved the fractional order cable equation in one and two dimensional spaces. Abd-Elhameed et al. (2015) presented a new operational matrix method to solve the various boundary value problems by using the collocation method and Petrov-Galerkin method. Bhrawy and Zaky (2017a) developed an exponential order accurate Jacobi-Gauss-Lobatto collocation method to find the solution of the fractional Schrodinger equations in one and two dimensions. Bhrawy and Zaky (2017b) have derived new operational matrices of the shifted Jacobi polynomials for the fractional derivatives of Caputo and Riemann-Liouville types. They also used this development to find the solution of the variable-order Schrodinger equations. The spectral methods used to solve the differential/integral equations are characterized by the representation of the function, to be known, by a truncated series of smooth functions like polynomials. In this expansion, the main concern is to determine the unknown expansion coefficients. Doha et al. (2018) presented an article to give an overview of numerical difficulties while determining these coefficients and proposed the rich variety of tools to resolve these difficulties. Recently, Zaky (2018) produced an efficient method based on the Legendre-tau approximation for fractional Rayleigh-Stokes problems for a generalized second-grade fluid. Zaky et al. (2018) established a Legendre spectral-collocation technique for numerical solution of the distributed order fractional initial value problems and also discussed the convergence analysis of the method.

2 Shifted Legendre polynomials and its properties

In the interval $[-1,1]$ , the well-known classical Legendre polynomials $\{L_i(x);i=0,1,...\}$ are defined and orthogonal. In order to utilize these polynomials on the interval $[a,b]$ , we define the shifted Legendre polynomials $L_i^{\ast }(x)$ as

In this study, the following results of the shifted Legendre polynomials (Abd-Elhameed et al., 2015) are used:

1. Let us first define a space

Now, the function $u(x)\in L_0^2[a,b]$ can be written as:

In numerical calculation, the series given in Eq. (12) for the function $u(x)$ can be approximated as

• The derivative of $\mathbf{\Phi }(x)$ in matrix form is given as:

In Eq. (17), $H_i$ and $H_j$ are harmonic numbers which are defined as

• The relation between second order derivative of $\mathbf{\Phi }(x)$ and the operational matrix D is given by

3 Formulation of the problem

Substituting the following transformation

into the Eqs. (1)–(5), we obtain

Now, let us consider the following similarity variables

and assume that

Substituting Eqs. (29) and (30) into Eqs. (24)–(27), we get the following system involving ordinary differential equations:

Now, we substitute the following transformation

4 Approximate solution

To solve Eq. (36), we take an approximation of $y(\eta )$ in terms of the shifted Legendre polynomials as:

As mentioned in Section 2, $y^{\prime }(\eta )$ and $y^{″}(\eta )$ can be approximated as

and

Substituting the considered approximations of $y(\eta )$ , $y^{\prime }(\eta )$ and $y^{″}(\eta )$ into the Eq. (36), we get the residual, denoted by $R_N\left(\eta \right),$ corresponding to the Eq. (36) which is given below:

According to the spectral collocation method (Abd-Elhameed et al., 2015), we impose $R_N(\eta )=0$ at the first $(N+1)$ roots of $L_{N+1}^{\ast }(\eta )$ which produces $(N+1)$ non-linear algebraic equations involving $\left(N+2\right)$ unknowns ( $c_0,c_1,\cdots ,c_N$ and λ ). Beside these $(N+1)$ non-linear algebraic equations, one additional algebraic equation can be obtained with the aid of Eq. (38) which is

The obtained system of $(N+2)$ algebraic equations can be solved by an appropriate numerical technique like Newton-Raphson method to get all the $(N+2)$ unknowns. From Eq. (39), the approximate solution of $y(\eta )$ can be found, hence the $f(x,t)$ can be obtained by back substitution. After getting $\lambda$ , the moving phase front $s(t)$ can also be achieved with the help of Eq. (30).

5 Exact solution

First, we take $\alpha =\beta$ in the problem (31)–(34), hence the ordinary differential Eq. (31) becomes:

According to Singh et al., (2018), the solution of the Eq. (45) with the conditions (32) and (34) is given by

The Eqs. (29) and (46) give rise to the following equation:

Now, the Eqs. (33) and (46) yield the following equation:

We can calculate the unknown $\lambda$ from transcendental Eq. (48) if it exists. After getting $\lambda$ , we can easily find analytical expressions of $s(t)$ and $f(x,t)$ with the aid of (30) and (47), respectively. The existence and uniqueness of $\lambda$ satisfying the transcendental Eq. (48) is deliberated in the next section.

6 Existence and uniqueness

In order to show the existence and uniqueness of the exact solution discussed in Section 5, we consider the following function:

To prove the uniqueness of solution to the considered problem, it is enough to show that there exists a unique value of $\lambda$ in $\left(0,\infty \right)$ which satisfies the Eq. (49). From the Eq. (49), it is obvious that the function $h(\lambda )$ is continuous and differentiable on the interval $(0,\infty )$ . Moreover, $\underset{\lambda \to 0^{+}}{lim}h(\lambda )=-(2+\beta )$ and $\underset{\lambda \to \infty }{lim}h(\lambda )=\infty$ for the positive parameters $\beta$ and $\mathit{Ste}$ .

Now, the derivative of $h(\lambda )$ is given as

It is also observed that $h^{\prime }\left(\lambda \right)>0$ on the interval $(0,\infty )$ for $\mathit{Ste}>0,Pe\le \sqrt{2}$ . Hence, $h(\lambda )$ is a strictly increasing function on the interval $(0,\infty )$ . This shows that the equation $h(\lambda )=0$ has exactly one positive root in the interval $(0,\infty )$ for the positive values of parameters.

7 Results and discussions

We first discuss about the correctness of the solution obtained by spectral collocation technique (deliberated in Section 4) and the results thus found are depicted through the Tables 1 and 2. In this study, Wolfram Research (8.0.0) software and the following matrices are used in the calculation:

The approximate dimensionless temperature $f_A(x,t)$ , exact dimensionless temperature $f_E(x,t)$ and absolute error between them are revealed in Table 1 at $v=1.5$ , $t=1$ and $N=2$ for different $\alpha ,\beta$ , $\mathit{Pe}$ and Stefan number. Table 2 portrays the assessment for the accuracy of the exact solution $s_E(t)$ and the approximate solution $s_A(t)$ of the moving phase front at a constant thermal diffusivity $v=1.5$ by considering the matrices given in (51). Both the Tables endorse that the proposed approximate solutions of $f(x,t)$ and $s(t)$ are sufficiently near to the analytical solution discussed in Section 5 for $\alpha =\beta$ . Therefore, spectral collocation approach is a useful procedure to solve the moving boundary problems associated with phase change phenomenon.

With the aid of the procedure discussed in Section 3 and considering the operational matrix of order three; the Figs. 1–6 are plotted. The dependency of dimensionless temperature $f(x,t)$ on $x$ for three different Peclet numbers $\left(\mathit{Pe}=1,2,3\right)$ at $v=1.5$ and $t=5$ is displayed in Fig. 1. The Fig. 2 demonstrates the variations of $f(x,t)$ versus $x$ for three different $\alpha$ $\left(\alpha =1,2,3\right)$ at the fixed thermal diffusivity ( $v=1.5$ ) and time ( $t=5$ ). These figures represent that the dimensionless temperature are zero at $x=0$ and continuously increases till the last point of domain, i.e. $x=1$ . It is also detected from Figs. 1 and 2 that the rate of change of dimensionless temperature with respect to $x$ decreases when we increase either the Peclet numbers or the parameter $\alpha$ . But, the temperature distribution is more affected with the variation of Peclet numbers than the parameter $\alpha$ . Fig. 3 shows the temperature distribution within the domain for three cases, i.e., $\beta =1,2$ and 3 at $v=1.5$ and $t=5$ . This figure presents that the rate of change of temperature with respect to $x$ increases with the enhancement in the value of $\beta$ till some points of $x$ and after that the reverse situation is observed that can be seen in the figure.

Close

Fig. 1

Plot of f(x,t) for different values of Pe at Ste = 0.5 and α = β = 1.

Export to PPT

Close

Fig. 2

Plot of f(x,t) for different values of at α at Ste = 0.5, β = 1 and Pe = 1.

Export to PPT

Close

Fig. 3

Plot of f(x,t) for different values of β at Ste = 0.5, α = 1 and Pe = 1.

Export to PPT

Close

Fig. 4

Plot of s(t) for different values of Pe at Ste = 0.5 and α = β = 1.

Export to PPT

Close

Fig. 5

Plot of s(t) for different values of α at Ste = 0.5, β = 1 and Pe = 1.

Export to PPT

Close

Fig. 6

Plot of s(t) for different values of β at Ste = 0.5, α = 1 and Pe = 1.

Export to PPT

Figs. 4–6 depict the dependency of trajectory of phase front $s(t)$ on $\mathit{Pe},\alpha$ and $\beta$ for the fixed thermal diffusivity ( $v=1.5$ ) and Stefan number $(Ste=0.5)$ . From Fig. 4, it is seen that the phase front $s(t)$ propagates faster in the direction of phase change material as we increase the Peclet numbers $\left(\mathit{Pe}=1,2,3\right)$ . Moreover, the similar observations are established from Figs. 5 and 6, i.e. the enhancement in the movement of phase front $s\left(t\right)$ is found if we increase the value of either $\alpha$ or $\beta$ or both. If phase front moves more quickly with the increment in the value of a parameter then this indicates that the material melts/solidifies faster when we increase or decrease the same parameter. Therefore, the advancement of melting process is detected by the improvement of either $\alpha$ or $\beta$ or Pe form the Figs. 4–6. It is also observed that the effect of Peclet numbers in the progression of tracking of the phase front $s(t)$ is more than the parameters $\alpha$ or $\beta$ .

To show the accuracy of the proposed numerical solution with increasing the number of terms, Figs. 7 and 8 are plotted according to Zaky et al. (2018). Fig. 7 demonstrates the plot of log₁₀|error| of the moving interface factor ( $\lambda$ ) for different approximating polynomials of degree (N + 2) at the value of $\alpha =\beta =0$ , $\mathit{Ste}=0.5\text{and}Pe=1.5$ . In Fig. 8, we plot the graph of log₁₀|error| of the obtained numerical solution of $\theta (\eta )$ for different approximating polynomials of degree (N + 2) at the value of $\alpha =\beta =0$ , $\eta =0.5$ , $\mathit{Ste}=0.5\text{and}Pe=1.5$ . From Figs. 7 and 8, it is clear that the proposed solution converges rapidly as the degree of approximating polynomials or N increases.

Close

Fig. 7

Convergence of approximate λ.

Export to PPT

Close

Fig. 8

Convergence of approximate θ(η).

Export to PPT

8 Conclusion

This article included a problem of melting process in which it is assumed that the thermal coefficients depend on temperature and the phase change material moves with a variable velocity. The spectral collocation method is successfully applied to get an approximate solution to the problem. It is found that the spectral collocation method is a simple and sufficiently accurate scheme to develop the solution of the phase change problems. Hence, spectral collocation approach is an effective tool to get the solution of the problems associated to phase change processes. Beside this solution, an exact solution to the problem is established for a particular case, and it is revealed that there occurs a unique solution to the problem when $\mathit{Pe}\le \sqrt{2}$ . Like classical phase change problems Gupta (2017), this problem also consists of phase front $s(t)$ proportional to $\sqrt{t}$ . This article also described that the melting process is dependent on Peclet number $\mathit{Pe}$ , $\alpha$ and $\beta$ ; and the melting process becomes rapid as the parameter $\mathit{Pe}$ or $\alpha$ or $\beta$ improves.

9 Declarations of interest

None.