An exact solution to a Stefan problem with variable thermal conductivity and a Robin boundary condition

TL;DR

This paper derives an exact similarity solution for a Stefan problem with temperature-dependent thermal conductivity and Robin boundary condition, introducing a generalized error function and analyzing its properties.

Contribution

It provides a novel exact solution framework for a Stefan problem with variable thermal conductivity and Robin boundary condition, including properties and approximations of the generalized error function.

Findings

Existence of a unique, bounded, analytic solution for small parameter values.

The generalized modified error function is concave and increasing.

Explicit approximations for the generalized error function are proposed.

Abstract

In this article it is proved the existence of similarity solutions for a one-phase Stefan problem with temperature-dependent thermal conductivity and a Robin condition at the fixed face. The temperature distribution is obtained through a generalized modified error function which is defined as the solution to a nonlinear ordinary differential problem of second order. It is proved that the latter has a unique non-negative bounded analytic solution when the parameter on which it depends assumes small positive values. Moreover, it is shown that the generalized modified error function is concave and increasing, and explicit approximations are proposed for it. Relation between the Stefan problem considered in this article with those with either constant thermal conductivity or a temperature boundary condition is also analysed.

Tables

Table 1. Table 1 : Approximative solutions β 1 ∗ subscript superscript 𝛽 1 \beta^{*}_{1} to the equation ( 13 ) for γ = 0.1 , 1 , 10 , 100 𝛾 0.1 1 10 100 \gamma=0.1,1,10,100 .

$\gamma$	$0.1$	$1$	$10$	$100$
${\beta }_1^{\ast }$	$1.55$	$3\times {10}^{-1}$	$3.65\times {10}^{-2}$	$3.75\times {10}^{-3}$

Equations

$$\rho c{T}_t(x,t)=(k(T(x,t))T_x(x,t){)}_x$$

$$s(0)=0$$

$$T(s(t),t)=T_f$$

$$k(T_f)T_x(s(t),t)=\rho l\dot{s}(t)$$

$$k(T(0,t))T_x(0,t)=\frac{h_0}{\sqrt{t}}(T(0,t)-T_{\infty })$$

$$k(T)=k_0\left(1+\beta \frac{T-T_{\infty }}{T_f-T_{\infty }}\right),$$

$$T(0,t)=T_{\infty }t>0.$$

$$\phi \left(\frac{x}{2\sqrt{{\alpha }_{0}t}}\right)=\frac{T(x,t)-T_{\infty }}{T_f-T_{\infty }}0<x<s(t),t>0,$$

$$T(x,t)=(T_f-T_{\infty })\phi \left(\frac{x}{2\sqrt{{\alpha }_{0}t}}\right)+T_{\infty }$$

$$s(t)=2\lambda \sqrt{{\alpha }_{0}t}$$

$$[(1+\beta y(\eta ))y^{\prime }(\eta ){]}^{\prime }+2\eta y^{\prime }(\eta )=0$$

$$y^{\prime }(0)+\beta y(0)y^{\prime }(0)-\gamma y(0)=0$$

$$y(\lambda )=1$$

$$\frac{{\phi }^{\prime }(\lambda )}{\lambda }=\frac{2}{(1+\beta ){\text{Ste}}_{\infty }},$$

$$||h|{|}_{\infty }=\sup \left\{|h(x)|:0\le x\le \lambda \right\}(h\in X).$$

$$K=\left\{h\in X:0\le h\text{ and }||h|{|}_{\infty }\le 1\right\}.$$

$$[{\Psi }_h(\eta )y^{\prime }(\eta ){]}^{\prime }+2\eta y^{\prime }(\eta )=0$$

$${\Psi }_h(0)y^{\prime }(0)-\gamma y(0)=0$$

$$y(\lambda )=1,$$

$$y(\eta )=D_h\left(\frac{1}{\gamma }+{\text{bigintsss}}_0^{\eta }\frac{\exp \left(-2{\int }_0^x\frac{\xi }{{\Psi }_h(\xi )}d\xi \right)}{{\Psi }_h(x)}dx\right)0<\eta <\lambda ,$$

$$D_h=\gamma {\left(1+\gamma {\text{bigintsss}}_0^{\lambda }\frac{\exp \left(-2{\int }_0^x\frac{\xi }{{\Psi }_h(\xi )}d\xi \right)}{{\Psi }_h(x)}dx\right)}^{-1}.$$

$$\left(\tau h\right)(\eta )=D_h\left(\frac{1}{\gamma }+{\text{bigintsss}}_0^{\eta }\frac{\exp \left(-2{\int }_0^x\frac{\xi }{{\Psi }_h(\xi )}d\xi \right)}{{\Psi }_h(x)}dx\right)0<\eta <\lambda ,(h\in K)$$

