A Stefan problem with temperature and time dependent thermal conductivity

1 Introduction

It is known that many processes like melting, freezing, sediment mass transport, tumour growth, etc. in the field of science and industry involve moving boundary/boundaries, and these problems are referred as moving boundary problems (or Stefan problems). Initially, the Stefan problems are restricted to heat-transfer problems and the formulations of these problems are developed for constant thermal properties (Crank, 1984). But, the Stefan problems are not only limited to heat-transfer problems with constant thermal properties. Some Stefan problems with different thermal properties and other diffusion controlled transport systems are discussed in Carslaw and Jaeger (1959), Hill (1986), Voller et al. (2004), Zhou and Li-jiang (2015).

From the literature (Cho and Sunderland, 1974; Oliver and Sunderland, 1987; Briozzo et al., 2007; Briozzo and Natale, 2015), it can be seen that moving boundary problems with temperature dependent thermal conductivity have been a fruitful research in the field of heat transfer. In 2017, Briozzo and Natale (2017) considered the temperature-dependent thermal conductivity in study of the supercooled one-phase Stefan problem for a semi-infinite material. Recently, Ceretani et al. (2018) discussed the similarity solutions for a one-phase Stefan problem with temperature-dependent thermal conductivity and a Robin condition at a fixed face. Voller and Falcini (2013) presented a one phase Stefan problem with diffusivity as a function of space and discussed an exact solution for it. In context of time dependent thermal conductivity, Hussein and Lesnic (2014) discussed the identification of time dependent thermal conductivity of an orthotropic rectangular conductor. Recently, Huntul and Lesnic (2017) also discussed an inverse problem of determining the time-dependent thermal conductivity and the transient temperature satisfying the heat equation with initial data. Motivated by these works, we consider a one phase Stefan problem with time and temperature dependent thermal conductivity of the form

Due to presence of moving boundary/boundaries or unknown domain, the moving boundary problems are nonlinear in nature even in its simplest form. If thermal conductivity is time and temperature dependent then the problem becomes more complicated to get its exact solution. In general, scaling invariance analysis and similarity variables (Briozzo et al., 2007; Ceretani et al., 2018; Fazio, 2013) play an important role for getting the exact solutions of these problems. In our study, we have also used the appropriate similarity variables to convert the governing system of partial differential equations into another system that includes ordinary differential equations with its conditions. After that, a shifted Chebyshev tau method based on Chebyshev operational matrix of differentiation is used to solve the transformed system. In Doha et al. (2011a, 2011b), the authors have discussed the shifted Chebyshev tau and collocation methods based on Chebyshev operational matrix of fractional derivatives for solving the linear multi-order fractional differential equations. Some other work related to shifted Chebyshev tau and collocation methods are reported in Ghoreishi and Yazdani (2011) and Vanani and Aminataei (2011).

2 The shifted Chebyshev polynomials and its operational matrix of differentiation

As we know that the Chebyshev polynomials $\{T_i(t)\text{;}i=0\text{,}1\text{,}\dots \}$ are defined on the interval $(-1\text{,}1)$ . In order to use these polynomials on the interval $x\in (0\text{,}L)$ , we introduce a new variable $t=\frac{2x}{L}-1$ in $T_i(t)$ which is called as shifted Chebyshev polynomials (Doha et al., 2011b; Ghoreishi and Yazdani, 2011).

Let the shifted Chebyshev polynomials $T_i\left(\frac{2x}{L}-1\right)$ be denoted by $T_{L\text{,}i}(x)$ , satisfying the following recurrence formula:

1. A square integrable function $u(x)$ in $(0\text{,}L)$ can be expressed in terms of the shifted Chebyshev polynomials as:

For practice purpose, only the first $(N+1)$ -terms shifted Chebyshev polynomials can be considered for the approximation of the function $u(x)$ . Hence, we can write

• The derivative of the vector $\varphi (x)$ is given by

