1 Introduction

Phase-change problems or Stefan problems encounter in many aspects of natural and industrial phenomena. These problems are particular cases of the moving boundary problems, where prior information about the location of the moving boundary is not known, but one has to calculate it as a part of the solution. Due to practical applications in the field of science, engineering and technology, Stefan problems have been attracting the interest of many researchers for over a century. The occurrence of Stefan problems can be seen in many specific fields of physical science such as melting or solidification process (Crank, 1987; Gupta, 2017), crystal growth process (Soni and Bharat, 1999), thermal energy storage process (Dincer and Rosen, 2002), metal casting (Poirier and Salcudean, 1988), shoreline problem (Trueba and Voller, 2010; Swenson et al., 2000; Voller et al., 2004; Rajeev and Kushwaha, 2013), and in many more areas.

In the classical Stefan problems, it has been assumed that the thermal coefficients are constants with respect to the temperature of the material. But it does not always happen with the many materials. Taking this fact into account, the researchers modified the Stefan problems in many ways to include new physical behaviour such as temperature-dependent thermo-physical property of the materials (Kumar et al., 2018,2019; Singh et al., 2011). Briozzo et al. (2007) considered the one-dimensional Stefan problem including temperature-dependent thermal coefficients and established the exact solution of the problem. Ceretani et al. (2018) assumed thermal conductivity as a function of temperature and Robin type boundary condition in the study of a one-phase Stefan problem and discussed the explicit solution to the problem. Many authors also considered these modifications and constructed the exact as well as numerical solutions of such type of problems (Ramos et al., 1994; Rogers and Broadbridge, 1988; Olive and Sunderland, 1987; Broadbridge and Pincombe, 1996; Briozzo and Natale, 2015). Mazzeo et al. (2015) presented an analytical solution of a Stefan problem in a finite PCM layer with time-dependent boundary condition in the steady periodic regime. In a steady periodic regime, Mazzeo and Oliveti (2017) discussed an approximation of the exact analytical solution to the Stefan problem in a finite PCM layer and a parametric study is also presented. Potvin and Gosselin (2009) discussed a numerical model to determine the thermal shielding of multilayer walls containing layers of phase change materials. In 2018, a numerical solution based on the finite difference scheme to the problem associated with thermal field and heat storage in a cyclic phase change process is presented by Mazzeo and Oliveti (2018).

In many physical processes, the phase change material is allowed to move when the phase change occurs. In the literature, the phase change problem with moving material is not adequately studied yet. Recently, Turkyilmazoglu (2018) presented some Stefan problems involving moving phase change material and discussed analytical solutions of the problems. This modification encourages us to explore a one dimensional Stefan problem which includes temperature-dependent thermal conductivity and specific heat, and movement of material during the phase change process simultaneously. We have taken thermal conductivity and specific heat, respectively in the following form:

and

In this paper, we have explored a one phase Stefan problem with variable thermal coefficients as mentioned in Eqs. (1) and (2) and it is also assumed that the phase change material is moving with unidirectional speed. The section wise description of this article is given below:

Section 2 presents some basic definitions of the Chebyshev polynomials and its properties. The mathematical model of the problem has been elaborated in detail in the section 3. Section 4 contains an approximate approach to the problem for all $\alpha$ and $\beta$ using spectral tau scheme. In the next section, the exact solution to the problem has been constructed for the case $\alpha =\beta$ and its uniqueness has also been established. Results and discussion section describes the findings of the study in detail and the outcomes have been presented by tables and figures. The last section, conclusion summarises the work a very lucid and concise manner.

2 Some basic definitions and properties

First, we define the first kind Chebyshev polynomials $T_n(t)$ recursively on the interval [−1,1] which is given below

Now we substitute $t=\frac{2x}{\lambda }-1$ in the Chebyshev polynomials to use these polynomials on the interval $\left[0,\lambda \right]$ . This substitution changes the Chebyshev polynomials to the polynomials $T_i(\frac{2x}{\lambda }-1)$ which are known as the shifted Chebyshev polynomial denoted by $T_{\lambda ,i}(x)$ . From the recurrence relation (3), we can find the desired number of shifted Chebyshev polynomials explicitly.

According to Doha et al. (2011), any square integrable function $f(x)$ on the interval $\left[0,\lambda \right]$ can be written as