$$g(x)=\frac{\sqrt{\pi }}{2}\gamma x(1+x{)}^{1/2}(3+x)x>0.$$

$$g(x)=1x>0$$

$$\begin{array}{rl}\left|(\tau h_1)(\eta )-(\tau h_2)(\eta )\right|& \le \gamma {\text{bigintsss}}_0^{\eta }\left|\frac{\exp \left(-2{\int }_0^x\frac{\xi }{{\Psi }_{h_1}(\xi )}d\xi \right)}{{\Psi }_{h_1}(x)}-\frac{\exp \left(-2{\int }_0^x\frac{\xi }{{\Psi }_{h_2}(\xi )}d\xi \right)}{{\Psi }_{h_2}(x)}\right|dx\\ & +\left|D_{h_1}-D_{h_2}\right|\left|\frac{1}{\gamma }+{\text{bigintsss}}_0^{\eta }\frac{\exp \left(-2{\int }_0^x\frac{\xi }{{\Psi }_{h_2}(\xi )}d\xi \right)}{{\Psi }_{h_2}(x)}dx\right|\\ & \le 2\gamma {\text{bigintsss}}_0^{\eta }\left|\frac{\exp \left(-2{\int }_0^x\frac{\xi }{{\Psi }_{h_1}(\xi )}d\xi \right)}{{\Psi }_{h_1}(x)}-\frac{\exp \left(-2{\int }_0^x\frac{\xi }{{\Psi }_{h_2}(\xi )}d\xi \right)}{{\Psi }_{h_2}(x)}\right|dx,\end{array}$$

$$\text{i) }0\le \phi (\eta )\le 1,\text{ii) }0<{\phi }^{\prime }(\eta ),\text{iii)}{\phi }^{\prime \prime }(\eta )<0\forall 0<\eta <\lambda .$$

$$(1+\beta y(\eta ))y^{\prime \prime }(\eta )+\beta (y^{\prime }(\eta ){)}^2+2\eta y^{\prime }(\eta )=00<\eta <\lambda ,$$

$$1+\beta \phi ({\eta }_0)>0,$$

$${\phi }^{\prime \prime }(\eta )=-\frac{\beta ({\phi }^{\prime }(\eta ){)}^2+2\eta {\phi }^{\prime }(\eta )}{1+\beta \phi (\eta )}0<\eta <\lambda .$$

$$\phi (\eta )=\sum _{n=0}^{\infty }{\beta }^n{\phi }_n(\eta )0<\eta <\lambda ,$$

$${\phi }^{(N)}(\eta )=\sum _{n=0}^N{\beta }^n{\phi }_n(\eta )0<\eta <\lambda ,$$

$$\sum _{n=1}^{\infty }\left(\sum _{k=1}^{n}a(\eta ,k-1,n-k)+b(\eta ,n)\right){\beta }^n+b(\eta ,0)=00<\eta <\lambda ,$$

$$a(\eta ,n,m)={\phi }_n^{\prime }(\eta ){\phi }_m^{\prime }(\eta )+{\phi }_n(\eta ){\phi }_m^{\prime \prime }(\eta )$$

$$b(\eta ,n)={\phi }_n^{\prime \prime }(\eta )+2\eta {\phi }_n^{\prime }(\eta )$$

$$\sum _{n=1}^{\infty }\left(\sum _{k=1}^n{a}_0(k-1,n-k)+b_0(n)\right)+b_0(0)=0,$$

Full text

Keywords: Stefan problems, exact solutions, temperature dependent thermal conductivity, convective boundary conditions, modified error function, phase-change processes.

1 Introduction

The understanding of phase-change processes has been inspiring scientists from the earlier 18th century. Already in 1831, Lamé and Clapeyron studied problems related to the solidification of the Earth planet [21]. Also the mathematical formulation of phase-change processes as free boundary problems dates from the 18th century, since it owes much to the ideas developed by Stefan in 1889 [29, 27, 28]. At present, their study is still an active area of research. Besides phase-change process are interesting in themselves, they attract interests because they are present in a wide variety of situations, both natural and industrial ones. Glass manufacturing and continuous casting of metals are examples of industrial activities involving them, some recent works in this area are [3, 19]. Controlling side-effects of certain industrial processes or preventing future problems derived from our energy-dependent lifestyle, are also examples of how phase-change processes arise as a subject of study [18, 15]. Permafrost phenomena or dynamics of snow avalanches are examples of natural situations whose study involves phase-change processes, some recent articles in these subjects are [20, 1, 6]. We refer the reader to [10, 32] and the references therein for a recent survey in applications and future challenges in free boundary problems. Other references can be seen in the last published Free Boundary Problems International Conference Proceedings [16].