If $N$ is even then we have $$D^{(1)}=\frac{2}{L}\left(\begin{array}{cccccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 4\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 3\hfill & \hfill 0\hfill & \hfill 6\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 8\hfill & \hfill 0\hfill & \hfill 8\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 5\hfill & \hfill 0\hfill & \hfill 10\hfill & \hfill 0\hfill & \hfill 10\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill \cdots \hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill N-1\hfill & \hfill 0\hfill & \hfill 2N-2\hfill & \hfill 0\hfill & \hfill 2N-2\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 2N\hfill & \hfill 0\hfill & \hfill 2N\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 2N\hfill & \hfill 0\hfill \end{array}\right)\text{.}$$

• The $\mathit{nth}$ order derivative of the vector $\varphi (x)$ is given by

3 Mathematical model

In this section, we consider the temperature and time dependent thermal conductivity as given in Eq. (1) and a mathematical model of one phase Stefan problem with nonlinear heat conduction is presented for melting process which is as follow:

By considering the following transformation:

4 Solution for the problem

Now, we consider the following similarity variables:

Substituting Eqs. (21) and (22) into Eqs. (16)–(19) which provide the following equations:

For the solution of Eqs. (23)–(26), finite numbers of terms i.e., the first $N+1$ terms of the series given in Eq. (3) are considered. Hence, the unknown function $f(\eta )$ is expressed in terms of the shifted Chebyshev polynomials as:

and $\varphi (\eta )={[T_{\lambda \text{,}0}(\eta )\text{,}{T}_{\lambda \text{,}1}(\eta )\text{,}\dots \text{,}{T}_{\lambda \text{,}N}(\eta )]}^T$ .

As given in Eq. (8), the derivatives are approximated as:

Using Eqs. (27) and (28), the residual $R_N(x)$ for Eq. (23) is defined as:

According to Tau method (Doha et al., 2011a, 2011b), we generate $(N-1)$ non-linear algebraic equations by using the condition

Also, by using Eqs. (27) and (28) in the Eqs. (24)–(26), we get

Eq. (30) generates $(N-1)$ equations and two more equations are generated by Eq. (31). Hence, we have $(N+1)$ equations in $(N+1)$ unknowns that can be easily solved and it gives the unknown coefficients of the vector $C$ . Consequently, $f(\eta )$ given in Eq. (27) can be calculated in terms of $\lambda$ which is still to be determined. In order to get the value of $\lambda$ , we use the calculated value of $f(\eta )$ in the interface condition given in Eq. (32).

5 Result and discussion

In this section, we discuss the accurateness of our obtained results as well as dependence of heat distribution and phase front on various parameters. By using the similarity transformation (given in Eqs. (21) and (22)), the analytical solution of Eqs. (23)–(26) is calculated for the constant thermal conductivity i.e., $\beta =0$ which is given as:

The location of phase front is given by

From Eq. (35), it is clear that for all $\lambda >0$ , $\mathit{dh}/d\lambda >0$ for positive values of $\alpha \text{,}\delta$ and $\mathit{Ste}$ . Moreover, $h(\lambda )\to -\infty$ if $\lambda \to 0$ and $h(\lambda )\to \infty$ as $\lambda \to \infty$ when $\alpha \text{,}\delta$ and $\mathit{Ste}$ are positive. Therefore, there exists one and only one positive value of $\lambda$ as the solution of Eq. (35). With the help of Eq. (34) and the obtained value of $\lambda$ from Eq. (35), the location of phase front $s(t)$ can be determined.

In order to show the accuracy, the comparisons between obtained results, exact results (given in Eqs. (33)–(35)) of temperature distribution and interface location at $\beta =0$ are depicted in Tables 1 and 2, respectively. Table 1 shows the obtained approximate values of temperature distribution ${\theta }_A$ for N = 3, 4, 5 and its exact value ${\theta }_E$ at $\mathit{Ste}=0.2$ , $\delta =1$ , $t=1.5$ and $\beta =0$ . Table 2 depicts the values of approximate position of phase front $s_A(t)$ for N = 3, 4, 5 on different time and its exact values $s_E(t)$ at $\mathit{Ste}=0.2\text{,}\delta =1$ and $\beta =0$ . From these tables, it is clear that our approximate results are near to exact value and accuracy increases as the order of operational matrix of differentiation increases.

When $\beta \ne 0$ , the obtained results are presented through Figs. 1–5 for the study of dependence of temperature distribution and location of phase front on various parameters. Fig. 1 demonstrates the variations of temperature distribution for different value of $\alpha (\alpha =2.0\text{,}1.0\text{,}0.2)$ at fixed values of $\beta =0.5$ , $\mathit{Ste}=0.2$ and $\delta =1.0$ . Fig. 2 depicts the variations of temperature distribution for different value of $\beta (\beta =1.5\text{,}1.0\text{,}0.3)$ at fixed values of $\alpha =0.5$ , $\mathit{Ste}=0.2$ and $\delta =1.0\text{.}$ From these figures, it is clear that temperature at $x=0$ is highest and decreases continuously to zero. It is also seen that the rate of change of temperature deceases as the parameters $\alpha$ and/or $\beta$ decrease.

Close

Fig. 1

Plot of temperature profile for different values of $\alpha$ at $\beta =0.5$ , $\mathit{Ste}=0.2$ and $\delta =1.0$ .

Export to PPT

Close

Fig. 2

Plot of temperature profile for different values of $\beta$ at $\alpha =0.5$ , $\mathit{Ste}=0.2$ and $\delta =1.0$ .

Export to PPT

In Fig. 3, the dependence of phase front on time for different $\alpha$ (i.e., exponent power of time) is presented at fixed value of $\beta =0.5$ , $\mathit{Ste}=0.2$ and $\delta =1.0\text{.}$ From this figure, it can be seen that the movement of phase front increases with the increment in the value of $\alpha (\alpha =0.2\text{,}1.0\text{,}2.0)$ . Consequently, the melting process becomes fast as we increase the value of $\alpha \text{.}$ Fig. 4 shows the trajectory of phase front for different $\beta (\beta =0.3\text{,}1.0\text{,}1.5)$ at fixed value of $\alpha =0.5$ , $\mathit{Ste}=0.2$ and $\delta =1.0$ . Fig. 5 demonstrates the trajectory of phase front for different Stefan numbers $(\mathit{Ste}=0.2\text{,}1.0\text{,}2.0)$ at fixed value of $\alpha =0.5$ , $\beta =0.5$ and $\delta =1.0$ . From Figs. 4 and 5, it is clear that the movement of phase front increases as the value of $\beta$ and/or $\mathit{Ste}$ increases. Hence, the melting process becomes fast if we increase the parameter $\beta$ and/or Stefan number $(\mathit{Ste})$ .

Close

Fig. 3

Plot of moving interface for different values of $\alpha$ at $\beta =0.5$ , $\mathit{Ste}=0.2$ and $\delta =1.0$ .

Export to PPT

Close

Fig. 4

Plot of moving interface for different values of $\beta$ at $\alpha =0.5$ , $\mathit{Ste}=0.2$ and $\delta =1.0$ .

Export to PPT

Close

Fig. 5

Plot of moving interface for different values of $\mathit{Ste}$ . at $\alpha =0.5$ , $\beta =0.5$ and $\delta =1.0$ .

Export to PPT

6 Conclusions

In this work, a special type of one phase Stefan problem with time and temperature dependent thermal conductivity is explored and its approximate solution is discussed by using similarity transformation method and shifted Chebyshev tau method based on Chebyshev operational matrix of differentiation. In order to check accuracy of our obtained results, an exact solution of the problem is also discussed for a particular case i.e., $\beta =0$ . From this study, it is seen that the proposed algorithm for the solution of Stefan problems is simple and accurate. Moreover, it is found that the rate of change of temperature increases as the power of time (i.e., $\alpha$ ) and/or $\beta$ increases and movement of moving interface increases if we increase the value of power of time (i.e., $\alpha$ ) or $\beta$ or $\mathit{Ste}$ . Consequently, the increment in the value of parameters $\alpha$ or $\beta$ or $\mathit{Ste}$ increases the rate of melting process. It is also observed that the variation of Stefan number is more pronounced than the parameters $\alpha$ and $\beta$ in the movement of interface.

Declarations of interest

None.