A moving boundary problem with variable specific heat and thermal conductivity

1 Introduction

Melting and freezing processes are encountered widely in nature and in many industrial processes, such as freezing of water, casting of melted alloys, thawing of food products, welding, thermal energy storage with phase change material, cryosurgery, production of steel and plastic products. During these processes, the material undergoes phase change includes a boundary that separates the two different phases. This boundary propagates in the material undergoing the phase change during the process. Mathematical formulation of the melting and freezing processes is governed by Stefan problems. Stefan problem (a moving boundary problem) describing the process of melting and freezing has been studied since eighteen century. These kinds of problems always attract interests due to the existence of one or more moving interfaces, inherent non-linear nature even in its simplest form and its wide applications in many natural/industrial processes. A detail discussion of various mathematical models related to the moving boundary problems and its analytical and approximate solutions is mentioned in the book of Crank (1984). The formulation of the problem with complicated boundary conditions can be seen in (Carslaw and Jaeger, 1959; Cho and Sunderland, 1974; Hill, 1986; Oliver and Sunderland, 1987; Petrova et al., 1994; Tritscher and Broadbridge, 1994).

From last one decade, the Stefan problem involving variable thermal coefficients (Briozzo et al., 2007; Briozzo and Natale, 2015; Briozzo and Natale, 2017; Kumar et al., 2020) has attracted great to Mathematicians as well as scientists because of its applicability and difficulty in getting its solution. Recently, Ceretani et al. (2018) considered a Stefan problem which involves thermal conductivity as a function of temperature and a Neumann type boundary condition at the left boundary and discussed the exact solution to the problem. A temperature-dependent thermal conductivity has been considered by Animasaun (2015) in his study of an incompressible electrically conducting Casson fluid flow along a vertical porous plate. Animasaun (2017) assumed temperature-dependent thermal conductivity and fluid viscosity in his study of a problem of steady mixed convection micropolar fluid flow towards stagnation point formed on horizontal linearly stretchable melting surface. Some more models involving temperature-dependent thermal conductivity can also be found in (Koriko and Animasaun, 2017; Makinde et al., 2018). Sandeep et al. (2017) presented a numerical exploration to examine the momentum, thermal and concentration boundary level behaviour of liquid-film flow of non-Newtonian nanofluids by assuming space and temperature dependent heat source/sink. Motivated by these works, we have discussed the following phase change problem in the domain $x>0$ that includes variable heat capacity and thermal conductivity:

To govern the position of moving interface, we need one additional condition on the boundary $x=s(t)$ which is known as the Stefan condition of the problem and it is given by

Besides conditions (2)–(4), an initial condition associated with the moving boundary is

In this paper, the heat capacity is assumed as