In this article we will focus on phase-change processes that are ensued from an external temperature imposed at some part of the fixed boundary of a homogeneous material. A classical simplification in modelling this sort of phenomena is to consider boundary conditions of Dirichlet type (temperature conditions). This is based on the assumption that heat is instantaneously transferred from the external advise through which a specific temperature is imposed to the material. In view that is physically unrealistic, several authors have suggested to consider conditions of Robin type (convective conditions) since they mimic the fact that the heat transfer at the boundary is proportional to the difference between the imposed temperature and the one the material presents at its boundary (see for example the books [2, 7]). Another classical simplification when modelling phase-change processes is to consider that thermophysical properties are constant. Though it is reasonable for most phenomena under moderate temperature variations [2], it is not what actually happens as a rule. In fact, this hypothesis has been removed in many works in the attempt to improve the mathematical model (see, for example [9, 4, 25, 34]). All this have encouraged us to look at phase-change processes with convective boundary conditions and non-constant physical properties.

In 1974, Cho and Sunderland studied a phase-change process for a one-dimensional semi-infinite material with temperature-dependent thermal conductivity [11]. The dependence was assumed to be linear, which is a quite good approximation of what actually happens with several materials (water, for example [2]). The phase-change process was assumed to be ensued from a constant temperature imposed at the fixed boundary of the body, what was modelled through a Dirirchlet condition. For the resulting Stefan problem, Cho and Sunderland have presented an exact similarity solution. The temperature was obtained through an auxiliary function $\Phi$ that they have called a Modified Error (ME) function and that was defined as the solution to a nonlinear ordinary differential problem of second order. Revisiting the work of Cho and Sunderland, a couple of curiosities have arised. On one hand, the existence of the ME function was not proved there. Despite of this lack of theoretical results, the ME function was widely used in the context of phase-change processes before their existence and uniqueness were proved in the recent article [8] (see, for example, [14, 35, 5, 17, 22, 23, 25, 26, 30, 13]). On the other hand, by following the arguments presented in [11] it is obtained that the ME function must satisfy a differential problem over a closed bounded interval $[0,\lambda]$ with $\Phi(0)=0$, $\Phi(\lambda)=1$. Nevertheless, in [11] it was considered a boundary value problem over $[0,+\infty)$ with $\Phi(0)=0$, $\Phi(+\infty)=1$. Although in this way it is clearer the relation between the modified and classical error functions (see [11, 8] for further details), the change made by Cho and Sunderland add some extra conditions on the temperature function.

In this article we consider a similar phase-change process to that studied in [11]. We are mainly motivated by: a) improving the modelling of the imposed temperature at the fixed boundary by considering a convective boundary condition, b) obtaining a solution of similarity type without any extra condition on the temperature distribution. We will study a solidification process, but a completely similar analysis can be done for the case of melting. Aiming for simplicity, we will restrict our presentation to a one-phase process. That is, the case in which the material is initially liquid at its freezing temperature.

The organisation of the paper is as follows. First (Sect. 2), we introduce the one-phase Stefan problem through which we will study the phase-change process. In this Section we also present a characterisation for any similarity solution to the Stefan problem in terms of a Generalized Modified Error (GME) function. This will be defined as the solution to a nonlinear boundary value problem of second order. Similarly to [11], this problem will depend on a positive parameter $\beta$ related to the slope of the thermal conductivity as a linear function of the temperature distribution. Next (Sect. 3), we analyse the existence of the similarity solutions given in Section 2. In particular, we prove that there exists a unique non-negative analytic GME function when $\beta$ assumes small positive values. We also prove that this GME function is concave and increasing, and explicit approximations are proposed for it. Finally (Sect. 4), we discuss about how the Stefan problem presented in Section 2 is related to those studied in [11] (Dirichlet condition at the fixed boundary) and [33] (constant thermal conductivity).

2 The Stefan problem.

The one-phase solidification process introduced in Section 1 will be studied through the following Stefan problem:

[TABLE]

In (1), the unknown functions are the temperature $T$ of the solid region and the free boundary $s$ separating the phases. The parameters $\rho>0$ (density), $c>0$ (specific heat), $l>0$ (latent heat per unit mass), $h_{0}>0$ (coefficient related to the heat transfer at $x=0$), $T_{f}\in\mathbb{R}$ (freezing temperature) and $T_{\infty}<T_{f}$ (constant temperature imposed in the neighbourhood of the boundary $x=0$) are all known constants. The function $k$ (thermal conductivity) is defined as:

[TABLE]

where $k_{0}>0$, $\beta>0$ are given constants.

Remark 1**.**

Let us assume for a moment that for each $h_{0}>0$ exists a solution to problem (1) such that $T(0,\cdot)$, $T_{x}(0,\cdot)$ admit bounds independent of $h_{0}$ (what actually happens in the most common physical situations). Then, by taking the limit when $h_{0}\to\infty$ for each fixed $t>0$ in (1e), we obtain:

[TABLE]

In other words, if we were able to consider an infinite value for the heat transfer coefficient $\frac{h_{0}}{\sqrt{t}}$ in the convective boundary condition (1e), then the temperature function given through problem (1) would satisfy the temperature boundary condition (1e*†*). Thus, the mathematical framework given by problem (1) agrees well with the physical ideas about temperature and convective boundary conditions discussed in Section 1 (see, for example, [2, 7] for a detailed explanation of physical interpretations of boundary conditions).

We are interested here in obtaining a similarity solution to problem (1). More precisely, one in which the temperature $T(x,t)$ can be written as a function of the single variable $\frac{x}{2\sqrt{\alpha_{0}t}}$, where $\alpha_{0}=\frac{k_{0}}{\rho c}>0$ (thermal diffusivity for $k_{0}$). Through the following change of variables in problem (1):

[TABLE]

and a few simple computations, the following theorem can be proved (we refer the reader to the proof of Theorem 2 in [31] for an illustrative example of this sort of demonstrations).

Theorem 2.1**.**

The Stefan problem (1) has the similarity solution $T$, $s$ given by:

[TABLE]

if and only if the function $\varphi$ and the parameter $\lambda>0$ satisfy the following differential problem:

[TABLE]

together with the following condition:

[TABLE]

where $\text{Ste}_{\infty}=\dfrac{c(T_{f}-T_{\infty})}{l}>0$ (Stefan number) and $\gamma=2\text{Bi}>0$ with $Bi=\dfrac{h_{0}\sqrt{\alpha_{0}}}{k_{0}}$ (generalized Biot number).

Any function $\varphi$ that satisfies problem (5) will be referred to as Generalized Modified Error (GME) function. In the following Section we will prove that such functions exist.

3 Similarity solutions.

In this Section we will analyse the existence and uniqueness of the similarity solution (4) to problem (1). By virtue of Theorem 2.1, it can be done through the analysis of problem (5)-(6). First, we will study the differential problem (5) by assuming that $\lambda$ is a positive given number. As we shall see shortly, it can be proved that problem (5) has a unique non-negative analytic solution $\varphi$ when $\beta$ assumes small positive values. Then, we will analyse the relation between the found solution $\varphi$ and the parameter $\lambda$ through the study of the equation (6).

3.1 Analysis of problem (5).

This Section is devoted to the GME function. First, we will present a result on its existence and uniqueness. Then, it will be proved that the GME function is increasing and concave, as the classical error function is. Finally, explicit approximations are proposed for the GME function and several plots are presented for different values of the parameters involved in the physical problem (1).

3.1.1 Existence and uniqueness of the GME function

The ideas that will be developed in the following are based on the recent article [8], where it was proved the existence of the ME function introduced in [11]. Through a fixed point strategy, we will prove the existence of the GME function, $\varphi$.

All throughout this Section we will consider $\lambda>0$, $\gamma>0$ and $\beta\geq 0$. We will denote with $X$ to the set of all bounded analytic functions $h:[0,\lambda]\to\mathbb{R}$. It is well known that $X$ is a Banach space with the supremum norm $||\cdot||_{\infty}\,$, which is defined by:

[TABLE]

The subset of $X$ given by all non-negative functions which are bounded by 1 will be referred to as $K$, that is:

[TABLE]

Note that $K$ is a non-empty closed subset in $(X,||\cdot{}||_{\infty})$. Finally, for each $h\in K$ we will write $\Psi_{h}=1+\beta h$.

Remark 2**.**

$1\leq\Psi_{h}\leq 1+\beta$* for any $h\in K$.*

The main idea in the analysis below is to study the nonlinear problem (5) through the linear problem given by:

[TABLE]

where $h$ is a known function which belongs to $K$. The advantage in considering problem (8) lies in the fact that the differential equation (8a) can be easily solved as a linear equation of first order in $y^{\prime}$.

Lemma 3.1**.**

Let $h\in K$. The only solution $y$ to problem (8) is given by:

[TABLE]

with $D_{h}$ defined by:

[TABLE]

Moreover, $y\in K$.

Proof.

By observing that the constant $D_{h}$ given by (10) is well defined since $\Psi_{h}$ is never zero, the proof follows easily by checking that the function $y$ given by (9) satisfies problem (8). ∎

The next result is a direct consequence of the previous Lemma.

Theorem 3.1**.**

Let $y\in K$. Then $y$ is a solution to problem (5) if and only if $y$ is a fixed point of the operator $\tau$ from $K$ to $X$ given by:

[TABLE]

where $D_{h}$ is defined by (10).

In the following we will focus on analysing the existence of fixed points of $\tau$.

Theorem 3.2**.**