$h_0=\pi ,h_j=\frac{\pi }{2},j=1,2,...$

and

$w_{\lambda }(x)=\frac{1}{\sqrt{\lambda x-x^2}},$

respectively.

To approximate the function $f(x)$ , we can take first $N+1$ terms of the series (4) and $f(x)$ may be denoted by $f_N(x)$ that is given by

where the vector $C$ and the shifted Chebyshev vector $\phi (x)$ are given by $$C=[c_0,c_1,...,c_N],$$

and $$\phi (x)={[T_{\lambda ,0}(x),T_{\lambda ,1}(x),...,T_{\lambda ,N}(x)]}^T.$$

Now, the relation of $\phi \prime (x)$ and the vector $\phi (x)$ is given by

The operational matrix of higher order derivatives can be given by

3 Mathematical model and its solution

We have considered a one phase Stefan problem with convective boundary condition (Ceretani et al. (2018)) in a semi-infinite domain governing the freezing process. Initially, it is assumed that the material is at its freezing temperature $T_f$ and the density does not change when phase change occurs. Apart from the classical Stefan problem, it is also considered that the phase change material is moving in the positive direction of $x$ -axis with a speed $u$ which depends on time and the considered Stefan problem (Turkyilmazoglu (2018)). The mentioned problem can be modelled as:

To solve the problem (10)–(14), we first use the following similarity transformation:

With the help of (15), the above Eq. (10) becomes the following ordinary differential equation (ODE):

From (12), (13) and (15), it can be found that the position of moving interface must be proportional to $\sqrt{\nu \mathit{t}},$ and hence can be presented as follows:

The similarity transformation (15) converts the boundary conditions (11), (12) and (13) respectively into the following forms:

To solve the Eq. (16) for the unknown function $\theta (\eta )$ along with the prescribed boundary conditions (18) and (19), we first approximate the unknown function $\theta (\eta )$ as given in (6) by

According to Eq. (9), the approximations of derivatives of $\theta (\eta )$ is considered as

Substituting the value $\theta (\eta )$ and its derivatives into the Eq. (16), we get the residual denoted by $R_N(\eta )$ of the differential Eq. (16) that is given below

As in tau method (Doha et al. (2011)), we can generate $(N-1)$ algebraic equations with $(N+2)$ unknowns from the following equation:

To complete the system of algebraic equations, we can find the remaining two equations from the boundary conditions (18) and (19), which are

In spite of above $N+1$ algebraic equations, one more algebraic equation can be found from the Eqs. (20) and (24) which is given below:

Eqs. (26) and (29) generate $N+2$ algebraic equations involving $N+2$ unknowns $\left(c_0,c_1,...,c_N,\lambda \right)$ . These obtained equations can be easily solved by Newton-Raphson method or any mathematical software for the unknowns and these results are required to get the solution to the problem.

4 Exact solution

In this section, we have constructed the exact solution of the problem (10)–(14) for the case $\alpha =\beta .$ With the help of similarity variable given in (15), the partial differential Eq. (10) becomes

Like the previous section, we assume again the location of the moving boundary $s(t)$ as

The general solution of the Eq. (30) can be given as

The similarity transformation (6) converts the boundary conditions (2), (3) and (4), respectively into following form:

where $\mathit{Ste}=\frac{c_0\left(T_f-T_{\infty }\right)}{L}$ is the Stefan number.

With the help of Eq. (32), the conditions (33) and (34) determine the unknown constants $C_1$ and $C_2$ in terms of unknown $\lambda$ as

where

The Eqs. (32), (36) and (37) gives the solution of the Eq. (30).

Now, we consider the Eqs. (8), (14) and (15) that give rise to the following transcendental equation:

From Eq. (39), we can find the unknown $\lambda ,$ if it exists, to get the solution. In the next section, we will show the existence and uniqueness of $\lambda$ which satisfies the transcendental Eq. (39).

5 Existence and uniqueness

To show the existence and uniqueness of solution constructed in the previous section, it is sufficient to show that there is a unique value of $\lambda$ in $\left(0,\infty \right)$ which satisfies the Eq. (29). To show this fact, we consider the following function:

The function $f(\lambda )$ is continuous and differentiable on the interval $(0,\infty )$ and ${\mathit{Lim}}_{\lambda \to 0^{+}}f(\lambda )$ is a negative real number for all the involved positive parameters. Moreover, $\underset{\lambda \to \infty }{\mathit{Lim}}f(\lambda )=\infty$ which shows that $f(\lambda )$ has a root in the interval $(0,\infty )$ that proves the existence of the solution the considered problem.

It can be seen that the derivative of $f(\lambda )$ is on the interval $(0,\infty )$ when we take the Peclet number $\mathit{Pe}⩽\sqrt{2}$ . Hence, $f(\lambda )$ is a strictly increasing function on the interval $(0,\infty )$ and this observation shows that $f(\lambda )$ has a unique root. It is also observed that when we take the value of the Peclet number $\mathit{Pe}>\sqrt{2}$ , the problem may have multiple solutions which agree with the findings explored in Turkyilmazoglu (2018).

6 Results and discussions

In this section, first, we discuss about the accuracy of the approximate solution given in section 3. In order to show the accuracy of the obtained approximate solution, two tables have been presented for $\alpha =\beta$ . Table 1 depicts the comparison of the location of interface for the exact solution $s{(t)}_E$ and the approximate solution $s_A(t)$ by taking operational matrices of order 3, 4 and 5. The numerical values in the Table 1 are given at $\alpha =\beta =1,{\alpha }_0=1,Pe=1$ and $\mathit{Ste}=0.5$ for different values of $\gamma$ in which ${\left.s_A(t)\right|}_{N=i}$ denotes the approximate value of the location of the interface for the operational matrix of order $i+1$ . Table 2 represents the assessment for the accuracy of the approximate value ${\theta }_A$ and the exact value ${\theta }_E$ of temperature profile for $\alpha =\beta =1,{\alpha }_0=1,t=1,Pe=1$ and $\mathit{Ste}=0.5.$ In table second, ${\left.{\theta }_A\right|}_{N=i},i=2,3,4$ represents the approximate values of temperature obtained by taking the operational matrices of order 3, 4 and 5, respectively. The Tables 1 and 2 clearly show that the approximate values are sufficiently accurate and are in good agreement with the exact solution for $\alpha =\beta$ .

Fig. 1 represents the plot between $s(t)$ and time $t$ for the parameters $\alpha =0.2,\beta =0.5,\gamma =0.5$ and $\mathit{Ste}=0.75$ . This figure clearly shows that the velocity of location of the interface increases rapidly with the increment of the values of Peclet number $\left(\mathit{Pe}\right)$ . This indicates that the material freezes fastly when we take the larger values of $\mathit{Pe}$ . The Fig. 2 shows the plot between the moving boundary $s(t)$ and time $t$ for different values of $\gamma$ at $\alpha =0.5,\beta =0.2,Pe=1.0$ and $\mathit{Ste}=0.75$ . This figure reveals that the moving interface propagates fast as we increase the parameter $\gamma$ and hence, the freezing process becomes fast.

Close

Fig. 1

Effect of $\mathit{Pe}$ on tracking of the moving interface $s(t)$ at $\alpha =0.2,\beta =0.5,\gamma =0.5$ and $\mathit{Ste}=0.75.$

Export to PPT

Close

Fig. 2

Effect of $\gamma$ on tracking of the moving interface $s(t)$ at $\alpha =0.5,\beta =0.2,Pe=1$ and $\mathit{Ste}=0.75.$

Export to PPT

7 Conclusion

In this article, a semi-analytical solution to a solidification problem with temperature-dependent thermal coefficients and moving phase change material is elaborated. The exact solution to the problem is also constructed for a particular case, i.e. $\alpha =\beta$ and it is observed that there exists a unique solution of the problem if we take $\mathit{Pe}⩽\sqrt{2}$ . In this model, it is seen that $s(t)$ is proportion to $\sqrt{t}$ which is the similar result as we have found in the literature for classical Stefan problems. Moreover, it is clear that the solidification process is affected by Peclet number ( $\mathit{Pe}$ ) and $\gamma$ , and the process becomes fast as the parameters $\mathit{Pe}$ and $\gamma$ increase. From this study, it is also observed that the proposed semi-analytical approach is simple, accurate and efficient and it will be helpful to the researchers of this area to handle Stefan-like problems.

Declaration of interests

None.