In this area, besides the mathematical model of the problem in different physical process, the establishment of solution of the mathematical model is also an exciting point of interest. Meek and Norbury (1984) presented a moving boundary problem which models the spreading of the viscous fluid under the gravitational force above a smooth horizontal plane and used the modified Keller box method to find a numerical solution of the problem. Therefore, many approximate, numerical and exact solutions of these problems have been reported in (Savovic and Caldwell, 2003; Natale and Tarzia, 2006; Rajeev et al., 2009; Słota and Zielonka, 2009; Rajeev, 2014; Fazio, 2013; Voller and Falcini, 2013; Zhou and Li-jiang, 2015). As far as author’s knowledge, exact solutions to the Stefan-type problems can be found by using similarity transformations only. In this study, the appropriate similarity variables are considered which allow us to convert the problem into an ordinary differential equation (ODE) along with boundary conditions. The exact solution to the proposed problem has been discussed for $m=n=1$ and $m=n=2$ . In order to discuss the solution for all positive integers $m$ and $n$ , the converted system of ODE is solved by using the shifted Chebyshev spectral technique. Parand and Razzaghi (2004) discussed to solve ordinary differential equations of higher order by the rational Chebyshev tau method. The approximate solution of ODE with the aid of shifted Chebyshev tau technique is described in (Doha et al. (2011a,b). An approximate solution to partial differential equations with fractional derivative by tau method is also discussed in (Vanani and Aminataei, (2011); Doha et al. (2011a,b). Ghoreishi and Yazdani (2011) discussed a generalization of the Tau method and presented its convergence analysis to numerical solution of multi-order fractional differential equations.

The paper has been arranged as follow: Section 2 demonstrates some properties related to the shifted Chebyshev polynomials. We have used this operational matrix of differentiation in our calculations. Next, the solution for all non-negative integers $m$ and $n$ is discussed in section 3 by applying a shifted Chebyshev tau method. Section 3 describes the exact solutions to the problem for two cases, i.e. $m=n=1$ and $m=n=2$ .

The existence and uniqueness of the exact solutions (obtained in Section 3) are discussed in Section 4. Finally, Section 5 contains the comparison of obtained approximate solution (given in Section 3) with exact solution for some cases. The dependent of temperature distribution and interface on m, n and Stefan number are also discussed in Section 5. The effect of Stefan number on the evolution of the moving boundary can be seen in the article of Savovic and Caldwell (2003).

2 Some prelimaries

It is well known that the first kind Chebyshev polynomials $\{T_i(t)\text{;}i=0,1,\dots \}$ are defined on $[-1,1]$ . We substitute a new variable $t=\frac{2x}{l}-1$ in $T_i(t)$ to use these polynomials on the interval $x\in [0,l]$ . The polynomials $T_i(\frac{2x}{l}-1)$ are called as first kind shifted Chebyshev polynomials. Let us denote shifted Chebyshev polynomials by $T_{l,i}(x)$ which satisfy the following recurrence formula (Doha et al. (2011a,b):

This work considers some properties of shifted Chebyshev polynomials (Doha et al. (2011a,b) which are given below:

1. In the interval $[0,l]$ , a function $g(x)$ which is square integrable can be expressed as:

For our calculation purpose, we consider the following partial sum of the series (9):

• The first derivative of $\varphi (x)$ and its relation with the operational matrix is described as:

• The nth derivative of $\varphi (x)$ in the terms of $D^{(1)}$ can written as

3 Solution for general case

First of all, we use the transformation defined as follows:

Now, we take the similarity variable defined as

Where $\lambda$ is a constant yet to be found.

Next, substituting the variables given in Eqs. (22), (23) into the Eqs. (17)–(20), we have the following system consisting of ODE:

Now, we can use the (N + 1)th partial sum of the series given in (11) for an approximate solution to the problem given in Eqs. (24)–(27). Therefore, the dependent variable $f(\eta )$ can be stated as:

As given in Eq. (14), the derivatives of dependent variable $f$ can be approximated as:

From Eqs. (28) and (29), the residual $R_N(x)$ corresponding to Eq. (24) is given as:

The $(N-1)$ algebraic equations can be found by the condition (Doha et al. (2011a,b) given below:

Moreover, by substituting the Eqs. (28) and (29) into the Eqs. (25)-(27), the following equations can be found:

Beside $(N-1)$ equations generated by Eq. (31), three more algebraic equations can be generated by Eqs. (32)–(34). Now, the system of $(N+2)$ algebraic equations with $(N+2)$ unknowns can easily be solved which determines the unknown vector $C$ and $\lambda$ . Consequently, the temperature distribution in liquid region $\theta (x,t)$ and $s(t)$ can be determined with the help of Eqs. (22), (23).

4 Exact solutions

In this section, we categorise the problem into two parts as:

Case 1: When $m=n=1$ and $\beta =\alpha$ then the Eqs. (24)–(26) can be written as:

The general solution of Eq. (35) is

After substituting the above values of $C_1$ and $C_2$ , the exact solution of the Eq. (35) along with boundary condition becomes:

In view of Eqs. (22)–(23) and (41), the solution of Eq. (17) at $m=n=1$ and $\alpha =\beta$ can be given by

Substituting Eqs. (41) into Eq. (37), we get the following transcendental equation:

The solution of Eq. (44) gives $\lambda$ and by substituting this value into (23), a tracking of the interface position $s(t)$ with time can be found.

Case 2: If $m=n=2$ and $\beta =\alpha$ then the Eqs. (24)–(27) become:

The solution of Eq. (45) with the boundary conditions (46) is given by

Consequently, the $\theta (x,t)$ at $m=n=2$ and $\alpha =\beta$ can be determined by substituting $\eta =x/2\sqrt{{\alpha }_{0}t}$ in the Eq. (48).

The Eqs. (47), (48) produce the following transcendental equation:

Solving (50) for $\lambda$ , will, on substitution into (23), provide the phase front $s(t)$ .

5 The existence and uniqueness

to validate the existence and uniqueness of solution established previously, we discuss as follows:

For case 1, we consider the transcendental equation given in Eq. (44) and suppose

It is obvious that $f_1(\lambda )$ is defined and continuous on $(0,\infty )$ and

From Eq. (52) and (53), it is clear that $f_1(\lambda )=0$ has at least one solution in $(0,\infty )$ .

Now, for all $\mathit{Ste}>0$ , it is clear that

Hence, $f_1(\lambda )$ is strictly increasing and this shows the uniqueness of $\lambda$ . Existence of unique $\lambda$ which satisfies the transcendental equation (44) assures the existence and uniqueness of solution to the problem (9)-(13) for $m=n=1$ and $\beta =\alpha$ .

For case 2, we define $f_2(\lambda )$ on $(0,\infty )$ with the help of transcendental Eq. (50) as

Clearly, $f_2(\lambda )$ is continuous on $(0,\infty )$ and

Hence $f_2(\lambda )=0$ has a solution in $(0,\infty )$ . Moreover, $f{\prime }_2(\lambda )>0$ on $(0,\infty )$ for all positive Stefan number which shows that $f_2(\lambda )$ is strictly monotonically increasing function. Hence, $f_2(\lambda )=0$ has a unique solution on $(0,\infty )$ . Consequently, there exists unique solution to the problem (9)–(13) for $m=n=2$ and $\beta =\alpha$ .

6 Comparisons and discussion

In this paper, all the computations for temperature distribution θ(x, t) and moving interface s(t) have been made with the help of Wolfram Research (8.0.0) software at fixed value of ${\alpha }_0=1.0$ .

We first present accurateness of the approximate solution described in section 3 through the figures for the proposed problem by considering the following matrices:

Tables 1 represents the comparisons of approximate temperature ${\theta }_A(x,t)$ and exact temperature distribution ${\theta }_E(x,t)$ in $0<x<s(t)$ at t = 1.0 for m = n = 1 and m = n = 2. The correctness of proposed approximate solution $s_A(t)$ for the moving interface at m = n = 1 and m = n = 2 is shown in Table 2. From these tables, it can be seen that our proposed approximate solutions are near to exact solutions $s_E(t)$ in the considered cases. Therefore, to explore the Stefan problem involving non-linear heat equation, this simple approach (stated in Section 3) can be useful to solve the problem.

With the help of proposed approximate solution, the variations of temperature distribution θ(x, t) and moving interface s(t) are shown in Figs. 1 and 2. In Fig. 1, the dependence of temperature distribution θ(x, t) on x is depicted at t = 1.0 and ${\alpha }_0=1.0$ for various values of m, n (m = n = 1, m = n = 2 and m = n =3) and $\beta$ $(\beta =0.5,1.5,2.5)$ . From this figure, it can be seen that the temperature is maximum at x = 0 and is continuously decreasing to zero at the moving interface. Moreover, it is clear that the temperature decreases in molten region as the value of m and/or n or $\beta$ decreases. Fig. 2 shows the trajectory of moving interface s(t) at ${\alpha }_0=1.0$ and $\alpha =0.5$ for different m, n and $\beta$ . This figure confirms that the velocity of moving interface s(t) improves when we increase either m and/or n or $\beta$ .This implies that the melting of material enhances when the parameters m or n or $\beta$ rises. The effect of Stefan number on the moving phase front is depicted in Figs. 3 and 4 for m = n = 1 and m = n = 2, respectively. These figures show that the larger values of Stefan numbers accelerate the movement of phase front which makes the process of melting fast. This is similar result as reported in the paper of Savovic and Caldwell (2003).

Close

Fig. 1

Plot of θ(x,t)vs. x at Ste = 0.5 and α = 0.5.

Export to PPT

Close

Fig. 2

Plot of s(t)vs. t at Ste = 0.5 and α = 0.5.

Export to PPT

Close

Fig. 3

Plot of phase front for different values of Ste at m = n = 1, α = β = 1 and α₀ = 1.

Export to PPT

Close

Fig. 4

Plot of phase front for different values of Ste at m = n = 2, α = β = 1 and α₀ = 1.

Export to PPT

7 Conclusion

In this study, the one-phase Stefan problem of melting process with variable thermal conductivity and heat capacity is discussed.Two exact solutions of the problem are presented for particular cases with the help of similarity variables method. Existence and uniqueness of exact solutions are also discussed. It is found that the movement of moving boundary s(t) is proportional to $\sqrt{t}$ in the proposed model and this result was well established earlier for the Stefan problem with $\alpha =\beta =0$ (Crank, 1984; Carslaw and Jaeger, 1959).

Besides exact solutions, an approximate approach based on similarity transformation and spectral tau method has been successfully applied to obtain the solution to the problem for general case. From section 6, it has been observed that the growth in the rate of change of temperature in molten region and the melting process are found if the value of m and/or n or $\beta$ increases. It is also observed that the proposed approximate approach is efficient, accurate and easy to apply on Stefan problems. The authors believe that this scheme is helpful for the researchers working in the field of moving boundary problem.

Declarations of interest

None.