$\tau(K)\subset K$.

Proof.

Let $h\in K$. We have that:

1. i)

$\tau h$ is an analytic function, since $h\in X$. 2. ii)

$0\leq\tau h$, since $0<D_{h}$. 3. iii)

$||\tau h||_{\infty}\leq 1$, since $\left|\left(\tau h\right)(\eta)\right|\leq\left(\tau h\right)(\lambda)=1$ for all $0<\eta<\lambda$.

Then, $\tau h\in K$. ∎

Lemma 3.2**.**

Let $h,h_{1},h_{2}\in K$ and $\eta\in[0,\lambda]$. We have:

1. a)

$\displaystyle\bigintsss_{0}^{\eta}\left|\frac{\exp\left(-2\displaystyle\int_{0}^{x}\frac{\xi}{\Psi_{h_{1}}(\xi)}d\xi\right)}{\Psi_{h_{1}}(x)}-\frac{\exp\left(-2\displaystyle\int_{0}^{x}\frac{\xi}{\Psi_{h_{2}}(\xi)}d\xi\right)}{\Psi_{h_{2}}(x)}\right|dx$

$\leq\frac{\sqrt{\pi}}{4}\beta(1+\beta)^{1/2}(3+\beta)||h_{1}-h_{2}||_{\infty}$** 2. b)

$0<D_{h}\leq\gamma$.

Proof.

2. a)

See [8, Lemma 2.1]. 3. b)

It is a direct consequence of the positivity of $\Psi_{h}$ (see Remark 2).

∎

Lemma 3.3**.**

Let $g$ be the real function defined by:

[TABLE]

The equation:

[TABLE]

has an only positive solution $\beta_{1}=\beta_{1}(\gamma)$.

Proof.

It follows from the fact that $g$ is an increasing function in $\mathbb{R}^{+}$ with $\displaystyle\lim_{x\to 0^{+}}g(x)=0$ and $\displaystyle\lim_{x\to+\infty}g(x)=+\infty$. ∎

Theorem 3.3**.**

Let $\beta_{1}$ the only positive solution to equation (13). If $0\leq\beta<\beta_{1}$, then $\tau$ is a contraction.

Proof.

Let $h_{1},h_{2}\in K$ and $\eta\in[0,\lambda]$. From Lemma 3.2 and:

[TABLE]

we have that that $||\tau h_{1}-\tau h_{2}||_{\infty}\leq g(\beta)||h_{1}-h_{2}||_{\infty}$. Recalling that $g$ is an increasing function, it follows that $\tau$ is a contraction when $\beta<\beta_{1}$. ∎

We are in a position now to formulate our main result.

Theorem 3.4**.**

Let $\beta_{1}$ as in Theorem 3.3. If $0\leq\beta<\beta_{1}$ then problem (5) has a unique non-negative analytic solution.

Proof.

It a direct consequence of Theorems 3.1, 3.2, 3.3 and the Banach Fixed Point Theorem. ∎

3.1.2 Some properties of the GME function

This Section is devoted to prove that the GME function $\varphi$ found in Section 3.1.1 shares the following properties with the classical error function erf:

[TABLE]

In the following we will consider $\gamma>0$, $\lambda>0$, $\beta_{1}$ the solution to equation (13) and $0\leq\beta<\beta_{1}$. The first property in (14) is an immediate consequence of the fact that $\varphi\in K$. In order to prove ii) in (14), it will be enough to show that $\varphi^{\prime}(\eta)\neq 0$ for all $0<\eta<\lambda$, since $\varphi(0)\leq 1$ and $\varphi(\lambda)=1$. Let us assume that there exists $0<\eta_{0}<\lambda$ such that $\varphi^{\prime}(\eta_{0})=0$. As we shall see shortly, we will reach a contradiction. Since equation (5a) can be written as:

[TABLE]

and we have that:

[TABLE]

we obtain that $\varphi^{\prime\prime}(\eta_{0})=0$. From this, by differentiating (15) and taking into account (16), it follows that $\varphi^{\prime\prime\prime}(\eta_{0})=0$. We continue in this fashion obtaining that $\varphi^{(n)}(\eta_{0})=0$ for all $n\in\mathbb{N}$. But this implies that $\varphi\equiv 0$ in $[0,\lambda]$, since $\varphi$ is an analytic function. This contradicts $\varphi(\lambda)=1$. Finally, we have that the last property in (14) is a direct consequence of i), ii) and the fact that $\varphi^{\prime\prime}$ is given by (see (15) and (16)):

[TABLE]

3.1.3 Approximation of the GME function.

The following is devoted to obtain explicit approximations for the GME function $\varphi$ found in Section 3.1.1.

Let $\lambda>0$, $\beta>0$, $\gamma>0$ be given. Based on the assumption that problem (5) has a solution $\varphi$ that can be represented as:

[TABLE]

where $\varphi_{n}$ are real functions that must be determined, we will propose approximations $\varphi^{(N)}$ of the GME function given as:

[TABLE]

with $N\in\mathbb{N}_{0}$.

If $\varphi$ is given by (18), equation (5a) is formally equivalent to:

[TABLE]

where:

[TABLE]

Similarly, we have that (5b) is formally equivalent to:

[TABLE]

where:

[TABLE]

Therefore, if the functions $\varphi_{n}$ are such that:

[TABLE]

then the function $\varphi$ given by (18) is a formal solution to (5). Thus, functions $\varphi_{n}$ might be determined through problem (23). Let us observe it states that $\varphi_{0}$ must be a solution to:

[TABLE]

while each $\varphi_{n}$, $n\in\mathbb{N}$, must satisfy:

[TABLE]

with:

[TABLE]

Remark 3**.**

Observe that problem (24) coincides with (5) when $\beta=0$.

In the following we will only work with the zero and first order approximations $\varphi^{(0)}$, $\varphi^{(1)}$. From elementary results in ordinary differential equations (see, for example [24]), it can be obtained that the solution $\varphi_{0}$ to problem (24) is the function given by:

[TABLE]

Having obtained $\varphi_{0}$, we can now compute $\varphi_{1}$ through the problem (25) (see (26)), and obtain [24]:

[TABLE]

with:

[TABLE]

Therefore, we have that $\varphi^{(0)}=\varphi_{0}$ and $\varphi^{(1)}=\varphi_{0}+\beta\varphi_{1}$, with $\varphi_{0}$ and $\varphi_{1}$ given by (27), (28) respectively. In order to analyse the relation between each of these approximations and the GME function $\varphi$, we define the error $\mathcal{E}^{(N)}$ as:

[TABLE]

In the following, we will show that the GME function $\varphi$ converges uniformly to the zero order approximation $\varphi^{(0)}$, when $\beta\to 0$. As we shall see below, this is closely related to how the problem (5) depends on the parameter $\beta$.

Definition 3.1**.**

Let $0<b<\beta_{1}$. We will say that problem (5) is Lipschitz continuous on the parameter $\beta$ over the interval $[0,b]$ if there exists $L>0$ such that for any $b_{1},b_{2}\in[0,b]$ the following inequality holds:

[TABLE]

where $\varphi_{b_{1}}$, $\varphi_{b_{2}}$ are the only solutions in $K$ to problem (5) with parameters $b_{1}$, $b_{2}$, respectively.

It follows from Definition 3.1 that if problem (5) is Lipschitz continuous on $\beta$ over some interval $[0,b]$ with $0<b<\beta_{1}$, then the GME function $\varphi$ converges uniformly on $0<\eta<\lambda$ to the function $\varphi_{0}$ given by (27), when $\beta\to 0$ (see Remark 3). In other words, that $\mathcal{E}^{(0)}\to 0$ when $\beta\to 0$ and therefore, that $\varphi_{0}$ is a good approximation of the GME function $\varphi$ when the positive parameter $\beta$ is small enough. In the following we will prove that problem (5) is in fact Lipschitz continuous on the parameter $\beta$ over $[0,b]$ for any choice of $0<b<\beta_{1}$.

Lemma 3.4**.**

Let $h\in K$, $\eta\in[0,\lambda]$ and $b_{1},b_{2}\in[0,\beta_{1})$. We have:

[TABLE]

Proof.

Let $f$ be the real function defined on $\mathbb{R}^{+}_{0}$ by $f(x)=\exp(-2x)$ and:

[TABLE]

It follows from the Mean Value Theorem applied to function $f$ that:

[TABLE]

where $u$ is a real number between $x_{1}$ and $x_{2}$. Without any lost of generality, we will assume that $b_{2}\leq b_{1}$. Then, $x_{1}\leq x_{2}$. We have:

[TABLE]

Therefore,

[TABLE]

Then, we have:

[TABLE]

The final bound is now obtained by integrating the last expression and by the definition of $\beta_{1}$. ∎

Theorem 3.5**.**

Let $0<b<\beta_{1}$. The problem (5) is Lipschitz continuous on the parameter $\beta$ over the interval $[0,b]$.

Proof.

Let $\eta\in[0,\lambda]$, $b_{1},b_{2}\in[0,b]$ and $\varphi_{b_{1}},\varphi_{b_{2}}\in K$ the solutions to problem (5) with parameters $b_{1}$, $b_{2}$, respectively.

Taking into consideration that $\varphi_{b_{1}}$ and $\varphi_{b_{2}}$ are the fixed points of the operator $\tau_{b_{1}}$ and $\tau_{b_{2}}$ defined by (11) for $\beta=b_{1}$ and $\beta=b_{2}$, respectively, we have that:

[TABLE]

Now, from Lemmas 3.2, 3.4 it follows that:

[TABLE]

Since $g$ is an increasing function, $g(\beta_{1})=1$ and $b_{1}\leq b<\beta_{1}$, we have that $g(b_{1})\leq g(b)<1$. Then:

[TABLE]

∎

We end this Section by presenting some comparisons between the GME function $\varphi$, and its zero and first order approximations $\varphi^{(0)}$, $\varphi^{(1)}$. Figures 1(b) to 4(b) show the evolution of the error $\mathcal{E}^{(N)}$, and the plots of the GME function $\varphi$ against the best approximation obtained between $\varphi^{(0)}$ and $\varphi^{(1)}$. Each Figure correspond to one value of $\gamma=0.1,1,10,100$. The GME function $\varphi$ was obtained after solving problem (5) through the bvodes routine implemented in Scilab. Numerical computations were made by considering $\lambda\in[0,10]$ and a uniform mesh of step $0.01$ for the interval $[0,\lambda]$. For each choice of the parameter $\gamma$, equation (13) was numerically solved. The approximative solutions $\beta_{1}^{*}$ for each value of $\gamma$ are presented in Table 1. According to the election of $\gamma$, parameter $\ beta_{1}$ was set as $\beta_{1}^{*}$ as Table 1 states. From Figures 1(b) to 4(b) it can be seen that good agreement between the GME function $\varphi$ and either the zero or first order approximations $\varphi^{(0)}$, $\varphi^{(1)}$ can be obtained. They also suggest that the election between $\varphi^{(0)}$ and $\varphi^{(1)}$ it is mediated by the value of $\gamma$.

3.2 Analysis of equation (6).

We will now investigate the relation between the solution $\varphi$ to problem (5) found in Section 3.1.1 and the upper bound of its domain, that is the parameter $\lambda>0$, in order to analyse the existence of solution to equation (6). Throughout this Section, $\beta_{1}$ will refer to the only positive solution to equation (13).

Lemma 3.5**.**

Let $\gamma>0$, $\lambda>0$, $0\leq\beta<\beta_{1}$. If $\varphi$ is the GME function which belongs to $K$, then we have:

[TABLE]

where $g_{1}(\lambda)$, $g_{2}(\lambda)$ are given by:

[TABLE]

Proof.

2. i)

From the definition of $\varphi$ as the unique fixed point of the operator $\tau$ given by (11), we have that:

[TABLE]

Now the proof follows from:

[TABLE]

the bounds for $\Psi_{\varphi}$ given in Remark 2 and elementary boundedness techniques. 3. ii)

It is a direct consequence of the second inequality in i).

∎

Remark 4**.**

The first part of Lemma 3.5 together with the Squeeze Theorem implies that $\displaystyle\lim_{\lambda\to 0^{+}}\varphi^{\prime}(\lambda)=\frac{\gamma}{1+\beta}$.

Lemma 3.6**.**

Let $\gamma>0$, $0\leq\beta<\beta_{1}$. If $\varphi$ is the GME error function which belongs to $K$, then $\varphi^{\prime}$ is continuous on the parameter $\lambda>0$.

Proof.

After the change of variables:

[TABLE]

we have that problem (5) is equivalent to:

[TABLE]

It follows from the fixed point expression of $\varphi$ that $\varphi^{\prime}$ belongs to $K$. Applying Theorem 7.5 of [12] to the previous differential problem in the space $(0,1)\times(K\times K)\times\mathbb{R}^{+}$, we have that its solution is $C^{1}$ on the parameter $\lambda>0$. ∎

Theorem 3.6**.**

Let $\gamma>0$. If $0\leq\beta<\beta_{1}$ then equation (6) has at least one solution.

Proof.

For any $\lambda>0$, let $\varphi$ be the only solution in $K$ to problem (5) on the domain $[0,\lambda]$. Let also be $H$ the real function defined by:

[TABLE]

Since $H$ is a continuous function (see Lemma 3.6) that satisfies (see Lemma 3.5 and Remark 4):

[TABLE]

the theorem follows by recalling that the RHS of equation (6) is a positive number. ∎

In Figure 5(d) we present some plots for the function $H$ defined by (32). To compute $H$, is was considered the same numerical framework described at the end of Section 3.1.3 with $\lambda\in[0,5]$. Plots in Figure 5(d) suggest that the solution found in Theorem 3.6 is the unique solution to equation (6).

4 Relation with other Stefan problems

4.1 Relation with the Stefan problem with constant thermal conductivity

Let $\lambda>0$, $\gamma=2\text{Bi}>0$ be given. As it was already noted in Section 3.1.3, the solution to problem (5) when $\beta=0$ (constant thermal conductivity) is the function $\varphi_{0}$ given by (27). From this, we have that condition (6) can be written as:

[TABLE]

Therefore, the Stefan problem (1) with $\beta=0$ has the similarity solution $T$, $s$ given by:

[TABLE]

and (4b) if and only if $\lambda$ satisfies (34). This result has been already obtained in [33], where it was studied the phase-change process considered here but with constant thermal conductivity.

4.2 Relation with the Stefan problem with Dirichlet condition

Let us consider now the Stefan problem (1) with the Dirichlet boundary condition (1e*†*) in place of the convective one given by (1e). We will refer to it as problem (1†).

By following the same steps that led us to Theorem 2.1, we obtain that problem (1†) has the similarity solution $T$, $s$ given by (4) if and only if the function $\varphi$ and the parameter $\lambda$ satisfy condition (6) and the differential problem given by equation (5a), condition (5c) and:

[TABLE]

We will refer to the function $\varphi$ and the parameter $\lambda$ associated to problem (1†) as $\varphi^{\dagger}$ and $\lambda^{\dagger}$, respectively. Problem (1†) was studied in [11], where it was obtained almost the same similarity solution than the one presented above. In [11], the function $\varphi^{\dagger}$ is defined over $\mathbb{R}_{0}^{+}$ through equation (5a), condition (5b*†*) and $y(+\infty)=1$ (that is, as the ME function $\Phi$). This last change adds some extra conditions that must be satisfied by the temperature distribution $T$. But it is avoidable, as we are showing here. Let (5†) be the problem given by (5a), (5b*†*), (5c). In [8] it was proved the existence and uniqueness of the ME function $\Phi$ for small positive values of $\beta$ through a fixed point strategy. By performing the same analysis for problem (5†), we obtain that it has a unique non-negative bounded analytic solution $\varphi^{\dagger}$ for any given $\lambda^{\dagger}>0$. Moreover, we have that $\varphi^{\dagger}$ is the unique fixed point of the operator $\tau^{\dagger}$ from $K$ to itself defined by:

[TABLE]

where $C_{h}$ is given by:

[TABLE]

From the definitons of $\tau$ and $\tau^{\dagger}$ given in (11) and (36), it follows that $\tau h\to\tau^{\dagger}h$ (pointwise) when $\gamma\to+\infty$ for any function $h\in K$. Then, $\varphi\to\varphi^{\dagger}$ (pointwise) when $\gamma\to+\infty$. When we consider $\gamma=2Bi$ (see Theorem 3.3), $\gamma\to+\infty$ is equivalent to $h_{0}\to+\infty$. Thus, the solution to problem (1†) can be obtained as the limit case of the solution to problem (1) when the coefficient $h_{0}$ that characterizes the heat transfer coefficient at $x=0$ goes to infinity. This agrees well with the physical interpretation of temperature and convective boundary conditions (see Remark 1, [2, 7]).

We end this Section with some plots for the GME function $\varphi$. From Figure 6(d) it can be seen that it converges pointwise to the solution $\varphi^{\dagger}$ to problem (5†). By an abuse of notation, we have also referred to $\varphi^{\dagger}$ as ME function. The plots for both GME and ME functions were obtained after solving problems (5) and (5†) for $\lambda=10$ in the same numerical framework described at the end of Section 3.1.3. Although it was considered $\lambda=10$, functions were drawn over the interval $[0,3]$ aiming at a better visualization.

5 Conclusions

In this article we have presented an exact solution of similarity type for a one-phase Stefan problem with temperature-dependent thermal conductivity and a Robin boundary condition. The temperature distribution was defined through a Generalized Modified Error (GME) function. This was defined as the solution to a nonlinear boundary problem of second order, for which it was proved a result on existence and uniqueness of solutions. From this, the existence of similarity solutions was proved. It was also shown that results for the Stefan problems with either constant thermal conductivity or Dirichlet boundary conditions can be obtained as particular o limit cases of the results presented in this article.

Since the GME function is only available from numerical computations, it was proposed an strategy to obtain explicit approximations for it. Several values from the parameters involved in the physical problem were considered in the analysis of errors between the GME function and the two proposed approximations. The analysis performed suggest that the choice of the best approximation between those presented here depends on the values of the parameters. Nevertheless, good agreement can be obtained with both of them. From these explicit approximations, those for the temperature distribution can be obtained since it linearly depends on the GME function. In order to give some properties of the temperature distribution, it were also investigated some properties of the GME function. It was proved that it is a non-negative bounded analytic function which is increasing and concave, just as the classical error function is.

Acknowledgments

This paper has been partially sponsored by the Project PIP No. 0275 from CONICET-UA (Rosario, Argentina) and AFOSR-SOARD Grant FA 9550-14-1-0122.