3.1. Derivation of the leading-order system

Since we have taken a tall, thin limit, we expect that the lower interface is flat to leading order. To demonstrate this, we note that the leading-order terms in (3.3i) and (3.3j) are

(3.4a,b)\begin{equation} -p \left( \frac{\partial a}{\partial r}\right)^{2}=0, \quad \eta \frac{\partial u_z}{\partial r}\left( \frac{\partial a}{\partial r}\right)^{2}=0. \end{equation}

If we did not take ${\partial a}/{\partial r}=0$, then in addition to zero pressure $p=0$, we would also need the vertical component of velocity to be radially independent at the interface. Due to the no-flux condition on the wall, this would require zero vertical velocity across the whole interface, and hence throughout the entire fluid region. Thus, since we want to include convection in our model, we must have that the interface is flat to leading order.

Further, letting $a=a_0(t) + \delta a_1(r,t) +{O}(\delta ^{2})$ we note that the next-order term in (3.3i) is given by

(3.5)\begin{equation} -p_0 \left(\frac{\partial a_1}{\partial r} \right)^{2}-p_0=0, \end{equation}

for leading-order pressure $p_0$. We thus obtain a zero pressure condition on the flat, lower interface

(3.6)\begin{equation} p_0=0, \quad \mbox{at } z=-a_0(t). \end{equation}

By similar reasoning with (3.3k) and (3.3l), we also obtain a flat upper interface to leading order, so $b=b_0(t)+{O}(\delta )$, and we find the pressure condition

(3.7)\begin{equation} p_0=8b_0(t), \quad \mbox{at } z=-b_0(t). \end{equation}

We now expand the temperature, velocity, pressure and viscosity variables. Due to the $\delta ^{-2}$ term in the heat equation (3.3d), we consider expansions of the form

(3.8a)\begin{gather} T= T_0 +\delta^{2} T_1 +{O}(\delta^{4}), \end{gather}

(3.8b)\begin{gather}\boldsymbol{u}=\boldsymbol{u}_0+\delta^{2} \boldsymbol{u}_1 +{O}(\delta^{4}), \end{gather}

(3.8c)\begin{gather}p=p_{0} +\delta^{2} p_{1} +{O}(\delta^{4}), \end{gather}

(3.8d)\begin{gather}\eta=\eta_0 + \delta^{2} \eta _1 + {O}(\delta^{4}). \end{gather}

The leading-order thermal equation is

(3.9)\begin{equation} 0=\frac{1}{St} \frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial T_0}{\partial r} \right), \end{equation}

which, combined with either the boundedness of $T_0$ as $r \to 0$ or the vanishing normal derivative on the wall $r=1$, gives $T_0=T_0(z,t)$ (for $St={O}(1)$ in $\delta$, which we have assumed). This, in turn, gives that the leading-order viscosity is also radially independent, with

(3.10)\begin{equation} \eta_0(z,t)=\exp(-\lambda (T_0(z,t)-1) ). \end{equation}

The leading-order fluid equations are

(3.11a)\begin{gather} \frac{\partial u_{z0}}{\partial z}+\frac{1}{r}\frac{\partial}{\partial r}\left(r u_{r0} \right)=0, \end{gather}

(3.11b)\begin{gather}0=-\frac{\partial p_{0}}{\partial z}+\frac{1}{r}\frac{\partial}{\partial r} \left(r \eta_0 \frac{\partial u_{z0}}{\partial r} \right)-8, \end{gather}

(3.11c)\begin{gather}0=-\frac{\partial p_{0}}{\partial r}, \end{gather}

(3.11d)\begin{gather}u_{r0}(1,z,t)=u_{z0}(1,z,t)=0, \end{gather}

(3.11e)\begin{gather}u_{r0}, \quad u_{z0} \mbox{ bounded as } r\to 0, \end{gather}

(3.11f)\begin{gather}p_0(r,-a_0(t),t)=0, \quad p(r,-b_0(t),t)=8b_0(t). \end{gather}

Since the leading-order pressure $p_0$ is also radially independent, (3.11b) can be integrated twice to give

(3.12)\begin{equation} u_{z0}=-\frac{1}{\eta_0(z,t)}\left(\frac{\partial p_0}{\partial z}+8 \right)\left(\frac{1-r^{2}}{4} \right). \end{equation}

Integrating the incompressibility condition (3.11a) radially gives

(3.13)\begin{equation} \int_{0}^{1} r\frac{\partial u_{z0}}{\partial z}\, \textrm{d} r+[r u_{r0} ]_{r=0}^{r=1}=0. \end{equation}

The boundary conditions on $u_{r0}$, combined with (3.12), give

(3.14)\begin{equation} \frac{\partial}{\partial z}\left( \frac{1}{\eta_0}\left( \frac{\partial p_0}{\partial z}+8\right) \right)=0,\end{equation}

and hence $u_{z0}$ is independent of height $z$. Substituting $u_{z0}$ into the incompressibility condition (3.11a), we find that $u_{r0}=0$. Thus, to leading order, the fluid moves downward with a velocity that is the same at each height, yet which varies radially and temporally. We have that

(3.15)\begin{equation} \boldsymbol{u}_0= -2 f(t)(1-r^{2})\boldsymbol{k}, \end{equation}

where $f(t)$ can be interpreted as the average speed of the fluid, since

(3.16)\begin{equation} \left|\frac{2\pi \displaystyle\int_{r=0}^{1} 2 f(t) (r^{2}-1) r \, \textrm{d} r}{2 \pi \displaystyle\int_{r=0}^{1} r \, \textrm{d} r}\right|=f(t). \end{equation}

We can find $f(t)$ by substituting (3.15) into (3.11b), which when rearranged gives

(3.17)\begin{equation} 8 f(t) \eta_0=8+\frac{\partial p_0}{\partial z}. \end{equation}

Integrating (3.17) vertically between $z=-a_0(t)$ and $z=-b_0(t)$, we have

(3.18)\begin{align} 8 f(t) \int_{-a_0(t)}^{-b_0(t)}\eta_0 (z,t) \, \textrm{d} z =[p_0 + 8 z ]_{z=-a_0(t)}^{z=-b_0(t)}=8b_0(t)+8(-b_0(t) + a_0(t)), \end{align}

and hence we obtain

(3.19)\begin{equation} f(t)=\frac{a_0(t)}{\displaystyle\int_{z=-a_0(t)}^{-b_0(t)}\eta_0(z,t) \, \textrm{d} z}. \end{equation}

Thus, $f(t)$ is a ratio of the gravitational force in the numerator driving the fluid downwards to the viscous force in the denominator resisting motion. Note that the weight of the whole charge is felt, but the resistance is due only to the viscosity of the lower region. We note that in the case of constant viscosity ($\lambda =0$), and no granular region ($T_b=1$), we have that $f(t)\equiv 1$, and we recover the usual Poiseuille flow results.

We determine an equation for $T_0(z,t)$ by considering the ${O}(1)$ heat equation, namely

(3.20)\begin{equation} \frac{\partial T_0}{\partial t}+\frac{Pe}{St}u_{z0} \frac{\partial T_0}{\partial z}=\frac{1}{St} \left(\frac{1}{r}\frac{\partial}{\partial r}\left(r \frac{\partial T_1}{\partial r} \right) +\frac{\partial^{2} T_0}{\partial z^{2}}\right). \end{equation}

Integrating (3.20) radially, noting that $T_0$ is independent of $r$, gives

(3.21)\begin{equation} \frac{1}{2} \frac{\partial T_0}{\partial t}-\frac{1}{2}\frac{Pe}{St} f(t) =\frac{1}{St} \int_{r=0}^{1} \frac{\partial}{\partial r}\left(r \frac{\partial T_1}{\partial r} \right)\, \textrm{d} r +\frac{1}{2}\frac{1}{St}\frac{\partial^{2} T_0}{\partial z^{2}}. \end{equation}

Since $T_1$ is bounded at $r \to 0$ and has a zero normal derivative on the wall $r=1$, we are left with

(3.22)\begin{equation} \frac{\partial T_0}{\partial t}-\frac{Pe}{St}f(t)\frac{\partial T_0}{\partial z}=\frac{1}{St}\frac{\partial^{2} T_0 }{\partial z^{2}}, \quad \mbox{in } -a_0(t) \leq z \leq 0, \end{equation}

with the corresponding initial and boundary conditions

(3.23)\begin{equation} T_0(z,0)=0,\quad T_0(-a_0(t),t)=1, \quad T_0(0,t)=0. \end{equation}

We note that the location of the interface $z=-b(t)$ is determined by

(3.24)\begin{equation} T_0(-b_0(t),t)=T_b. \end{equation}

We also need to determine the leading-order Stefan condition. Considering the leading-order terms in (3.3m) gives

(3.25)\begin{equation} -\frac{\textrm{d} a_0}{\textrm{d} t}+\frac{Pe}{St} 2f(t)(1-r^{2})=\gamma +\frac{\partial T_0}{\partial z}, \quad \mbox{at } z=-a_0(t), \end{equation}

which cannot be satisfied by the leading-order solution, because $a_0$ and $T_0$ are independent of $r$ and $f(t) >0$, so nothing can balance the single term with radial dependence. The problem arises because we have shown that the lowest-order interface $a_0(t)$ is flat, and that the lowest-order temperature $T_0$ is radially independent. However, at higher orders, the interface is not flat, and varies radially. We thus we must consider an inner region adjacent to the interface where temperature could potentially vary radially. We expect this region to be of thickness ${O}(\delta )$, and hence in a domain that is not tall and thin. We now consider this boundary region to find that we should impose a radially averaged boundary condition on the lowest-order outer problem.

3.2. Averaged velocity from inner region

We consider the small inner region at the bottom of the lower, viscous region by defining the inner variables $a_I(r,t), \hat {z}$, and $T_I(\hat {z},r,t)$ as

(3.26a–c)\begin{equation} a=a_0(t) + \delta a_I(r,t), \quad T=1+\delta T_I(\hat{z},r,t), \quad z=-a_0(t) +\delta \hat{z}, \end{equation}

where $r, t \geq 0$ and $\hat {z} \geq -a_I$. We have assumed that the correction term lies below the leading-order interface, i.e. $a_I >0$. A schematic diagram of the inner region is given in figure 3.

Figure 3. Schematic of the inner region near $z=-a_0(t)$.

The leading-order temperature equation in the inner region is given by

(3.27)\begin{equation} \hat{\nabla}^{2} T_I=0, \quad \hat{z} \geq -a_I(r,t), \end{equation}

where the hat on the Laplacian denotes derivatives with respect to $r$ and $\hat {z}$. We have the inner velocity profile, $\boldsymbol {u}_I$, which scales with the outer velocity and satisfies the incompressibility condition $\hat {\boldsymbol {\nabla }} \boldsymbol{\cdot} \boldsymbol {u}_I=0$.

The leading-order Stefan condition in the inner region is given by

(3.28)\begin{equation} -\frac{\textrm{d} a_0(t)}{\textrm{d} t}+\frac{Pe}{St} \boldsymbol{u}_I \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} (-a_I) =\gamma-\hat{\boldsymbol{\nabla}} T_I \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}}(- a_I), \quad \mbox{at } \hat{z}=-a_I, \end{equation}

or, alternatively, by

(3.29)\begin{equation} -\frac{1}{\| \hat{\boldsymbol{\nabla}} (-a_I) \|}\frac{\textrm{d} a_0(t)}{\textrm{d} t}+\frac{Pe}{St} \boldsymbol{u}_I \boldsymbol{\cdot} \boldsymbol{n}_I =\frac{\gamma}{\| \hat{\boldsymbol{\nabla}} (-a_I) \|}-\hat{\boldsymbol{\nabla}} T_I \boldsymbol{\cdot} \boldsymbol{n}_I, \quad \mbox{at } \hat{z}=-a_I, \end{equation}

where we have defined the downwards pointing inner normal $\boldsymbol {n}_I=\hat {\boldsymbol {\nabla }} (-a_I)/\|\hat {\boldsymbol {\nabla }}(-a_I)\|$. We also have that

(3.30)\begin{equation} T_I=0, \quad \mbox{at } \hat{z}=a_I, \quad \frac{\partial T_I}{\partial r}=0=u_{Iz}=u_{Ir}, \quad \mbox{ at } r=1, \end{equation}

and the matching conditions with the leading-order outer variables $T_0$ and $\boldsymbol {u}_0$. Using Van Dyke's matching principle, noting that $T_{0z}(-a_0(t) + \delta \hat {z},t) = T_{0z}(-a_0(t),t)+{O}(\delta )$, we have

(3.31)\begin{equation} \frac{\partial T_I}{\partial \hat{z}}(\infty,r,t) =\frac{\partial T_0}{\partial z}(-a_0(t),t), \quad \boldsymbol{u}_I(\infty,r,t)=\boldsymbol{u_0}(-a_0(t),t). \end{equation}

Integrating (3.29) along the surface $\hat {z}=-a_I(r,t)$, we obtain

(3.32)\begin{align} \iint_{\hat{z}=-a_I}\left( - \frac{1}{\| \hat{\boldsymbol{\nabla}} (-a_I) \|}\frac{\textrm{d} a_0(t)}{\textrm{d} t}+\frac{Pe}{St} \boldsymbol{u}_I \boldsymbol{\cdot} \boldsymbol{n}_I \right)\, \textrm{d} S =\iint_{\hat{z}=-a_I} \left(\frac{\gamma}{\| \hat{\boldsymbol{\nabla}} (-a_I) \|}-\hat{\boldsymbol{\nabla}} T_I \boldsymbol{\cdot} \boldsymbol{n}_I \right) \, \textrm{d} S. \end{align}

Using the divergence theorem in the region bounded by $\hat {z}=a_I, \hat {z}=\infty$ and the cylindrical wall $r=1$, and noting that $\hat {\nabla }^{2} T_I=0$, we find

(3.33)\begin{equation} \iint_{\hat{z}=-a_I}\hat{\boldsymbol{\nabla}} T_I \boldsymbol{\cdot} \boldsymbol{n}_I \, \textrm{d} S = -\iint_{\hat{z}=\infty}\hat{\boldsymbol{\nabla}}T_I \boldsymbol{\cdot} \boldsymbol{k} \, \textrm{d} S. \end{equation}

Since $\hat {\boldsymbol {\nabla }}T_I \to \boldsymbol {\nabla } T_0 (-a_0(t),t)$ as $\hat {z}\to \infty$, we obtain

(3.34)\begin{equation} \iint_{\hat{z}=-a_I}\hat{\boldsymbol{\nabla}}T_I \boldsymbol{\cdot} \boldsymbol{n}_I \, \textrm{d} S =-\int_{\theta=0}^{2\pi} \int_{r=0}^{1} r\frac{\partial T_0}{\partial z}(-a_0(t),t) \, \textrm{d} r \, \textrm{d} \theta= -\pi \frac{\partial T_0}{\partial z}(-a_0(t),t). \end{equation}

Similarly, since $\boldsymbol {u}_I$ is incompressible and there are zero flux conditions on the walls, we can use the divergence theorem to write

(3.35)\begin{equation} \iint_{\hat{z}=-a_I}\boldsymbol{u}_I \boldsymbol{\cdot} \boldsymbol{n}_I \, \textrm{d} S =-\iint_{\hat{z}=\infty}\boldsymbol{u}_I \boldsymbol{\cdot} \boldsymbol{k} \, \textrm{d} S = -\int_{\theta=0}^{2\pi} \int_{r=0}^{1} r u_{z0}(-a_0(t),t) \, \textrm{d} r \, \textrm{d} \theta= \pi f(t), \end{equation}

where we have used $\boldsymbol {u}_0=-2f(t)(1-r^{2})\boldsymbol {k}$ from the outer problem (3.15).

Noting that

(3.36)\begin{equation} \iint_{\hat{z}=-a_I} \frac{1}{\| \hat{\boldsymbol{\nabla}}(-a_I) \|} \, \textrm{d} S = \iint_{\hat{z}=-a_I} -\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{n}_I \, \textrm{d} S = \pi, \end{equation}

is the projected area of the surface onto the plane orthogonal to $\boldsymbol {k}$, and that ${\textrm {d} a_0}/{\textrm {d} t}$ and $\gamma$ are constant with respect to the integrals, we find

(3.37)\begin{equation} -\frac{\textrm{d} a_0(t)}{\textrm{d} t}+\frac{Pe}{St}f(t)=\gamma + \frac{\partial T_0}{\partial z}, \quad \mbox{at } z=-a_0(t). \end{equation}

Thus our analysis of the inner problem has identified that we should impose a radially averaged boundary condition (3.37), instead of (3.25), on the flat outer position of the interface for the lowest-order outer variables. Equation (3.37) is an averaged Stefan condition and closes the leading-order outer problem.

3.3. Statement of the model for a tall, thin cylinder

We simplify the notation used in the previous sections when we derived the lowest-order behaviour by dropping the $0$ subscripts, so that the problem we seek to study is given by

(3.38a)\begin{gather} \frac{\partial T}{\partial t}-\frac{f(t) Pe}{ St}\frac{\partial T}{\partial z}=\frac{1}{St}\frac{\partial^{2} T}{\partial z^{2}},\quad \mbox{for } -a(t) \leq z \leq 0, \end{gather}

(3.38b)\begin{gather}-\frac{\textrm{d} a}{\textrm{d} t}+\frac{f(t) Pe}{St}=\gamma + \frac{\partial T}{\partial z}, \quad \mbox{at } z=-a(t), \end{gather}

(3.38c)\begin{gather}T(0,t)=0, \end{gather}

(3.38d)\begin{gather}T(-b(t),t)=T_b, \end{gather}

(3.38e)\begin{gather}T(-a(t),t)=1, \end{gather}

(3.38f)\begin{gather}T(z,0)=0, \end{gather}

(3.38g)\begin{gather}a(0)=1, \quad b(0)=1, \end{gather}

where

(3.38h)\begin{equation} f(t)=\frac{a(t)}{\displaystyle\int_{z=-a(t)}^{-b(t)}\exp(-\lambda(T(z,t)-1)) \, \textrm{d} z}. \end{equation}

Upon solving this system for $T(z,t), a(t), b(t)$ and $f(t)$, we may determine the fluid velocity, which is given by

(3.39)\begin{equation} \boldsymbol{u}=-2f(t)(1-r^{2})\boldsymbol{k}, \end{equation}

as well as the volumetric flow rate of material across the interface $z=-a(t)$, which is given by

(3.40)\begin{align} Q(t) &= -\int_{\theta=0}^{2\pi} \int_{r=0}^{1}\left(\frac{\textrm{d} a}{\textrm{d} t} -\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{n}\right) r \, \textrm{d} r \, \textrm{d} \theta \nonumber\\ &= -2\pi \int_{\theta=0}^{\pi} \int_{r=0}^{1} \left( \frac{\textrm{d} a}{\textrm{d} t}-2f(t)(1-r^{2})\right)r \,\mathrm{d}r \nonumber\\ &= \pi\left( f(t)-\frac{\textrm{d} a}{\textrm{d} t}\right). \end{align}

Here we have chosen to define $Q>0$ in (3.40) when material is dripping from the base of the fluid. A schematic diagram of the tall, thin model is given in figure 4.

Figure 4. Schematic of the problem geometry for the tall, thin cylinder model (3.38). The viscous region (blue) is separated from the granular region (red) above by the isotherm $T=T_b$ at $z=-b(t)$ and is bounded below by the isotherm $T=1$ at $z=-a(t)$. The top of the granular region is cold, $T=0$, and at a fixed height $z=0$. The moving boundaries have initial conditions $a(0)=b(0)=1$.

4.1. Steady-state solutions

For our model, we find that steady-state solutions cannot be written in an explicit closed form, but depend upon an implicit equation for $f$ that arises from (3.38h), which must be solved numerically. We find two roots to this equation, corresponding to two steady states. We show numerical evidence suggesting that one of these steady states is stable, while the other steady state is unstable. Subsequently, in § 5.2, we plot bifurcation diagrams of the model behaviour in parameter space.

Suppose we look for a steady-state solution, so that $a(t)=a^{*}, b(t)=b^{*}, T(z,t)=T^{*}(z)$ and $f(t)=f^{*}$. The system (3.38) then becomes

(4.1a)\begin{gather} -f^{*} Pe \frac{\textrm{d} T^{*}}{\textrm{d} z} =\frac{\textrm{d}^{2} T^{*}}{\textrm{d} z^{2}}, \quad \mbox{in } -a^{*} \leq z \leq 0, \end{gather}

(4.1b)\begin{gather}\frac{f^{*} Pe}{St}=\gamma+\frac{\partial T^{*}}{\partial z}, \quad \mbox{at } z=-a^{*}, \end{gather}

(4.1c)\begin{gather}T^{*}(-a^{*})=1, \quad T^{*}(0)=0, \end{gather}

(4.1d)\begin{gather}T^{*}(-b^{*})=T_b, \end{gather}

(4.1e)\begin{gather}f^{*}=\frac{a^{*}}{\displaystyle\int_{z=-a^{*}}^{-b^{*}}\exp\left(-\lambda(T^{*}-1) \right)\, \textrm{d} z}. \end{gather}

The equation (4.1a) combined with the boundary conditions (4.1c) can be solved to find the steady-state temperature profile

(4.2)\begin{equation} T^{*}(z)=\frac{\exp(-f^{*}Pe\, z)-1}{\exp(\,f^{*}Pe\, a^{*})-1}. \end{equation}

We can then use the Stefan condition (4.1b) to find an equation for the position of the free boundary, which reads

(4.3)\begin{equation} \frac{f^{*} Pe}{ St}=\gamma-f^{*} Pe\frac{\exp(\,f^{*} Pe\, a^{*})}{\exp(\,f^{*} Pe\, a^{*})-1}. \end{equation}

This can be rearranged to determine that the lower interface is at

(4.4)\begin{equation} a^{*}=\frac{1}{f^{*} Pe}\ln\left(\frac{\dfrac{\gamma }{f^{*} Pe}-\dfrac{1}{St}}{\dfrac{\gamma }{f^{*} Pe}-\dfrac{1}{St}-1}\right).\end{equation}

We use the isotherm (4.1d) to find that the upper interface is at

(4.5)\begin{equation} b^{*}=\frac{1}{f^{*} Pe}\ln(1+T_b(\exp(\,f^{*} Pe\, a^{*})-1 )). \end{equation}

When $T_b=0$, the upper interface is at the top of the domain ($b^{*}=0$). When $T_b=1$, we find that (4.5) simplifies to give $b^{*}=a^{*}$, which is also expected. Note that we have not yet found $a^{*}$ and $b^{*}$ explicitly as functions of the parameters $Pe, St, \gamma , T_b$, and $\lambda$, since (4.4) and (4.5) still have a dependence on $f^{*}$. The crucial equation we have not yet used is (4.1e), which we write as the implicit integral equation for $f^{*}$

(4.6)\begin{equation} f^{*}=\frac{a^{*}}{\displaystyle\int_{z=-a^{*}}^{-b^{*}}\exp\left(-\lambda\left(\frac{\exp\left(-f^{*} Pe\, z\right)-1}{\exp\left(\,f^{*} Pe\, a^{*}\right)-1}-1\right) \right)\, \textrm{d} z}.\end{equation}

We thus seek the solutions $a^{*}, b^{*}$, and $f^{*}$ to the coupled equations (4.4)–(4.6), noting that if we can find $f^{*}$ we can subsequently obtain $a^{*}$ and $b^{*}$ through substitution.

To help solve the problem (4.4)–(4.6) we make the change of variable

(4.7)\begin{equation} \theta =\frac{\gamma }{f^{*} Pe}-\frac{1}{St}. \end{equation}

This allows us to write (4.4), (4.5) and (4.6), respectively, as

(4.8)\begin{gather} a^{*}=\frac{1}{f^{*}Pe}\ln\left(\frac{\theta}{\theta-1}\right)=\frac{1}{\gamma}\left(\theta+\frac{1}{St} \right)\ln\left(\frac{\theta}{\theta-1}\right), \end{gather}

(4.9)\begin{gather}b^{*}=\frac{1}{f^{*}Pe}\ln \left(1+\frac{T_b}{\theta -1} \right)=\frac{1}{\gamma}\left(\theta+\frac{1}{St} \right)\ln \left(1+\frac{T_b}{\theta -1} \right), \end{gather}

(4.10)\begin{gather}f^{*} \textrm{e}^{\lambda} \int_{z=-a^{*}}^{-b^{*}}\exp \left(-\lambda (\theta-1) \left(\exp\left(-f^{*} Pe\, z\right)-1 \right) \right)\, \textrm{d} z=\frac{1}{\gamma}\left(\theta+\frac{1}{St} \right)\ln\left(\frac{\theta}{\theta-1}\right). \end{gather}

We now look for the solutions $a^{*}, b^{*}$ and $\theta$ of the system (4.8)–(4.10).

We must take care to ensure the logarithms in (4.8)–(4.10) are well-defined. The parameters $St, Pe$ and $\gamma$ are all positive, $T_b \in [0,1]$ and $f^{*}$ is non-negative. The logarithm in (4.8) and (4.10) is not well defined for $\theta \in [0,1]$. In the case that $\theta <0$, we have $a^{*}<0$ and thus a steady state at $z=-a^{*}$ which is above the top of the domain, which is evidently not physically attainable. The logarithm in (4.9) is not defined for $\theta \in [1-T_b,1]$. For $\theta <1-T_b$ we have $b^{*}<0$, which is also not physically attainable. Hence, we must have $\theta >1$ for physically relevant steady states to exist.

We seek to simplify (4.10), by making the substitution

(4.11)\begin{equation} w=\lambda(\theta-1)\exp(-f^{*} Pe\, z),\end{equation}

which transforms (4.10) to

(4.12)\begin{equation} -\frac{ \textrm{e}^{\lambda \theta}}{Pe} \int_{w=\lambda \theta}^{\lambda(\theta+T_b-1)} \frac{\textrm{e}^{-w} }{w}\, \textrm{d} w=\frac{1}{\gamma}\left(\theta+\frac{1}{St} \right)\ln\left(\frac{\theta}{\theta-1}\right).\end{equation}

Using the exponential integral $E_1(x)$, given by

(4.13)\begin{equation} E_1(x)=\int_{x}^{\infty} \frac{\textrm{e}^{-t}}{t} \, \textrm{d} t, \end{equation}

we can rewrite (4.12) as

(4.14)\begin{equation} \textrm{e}^{\lambda \theta}( E_1(\lambda(\theta-1+T_b))-E_1( \lambda \theta) )=\frac{Pe}{\gamma}\left(\theta+\frac{1}{St}\right)\ln\left(\frac{\theta}{\theta-1}\right). \end{equation}

Since there is a one-to-one mapping between $f^{*}$ and $\theta$ (defined in (4.7)), each root of (4.14) for $\theta$ corresponds to a steady state $f^{*}$. Note also that $Pe/\gamma$ acts as one parameter in (4.14), and so there are only four parameter groupings in (4.14), being $St, T_b, \lambda$ and $Pe/\gamma$. In exploring this parameter regime we recall we require $\theta >1$, which implies $f^{*} < {\gamma }/{Pe (1+1/St)}$.

Numerical evidence suggests the existence of at most one stable steady state for given parameters, although we shall now show that multiple steady states are possible, with the additional steady state seemingly unstable. We introduce the function

(4.15)\begin{equation} \phi(\theta)= \textrm{e}^{\lambda \theta}\left( E_1\left(\lambda(\theta-1+T_b)\right)-E_1\left( \lambda \theta\right) \right)-\frac{Pe}{\gamma}\left(\theta+\frac{1}{St}\right)\ln\left(\frac{\theta}{\theta-1}\right), \end{equation}

so that each zero of (4.15) corresponds to a root of (4.14). If the partial derivatives of $\phi$ were all single-signed then it could have at most one root. This is true for the derivatives with respect to $St, Pe/\gamma$, and $T_b$, but unfortunately not for $\lambda$ and $\theta$. Instead, to determine the number of roots of (4.14), we plot contours of roots of $\phi$ for $\lambda$ against $\theta$, for different values of the other parameters. If the contours are monotonic, then only one steady state is possible. However, if several possible $\theta$ exist for each $\lambda$, then multiple steady states are possible. In the parameter space explored, we find curves exhibiting the existence of two steady states, and give an example of this in figure 5, circling the location of two roots for $T_b=0.8$ and $\lambda =1$. We plot values of $\lambda$ within the range $\lambda \in [-10,10]$ for demonstrative purposes, although in practice we expect the viscosity of the fluid to increase as temperature decreases, and hence are only interested in $\lambda \geq 0$. We shall later show through numerical simulations that solutions tend toward only one of these steady state, with the other being unstable.

Figure 5. Contour plots of roots of (4.15) for $\theta$, given a value of $\lambda$ and $T_b$, for $Pe=1, St=2$ and $\gamma =30$. The multiple roots for $T_b=0.8$ and $\lambda =1$ are circled.

For each fixed value of $T_b$ in figure 5, there are values of $\lambda$ for which one root occurs for $\theta$ close to $1$, and the other for large $\theta$. We now examine each of these limits.

For $\theta$ near 1, we note that

(4.16)\begin{equation} \frac{\gamma}{f^{*} Pe}-\dfrac{1}{St} \sim 1. \end{equation}

In particular, for the large Stefan number limit where there is small latent heat, we see that this corresponds to heating and convection being in approximate balance. To analyse the case of $\theta \sim 1$ more formally, we set $\theta =1+\theta _{\epsilon }$, where $0< \theta _{\epsilon } \ll 1$. From (4.7) we have

(4.17)\begin{equation} f^{*}=\frac{\gamma}{Pe}\left(\frac{1}{1+\dfrac{1}{St}}\right)+{O}(\theta_{\epsilon}). \end{equation}

Substituting $\theta =1+\theta _{\epsilon }$ into (4.8) and (4.9), and Taylor expanding the relevant quantities, we find

(4.18)\begin{align} \left.\begin{array}{c@{}} \displaystyle a^{*}=\dfrac{1}{\gamma}\left( 1+ \theta_{\epsilon} + \dfrac{1}{St}\right) \ln \left(\dfrac{1+\theta_{\epsilon}}{1+\theta_{\epsilon}-1} \right)=-\dfrac{1}{\gamma}\left(1+\dfrac{1}{St}\right) \ln \theta_{\epsilon}+ {O}\left(\theta_{\epsilon}\ln \theta_{\epsilon} \right), \\ \displaystyle b^{*}=\dfrac{1}{\gamma}\left( 1+ \theta_{\epsilon} + \dfrac{1}{St}\right) \ln \left(1+\dfrac{T_b}{1+\theta_{\epsilon}-1} \right)=-\dfrac{1}{\gamma}\left(1+\dfrac{1}{St}\right) \ln \theta_{\epsilon}+ {O}\left(1 \right). \end{array}\right\}\end{align}

We can find the value of $\theta _{\epsilon }$ in terms of the dimensionless parameters by expanding (4.14) in $\theta _{\epsilon }$, namely

(4.19)\begin{equation} \textrm{e}^{\lambda}(E_1(\lambda T_b)-E_1(\lambda))=\frac{Pe}{\gamma}\left(1+\frac{1}{St}\right)\ln \left(\frac{1}{\theta_{\epsilon}} \right)+{O}(\theta_{\epsilon}), \end{equation}

to give the approximation

(4.20)\begin{equation} \theta_{\epsilon}\approx \exp \left[ -\left(\frac{\textrm{e}^{\lambda}(E_1(\lambda T_b)-E_1(\lambda))}{\dfrac{Pe}{\gamma}\left(1+\dfrac{1}{St}\right) }\right) \right], \end{equation}

and hence the approximation (4.18) is valid if the parameter values $\lambda , T_b, St, \gamma$ and $Pe$ ensure that $\theta _{\epsilon } \ll 1$.

We note from (4.18) that $a^{*}$ and $b^{*}$ can become arbitrarily large as $\theta _{\epsilon }\searrow 0$, corresponding to $\theta \searrow 1$. Hence, reducing the heating, or increasing the convective term, one can lower the steady states. If the heating is too low, however, so that $\theta <1$, then no steady state exists and the interfaces fall indefinitely. We later show that numerical simulations suggest that the root close to $\theta =1$ corresponds to the unstable steady state.

In figure 5 we observe there is also a root for large $\theta$, so

(4.21)\begin{equation} \frac{\gamma}{f^{*} Pe}-\frac{1}{St} \gg 1. \end{equation}

This corresponds physically to large heating or small convection. To find an approximation to this root we consider $1/\theta \to 0$, and use (4.7) to find

(4.22)\begin{equation} f^{*}=\frac{\gamma}{Pe}\frac{1}{\theta}+{O}\left(\frac{1}{\theta^{2}}\right). \end{equation}

We can then expand the logarithms in (4.8) and (4.9), to obtain

(4.23)\begin{equation} \left.\begin{array}{c@{}} a^{*} =\dfrac{1}{\gamma}+{O}\left(\dfrac{1}{\theta}\right),\\ b^{*}=\dfrac{T_b}{\gamma}+{O}\left(\dfrac{1}{\theta}\right). \end{array}\right\}\end{equation}

Continuing on, we expand (4.14) in ${1}/{\theta }$, namely

(4.24)\begin{equation} \frac{1}{\theta} \left(\frac{\textrm{e}^{\lambda(1-T_b)}-1}{\lambda}\right) =\frac{Pe}{\gamma}+\frac{1}{\theta}\frac{Pe}{\gamma}\left(\frac{1}{2}+ \frac{1}{St}\right)+{O}\left({\frac{1}{\theta^{2}}}\right), \end{equation}

to find

(4.25)\begin{equation} \theta\approx \frac{\gamma}{Pe}\left(\frac{\textrm{e}^{\lambda(1-T_b)}-1}{\lambda}\right)-\left(\frac{1}{2}+\frac{1}{St}\right). \end{equation}

From (4.23) we see that $a^{*} \to 1/\gamma , b^{*} \to T_b/ \gamma$ as $\theta \to \infty$, so we can determine how close to the top of the domain the interfaces will be if input parameters are chosen to have very large heating or small convection. Note also that (4.23) gives us a possible scaling for the dimensional thickness of $a^{*}$ and $b^{*}$ as

(4.26a)\begin{gather} \mbox{dimensional thickness of }a^{*}=\frac{k(T_d-T_c)}{\epsilon_R \sigma (T_h^{4}-T_d^{4})}, \end{gather}

(4.26b)\begin{gather}\mbox{dimensional thickness of }b^{*}=\frac{k(T_k-T_c)}{\epsilon_R \sigma (T_h^{4}-T_d^{4})}. \end{gather}

We observe that the initial dimensional height $h$ has scaled out of the problem, which is reassuring since this is artificially imposed. Although in practice there is a floor where the dripped material accumulates, there is nothing in this model to prevent the interfaces from moving below the dimensional initial condition of $z=-h$. We suspect that the root corresponding to large $\theta$ corresponds to the stable steady state.

4.2. Asymptotic limits of the tall, thin model

We now consider special limits of the tall, thin model by exploiting the possible sizes of the dimensionless parameters $T_b, Pe, St, \gamma$ and $\lambda$. In each case, we assume that the aspect ratio $\delta$, used to derive the tall, thin model in § 3, is smaller than any other parameter we subsequently take to be small. We list some possible dominant balances which admit steady states in table 1, indicating which parameters we take to be small in each case. These distinguished limits are as follows:

• • Case A: the material is falling slowly, since the Péclet number is small. This could be due to high viscosity or low density, or alternatively high thermal diffusivity.

• • Case B: the radiative heating from underneath the charge is very large.

• • Case C: the radiative heating is large and the Stefan number is small. This could correspond to a small temperature difference between the top and bottom of the charge.

• • Case D: the Stefan number is large. This could be caused by very small latent heat.

• • Case E: the radiative heating is large and the Péclet number is small, with the Stefan number being extremely small, balancing the ratio of convection to heating.

• • Case F: the convective term is small and does not balance with other parameters, so this gives a limit without convection.

Table 1. Table of dominant balances which give steady states in the model (3.38).

We examine each of these limits in turn. The ratio $(a-b)/a$ gives the relative thickness of the crust to the total size of the charge, and so we start by analysing balances A and B, as these lead to a thin crust. Then, we analyse balances C, D, E and F, as these cases give rise to an ${O}(1)$ crust thickness.

To examine cases A and B we introduce a small parameter, $\epsilon$, defined by

(4.27)\begin{equation} \epsilon = 1-T_b, \end{equation}

and consider the limit $0 < \epsilon \ll 1$. Since we have already taken the aspect ratio $\delta$ to be small, we note that we are taking the limit $0 < \delta \ll \epsilon \ll 1$, so that the thin crust region is still much taller than it is wide. We find it easier to introduce the small parameter in terms of temperature rather than height, but we assume that for ${O}(1)$ time, $b(t)-a(t)={O}(\epsilon )$, and so we introduce the ${O}(1)$ variable $c(t)>0$ such that

(4.28)\begin{equation} b(t)=a(t)-\epsilon c(t). \end{equation}

If we substitute $T_b=1-\epsilon$ into (4.14), and expand in terms of $\epsilon$, assuming all other parameters are ${O}(1)$, we obtain

(4.29)\begin{equation} \frac{\epsilon}{\theta}=\frac{Pe}{\gamma}\left(\theta+\frac{1}{St}\right) \ln\left(\frac{\theta}{\theta-1}\right)+{O}(\epsilon^{2}).\end{equation}

We need to find a balance for the $\epsilon$ term, which cannot be achieved if $Pe, St, \gamma , \theta$ and $\ln (\theta /(\theta -1))$ are all ${O}(1)$ parameters. We note that $Pe, \gamma$ and $St$ are all positive, and that in previous analysis we showed that $\theta >1$ is required for physically realistic steady states (i.e. $a^{*},b^{*} \geq 0$). There are several limits in which a necessary balance can be found, and we will examine each of these in turn.

After considering cases A and B, we will consider dominant balances for which the relative size of the crust to the total charge is not small, but order one. These are cases C, D, E and F.

4.2.1. Case A: small Péclet number, $Pe={O}(\epsilon )$

If the Péclet number is small, then thermal convection is small compared to thermal conduction, and so the material is falling relatively slowly. We scale $Pe=\epsilon \widehat {Pe}$, where $\widehat {Pe}$ is an ${O}(1)$ constant, and treat $St, \gamma$ and $\lambda$ as ${O}(1)$ constants.

We seek to write the system (3.38) in terms of $\epsilon$. We first consider (3.38h), use the integral mean value theorem and assume that $a(t)={O}(1)$ to obtain

(4.30)\begin{equation} f(t)=\frac{a(t)}{\displaystyle\int_{z=-a(t)}^{-a(t)+\epsilon c(t)} \exp(-\lambda(T(z,t)-1)) \, \textrm{d} z}=\frac{1}{\epsilon}\left( \frac{a(t)}{c(t)}+{O}(\epsilon)\right). \end{equation}

At the upper interface, we have

(4.31)\begin{equation} T(-a(t)+\epsilon c(t),t)=1-\epsilon, \end{equation}

which can be expanded to give

(4.32)\begin{equation} T(-a(t),t)+\epsilon c(t) T_z(-a(t),t)=1-\epsilon +{O}(\epsilon^{2}), \end{equation}

and hence

(4.33)\begin{equation} c(t)=-\frac{1}{T_z(-a(t),t)}, \end{equation}

providing $T_z(-a(t),t)$ is negative and ${O}(1)$. We thus have the leading-order system

(4.34a)\begin{gather} \frac{\partial T}{\partial t}+T_z(-a(t),t)a(t)\frac{ \widehat{Pe}}{ St}\frac{\partial T}{\partial z}=\frac{1}{St}\frac{\partial^{2} T}{\partial z^{2}},\quad \mbox{for } -a(t) \leq z \leq 0, \end{gather}

(4.34b)\begin{gather}-\frac{\textrm{d} a}{\textrm{d} t}-\frac{\partial T}{\partial z}(-a(t),t)a(t)\frac{\widehat{Pe}}{St}=\gamma + \frac{\partial T}{\partial z}, \quad \mbox{at } z=a(t), \end{gather}

(4.34c)\begin{gather}T(-a(t),t)=1, \quad T(0,t)=0, \end{gather}

(4.34d)\begin{gather}T(z,0)=0, \quad a(0)=0. \end{gather}

We look for steady states as before, finding

(4.35)\begin{equation} T^{*}(z)=\frac{\exp(-T^{*}_z(-a^{*})a^{*}\,\widehat{Pe} z)-1}{\exp(T_z(-a^{*}){a^{*}}^{2}\widehat{Pe}^{\vphantom{A^A}})-1}, \end{equation}

which we can differentiate to find the implicit relationship for $T^{*}_z(a^{*})$

(4.36)\begin{equation} T^{*}_z(-a^{*})=-T^{*}_z(-a^{*})a^{*}\widehat{Pe} \frac{\exp(T^{*}_z(-a^{*}){a^{*}}^{2}\widehat{Pe})}{\exp(T_z(-a^{*}){a^{*}}^{2}\widehat{Pe}^{\vphantom{A^A}})-1}.\end{equation}

Since $T^{*}_z(a^{*}) \neq 0$ by assumption, we can rearrange (4.36) to obtain

(4.37)\begin{equation} T^{*}_z(-a^{*})=-\frac{\ln(1+a^{*}\widehat{Pe} )}{{a^{*}}^{2}\widehat{Pe}^{\vphantom{A^A}}}. \end{equation}

The Stefan condition (4.34b) gives

(4.38)\begin{equation} -T^{*}_z(-a^{*})a^{*}\frac{\widehat{Pe}}{St}=\gamma+T^{*}_z(-a^{*}), \end{equation}

and hence we can substitute in $T^{*}_z(a^{*})$ from (4.37) to obtain the implicit equation for $a^{*}(\gamma ,St,\widehat {Pe})$, i.e.

(4.39)\begin{equation} \frac{\ln (1+a^{*}\widehat{Pe})}{a^{*}\widehat{Pe}^{\vphantom{A^A}}}=\frac{\gamma}{1+a^{*}\dfrac{\widehat{Pe}^{\vphantom{A^A}}}{St}}. \end{equation}

In our tall, thin model the averaged convection comes from $Pe \,f(t)$, with the downwards motion under gravity caused by the weight of the fluid, resisted by the viscous layer. When the viscous layer is thin, i.e. there is a thin crust, the weight of the fluid will be too strong and the charge will undergo free fall. However, if the thermal convection is slow, due to a small Péclet number, then the gravitational and viscous effects are in balance, $Pe \,f(t) ={O}(1)$, and the system is able to run to a steady state.

4.2.2. Case B: large heating, $\gamma \gg 1$

Suppose instead that $Pe={O}(1)$, but that $\gamma \gg 1$, meaning that the radiative heating into the fluid is very large. Again we keep $\lambda$ and $St$ as ${O}(1)$ constants. We look for balances such that the system tends to a steady state, and so examine the steady state system (4.1). A possible balance from the Stefan condition (4.1b) is $f={O}(\gamma )$. However, if $f \gg 1$, then (4.1a) and (4.1c) cannot be satisfied unless there is a boundary layer of size ${O}(1/f)={O}(1/\gamma )$. In order to allow matching with the outer solution, the boundary layer has to be at the end near $z=-a^{*}$. In this boundary layer $T_z={O}(\gamma )$, and hence all three terms balance in the Stefan condition. We thus consider the scalings

(4.40)\begin{equation} f^{*}=\gamma \hat{f}, \quad z= -a^{*} + \frac{1}{\gamma}y, \end{equation}

so that in the boundary layer we have

(4.41)\begin{equation} T_{inner} (y)=\exp(-\hat{f}Pe\, y ), \end{equation}

where we have applied the boundary condition $T(-a^{*})=1$ and matched with the outer solution $T_{\textrm {outer}}\equiv 0$. Note that the inner solution is thus the composite solution for the whole region, so

(4.42)\begin{equation} T(z)=\exp(-\hat{f}Pe\, \gamma (z+a^{*}) ). \end{equation}

The Stefan condition (4.1b) thus gives

(4.43)\begin{equation} \hat{f} \gamma\frac{Pe}{St}=\gamma -\hat{f}Pe\, \gamma, \end{equation}

so that

(4.44)\begin{equation} \hat{f}=\frac{1}{Pe}\frac{St}{1+St}. \end{equation}

Applying the boundary condition $T(-b^{*})=1-\epsilon$, we find

(4.45)\begin{equation} a^{*}-b^{*}=-\left(1+\frac{1}{St}\right)\frac{1}{\gamma}\ln(1-\epsilon).\end{equation}

Since $\lambda ={O}(1)$, we have that

(4.46)\begin{equation} f^{*}=\gamma \hat{f}=\frac{a^{*}}{a^{*}-b^{*}}, \end{equation}

which is the reciprocal of the ratio of the thickness of the viscous crust to the total thickness of the charge (including the crust), and hence, combining (4.44) and (4.45) with (4.46), the thickness of the total charge is

(4.47)\begin{equation} a^{*}=-\frac{1}{Pe}\ln(1-\epsilon)=\frac{\epsilon}{Pe}+{O}(\epsilon^{2}), \end{equation}

while the relative thickness of the crust to the total charge is

(4.48)\begin{equation} \frac{a^{*}-b^{*}}{a^{*}}=\frac{1}{\gamma}Pe \left(1+\frac{1}{St}\right), \end{equation}

which is size ${O}(1/\gamma )$. Thus, we have a thin crust for all $\gamma \gg 1$. Notice that these results have been derived with $\gamma$ independent of $\epsilon$, so hold regardless of the relative sizes of $\epsilon$ and $\gamma ^{-1}$.

4.2.3. Case C: small Stefan number, $St \ll 1$

Suppose that the Stefan number is small, so that we have $0 < St \ll 1$. This could correspond physically to very large latent heat, or to very small specific heat. We rescale our problem onto the fast time scale $t=St \hat {t}$, and for a balance we must have large heating, i.e. $\gamma =\hat {\gamma }/St$. With these scalings we obtain the system

(4.49a)\begin{gather} \frac{\partial T}{\partial \hat{t}}-{f(\hat{t}) Pe}\frac{\partial T}{\partial z}=\frac{\partial^{2} T}{\partial z^{2}}, \quad \mbox{in } -a(\hat{t}) \leq z \leq 0, \end{gather}

(4.49b)\begin{gather}-\frac{\textrm{d} a}{\textrm{d} \hat{t}}+{f(\hat{t}) Pe}={\hat{\gamma}} + St \frac{\partial T}{\partial z}, \quad \mbox{at } z=-a(\hat{t}), \end{gather}

(4.49c)\begin{gather}T(-a(\hat{t}),\hat{t})=1, \quad T(-b(\hat{t}),\hat{t})=T_b, \quad T(1,\hat{t})=0, \end{gather}

(4.49d)\begin{gather}T(z,0)=0, \end{gather}

(4.49e)\begin{gather}f(\hat{t})=\frac{a(\hat{t})}{\displaystyle\int_{z=-a(\hat{t})}^{-b(\hat{t})}\exp\left(-\lambda(T(z,\hat{t})-1)\right) \, \textrm{d} z}, \end{gather}

(4.49f)\begin{gather}a(0)=1, \quad b(0)=1. \end{gather}

In the limit $St \to 0$, we obtain the usual system, but without the conductive flux in the Stefan condition (4.49b). Looking for a steady-state solution, we find that from (4.49b) we must have $\hat {\gamma }=Pe \,f^{*}$, and hence the heating from below and the convective term are identically in balance. The intuition behind this statement is that if the heating is too high, or convection too small, all the charge will be melted and the domain will shrink. Conversely, if the heating is too small, or convection too fast, the charge will undergo free fall.

We find the steady-state temperature

(4.50)\begin{equation} T^{*}(z)=\frac{\exp(-f^{*}Pe\, z)-1}{\exp(\,f^{*}Pe\, a^{*})-1}=\frac{\exp(-\hat{\gamma} z)-1}{\exp(\hat{\gamma} a^{*})-1}, \end{equation}

and thus obtain the integral equation for $f^{*}$, as follows:

(4.51)\begin{equation} \frac{a^{*}}{\displaystyle\int_{z=-a^{*}}^{-b^{*}}\exp( -\lambda (T^{*}(z)-1)) \, \textrm{d} z}=\frac{\hat{\gamma}}{Pe}.\end{equation}

Note that we cannot determine $a^{*}$ from the Stefan condition (4.49b), since the conductive flux does not feature. Instead, we can rearrange (4.51) using algebraic steps analogous to those used in § 4.1. Changing the variable $a^{*}$ to

(4.52)\begin{equation} \varTheta = \frac{\exp(\hat{\gamma}a^{*})}{\exp(\hat{\gamma}a^{*})-1},\end{equation}

we obtain

(4.53)\begin{equation} \textrm{e}^{\lambda \varTheta}(E_1( \lambda \varTheta) - E_1 (\lambda (\varTheta - 1+T_b) ) )=\frac{Pe}{\hat{\gamma}}\ln \left(\frac{\varTheta}{\varTheta-1} \right). \end{equation}

Upon finding a root of (4.53), we then obtain

(4.54a)\begin{gather} a^{*}=\frac{1}{\hat{\gamma}}\ln\left(\frac{\varTheta}{\varTheta-1} \right), \end{gather}

(4.54b)\begin{gather}b^{*}=\frac{1}{\hat{\gamma}}\ln\left(1+\frac{T_b}{\varTheta-1} \right). \end{gather}

We remark that while this approach may appear similar to the method used to find general steady states in § 4.1, the ordering of algebraic steps is different. This is because we cannot define the change of variables as in (4.7), since writing $\theta =(1/St)({\hat {\gamma }}/{f^{*} Pe}-1)$ is not compatible with the condition $\hat {\gamma }=f^{*} Pe$. Instead, we first take the limit $St \to 0$, and then subsequently make the change of variables (4.52).

4.2.4. Case D: large Stefan number, $St \gg 1$

The case of a large Stefan number, $St \gg 1$, could correspond to the limit of small latent heat, where the fluid drips off without any vaporisation due to phase-change or chemical reactions.

If a steady state exists, then it can be explored by taking the limit $St \to \infty$ in (4.4) and (4.5), giving

(4.55a)\begin{gather} a^{*}=\frac{1}{f^{*}Pe}\ln\left(\frac{\gamma/f^{*}Pe}{\gamma/f^{*}Pe-1}\right), \end{gather}

(4.55b)\begin{gather}b^{*}=\frac{1}{f^{*}Pe}\ln \left(1+\frac{T_b}{\gamma/f^{*}Pe -1} \right). \end{gather}

Here, $f^{*}$ can be found from the root of (4.15), with $St \to \infty$ and $\theta \to \gamma /f^{*}Pe$.

We note from (4.55b) that we require $f^{*}Pe/\gamma <1$ for a steady state to exist. Since the integral term $f^{*} \geq 1$, no steady states exist if $Pe/\gamma >1$. For the simplified case of $\lambda =0=T_b$, on the bifurcation plot in figure 10 we can see that for large $St$ all that matters is the value of $Pe/\gamma$. This ratio is given in terms of physical dimensional parameters as

(4.56)\begin{equation} \frac{Pe}{\gamma}=\frac{{Uh}/{\kappa}}{{\epsilon_R \sigma (T_h^{4}-T_d^{4})}h/{k(T_d-T_c)} }=\frac{U \rho C_p (T_d - T_c)}{\epsilon_R \sigma (T_h^{4}-T_d^{4})}. \end{equation}

Using the Poiseuille velocity scaling $U=({\rho g}/{8 \eta _s})R^{2}$, we can further write

(4.57)\begin{equation} \frac{Pe}{\gamma}=\frac{\rho^{2} g R^{2} C_p (T_d - T_c)}{8 \eta_s \epsilon_R \sigma (T_h^{4}-T_d^{4})}. \end{equation}

We see that a higher density or lower viscosity will lead to a higher value of $Pe/\gamma$ in (4.57), if all other quantities remain the same. This could explain why a continuous downwards flow is often experienced in ferro-silicon, but not silicon, production in practical furnace operations. Of course, the other parameters will not be the same in both settings, and so specific values should be used when comparing these processes.

4.2.5. Case E: large heating and small convection

Suppose $\gamma \gg 1, Pe \ll 1, St=(Pe/\gamma ) \widehat {St}$, with $\lambda ={O}(1)$ and $1-T_b={O}(1)$. Then we have the steady state $T^{*}(z)=-z/a^{*}$ and $b^{*}=T_b a^{*}$, and can integrate the denominator in $f^{*}$ to find $f^{*}=\lambda /(\exp (\lambda (1-T_b))-1 )$. There is no way to find $a^{*}$ from the Stefan condition, unless $a^{*}={O}(\gamma ^{-1})$. We thus set $a^{*}=\gamma ^{-1} \hat {a}$, and from the steady-state Stefan condition obtain $\hat {a}=\widehat {St}/(\widehat {St}-f^{*})$, and hence

(4.58)\begin{equation} a^{*}=\frac{1}{\gamma}\frac{\widehat{St}}{\widehat{St}^{\vphantom{A^A}}-\lambda/\exp(\lambda(1-T_b) )}, \end{equation}

with $a^{*}-b^{*}=(1-T_b)a^{*}$. This corresponds physically to the case of a large latent heat, $L$, and a high source temperature $T_h$.

4.2.6. Case F: no convection

Suppose $f\,Pe \ll 1$, for example if the material is very light or highly thermally conductive. Then, for all other parameters being ${O}(1)$, we recover a problem similar to the canonical Stefan problem, without convection, but with a radiative flux, namely

(4.59a)\begin{gather} \frac{\partial T}{\partial t}=\frac{1}{St}\frac{\partial^{2} T}{\partial z^{2}},\quad \mbox{for } -a(t) \leq z \leq 0, \end{gather}

(4.59b)\begin{gather}-\frac{\textrm{d} a}{\textrm{d} t}=\gamma + \frac{\partial T}{\partial z}, \quad \mbox{at } z=-a(t), \end{gather}

(4.59c)\begin{gather}T(-a(t),t)=1, \quad T(-b(t),t)=T_b, \quad T(0,t)=0, \end{gather}

(4.59d)\begin{gather}T(z,0)=0, \end{gather}

(4.59e)\begin{gather}a(0)=1, \quad b(0)=1. \end{gather}

The system (4.59) admits the steady-state solution

(4.60)\begin{equation} T^{*}(z)=-\frac{z}{a^{*}}, \quad a^{*}=\frac{1}{\gamma}, \quad b^{*}=\frac{T_b}{\gamma}, \end{equation}

and hence the ratio of the thickness of the crust to the total charge is given by

(4.61)\begin{equation} \frac{a^{*}-b^{*}}{a^{*}}=1-T_b. \end{equation}

Since convection did not enter the problem, the presence of the upper interface $z=-b(t)$ is artificial, since it has no influence on the system, but merely marks where we believe the crust meets the rest of the charge.

5.1. Numerical simulations in time

It is natural to use a finite volume scheme to obtain numerical solutions of our one-dimensional problem. To account for the moving lower interface, we map the problem from the moving domain $z\in [-a(t),0]$ onto the fixed domain $x\in [0,1]$, with the transformation $x=({z}/{a(t)})+1, \tau =t$. In addition to the moving interfaces, potential difficulties come from the non-local convective term $f(t)$, and to perform the necessary quadrature we use an in-built trapezium rule implementation in Matlab. When finding the temperature from the discretised version of (3.38a) in the transformed coordinates, we use a scheme that is implicit in $T$ but explicit in the coefficients of $T$ derivatives. Details of our numerical scheme are provided in appendix A. We use a spatial step of ${\rm \Delta} x=4\times 10^{-3}$ and a time step of ${\rm \Delta} t=10^{-4}$, unless otherwise stated.

In order to demonstrate a sample dynamics, we plot the evolution of temperature and velocity in the original domain for one set of parameters, in figure 6. For the parameter values used, the system tends towards a steady state.

Figure 6. Heat maps of temperature profiles $T$ (a–c) and vertical velocity component $u_z$ (d–f) at times $t=0.01$ (a,d), $t=0.1$ (b,e) and $t=1$ (c,f). In the temperature heat maps we illustrate the location of the interfaces $z=-a(t)$ and $z=-b(t)$ by solid and dashed lines, respectively. Parameter values used are $St=0.1, Pe=1.2, \lambda =0.1, T_b=0.1$ and $\gamma =20$. (a) $T$ at $t=0.01$, (b) $T$ at $t=0.1$, (c) $T$ at $t=1$, (d) $u_z$ at $t=0.01$, (e) $u_z$ at $t=0.1$ and (f) $u_z$ at $t=1$.

We are mostly interested in the position of the interfaces, with temperature being an intermediate variable to this end. We thus show some example trajectories of the distances $a(t)$ and $b(t)$ from the top of the domain, where we vary one parameter but keep the other four fixed. In figure 7 we vary the dimensionless condensation temperature, $T_b$. Note that for some parameter values the system tends towards a steady state, but for larger $T_b$ both interfaces move downwards to $z=-\infty$. The relative crust thickness, $(a-b)/a$ tends to a positive value when the system runs to a stable steady state, but when the system is in free fall it tends to zero, since $a$ and $b$ become unbounded. For early time $Q(t)$ has a large negative spike due to the initial downwards motion of the interface $z=-a(t)$, caused by the discontinuity between the initial condition and the lower boundary condition for the temperature. As the temperature at the base increases, due to radiative heating from below, the interface rises. In the case of the steady-state cases $T_b=0$ and $T_b=0.1$ the flux tends to $Q^{*}=\pi f^{*}$, as given in (3.40). The pilot furnace experiments discussed in Sloman et al. (Reference Sloman, Please, Van Gorder, Valderhaug, Birkeland and Wegge2017) take 80 minutes, which roughly corresponds to the dimensionless time for the interfaces to reach steady state in figure 7.

Figure 7. Plots showing the influence of varying the parameter $T_b$ for (a) the position of the interfaces $z=-a(t)$ in solid lines and $z=-b(t)$ in dashed lines, (b) the relative crust thickness $(a(t)-b(t))/a(t)$, (c) the volumetric material flux $Q(t)$, given by (3.40), all over time. Other parameter values used are $St=0.1, Pe=1.2, \lambda =0.1$ and $\gamma =20$.

In figure 8 we vary the Péclet number $Pe$. For small values of the Péclet number, the system tends towards a steady state, but the interfaces undergo free fall for larger values of the Péclet number, since thermal convection is strong enough to overcome the radiative heating. As with figure 7, we find that the relative crust thickness tends to a positive value when the system is in steady state, but tends to zero when the system is in free fall. Similarly, we get a positive flux $Q(t)$ after an initial transient period, when the system runs to steady state. In the model we have set initial conditions for the interfaces as $a(0)=b(0)=1$. Despite the potential numerical problem from having an initially singular viscous domain, the simulations seem to cope, with the interfaces separating in a very fast transient period, as demonstrated in figures 7 and 8.

Figure 8. Plots showing the influence of varying the parameter $Pe$ for (a) the position of the interfaces $z=-a(t)$ in solid lines and $z=-b(t)$ in dashed lines, (b) the relative crust thickness $(a(t)-b(t))/a(t)$, (c) the volumetric material flux $Q(t)$, given by (3.40), all over time. Other parameter values used are $St=0.1, \lambda =0.1, \gamma =20$ and $T_b=0.1$.

Our numerical simulations suggest that the steady-state solution corresponding to the larger $\theta$ root from (4.14) is stable, while the steady state corresponding to the smaller $\theta$ root is unstable. In figure 9 we plot trajectories starting from steady states corresponding to the circled points shown in figure 5. In each case we set the initial temperature profile to be $T^{*}(z)$ and $f(0)=f^{*}$. All values used were to five decimal places. The upper trajectories (solid lines) return to the steady state, but the lower trajectories (dashed lines) move away from their respective steady states, eventually undergoing free fall. This explains why at most one steady state was found from numerical solutions of the model in § 5.1. Whenever we have two steady states, we find from numerical simulations that it is the steady state with larger $a$ and $b$ which is unstable and drops.

Figure 9. Plot showing the evolution of trajectories from a stable steady state (solid) and unstable steady state (dashed) for $z=-a(t)$ (red) and $z=-b(t)$ (blue). Parameter values used are $Pe=1, St=2, \gamma =30, T_b=0.8$ and $\lambda =1$. A spatial step size of ${\rm \Delta} x= 4 \times 10^{-3}$ and a time step of ${\rm \Delta} \tau =10^{-6}$ are used.

Numerical simulations of the model thus show that either the interfaces and temperature tend to a steady state, or the interfaces drop indefinitely, and so the charge experiences free fall. When there is free fall, this is not completed in finite time, but the interfaces eventually move at an approximately constant velocity. The steady states shown in figures 6–8 appear to be stable, in as much they can be obtained numerically for different initial conditions, within a given region of attraction. We now turn our attention to examining how the steady states vary with the model parameters.

5.2. Long-time steady-state dynamics

In figures 7 and 8 we showed example trajectories that tend to a stable steady state, and others which underwent free fall. It is of interest to determine where in the parameter space stable steady states will occur, and so here we present bifurcation diagrams illustrating where stable steady states can be found.

We first consider the simplified case $\lambda =0$ and $T_b=0$. This corresponds to the upper interface fixed at the top of the domain, $z=0$, and the whole charge being filled with a viscous fluid of constant viscosity. In this limit, we have $f^{*}=1$, with the gravity and viscous terms in balance, and the material flux given by $Q^{*}=\pi$ from (3.40). For any values of $Pe, St$, and $\gamma$ the location of $a^{*}$ can be determined from (4.4). We plot how $a^{*}$ depends on $Pe$ and $St$ for fixed $\gamma$, in figure 10. Where a steady state exists, its value is given by the heat map provided. The region where steady states cannot exist is labelled free fall. Note that there is a thin boundary layer where the steady state $a^{*}$ can be arbitrarily large, although we have just plotted steady states to a value of $a^{*}=2$, due to difficulties in numerical resolution. The interface between steady state and free fall in parameter space is given by $\theta =1$. Using (4.7), with $f^{*}=1$, we find that this corresponds to $Pe(1+({1}/{St}))=\gamma$, which asymptotes to $Pe=\gamma$ as $St$ becomes large. Thus, for large $St$, if the Péclet number is sufficiently large then the charge will keep on falling under gravity.

Figure 10. Heat maps of the value of the steady state $a^{*}$, comparing $Pe$ and $St$, for fixed $\gamma =8, \lambda =0$ and $T_b=0$.

We also consider the general case $\lambda >0$ and $T_b>0$. To find steady-state solutions we must first find zeros of (4.15) and then substitute the solution $\theta ^{*}$ into (4.8). In figure 11 we show heat maps of $a^{*}$, the material flux $Q^{*}$, the crust thickness $a^{*}-b^{*}$ and the relative crust thickness $(a^{*}-b^{*})/a^{*}$. The necessary root finding was achieved using a inbuilt solver in Matlab. We see that the steady state $a^{*}$ is greatest away from the interface between steady state and free fall. There is a narrow region for small $Pe$ where the flux $Q^{*}$ is very large. Although the fluid is falling slowly, since $a^{*}$ is small the top of the granular region has to be refilled rapidly, causing a high flux through the interface $z=-a^{*}$. Both the crust thickness and relative crust thickness are fairly constant in the parameter regime considered, although the crust thickness decreases as the Stefan number increases and the relative crust thickness decreases as the Péclet number increases. As with figure 10, there is difficulty with numerical resolution close to the boundary in parameter space between steady state and free fall.

Figure 11. Heat maps of the steady-state values of $a^{*}$ (a), $Q^{*}$ (b), $a^{*}-b^{*}$ (c) and $(a^{*}-b^{*})/a^{*}$ (d). Parameter values used are $\gamma =20, \lambda =0.2$ and $T_b=0.8$. (a) Interface position $a^{*}$. (b) Volumetric material flux $Q^{*}$. (c) Crust thickness $a^{*}-b^{*}$. (d) Relative crust thickness $(a^{*}-b^{*})/a^{*}$.

We can also determine where the steady states $a^{*}$ and $b^{*}$ occur for general $T_b$ and $\lambda$, and thus where the system undergoes free fall. We know that steady states only exist for $\theta >1$, and that as $\theta \searrow 1$ the value of $a^{*}$ becomes arbitrarily large. Since the interface is given by $\theta =1$, we numerically find the root of (4.15) with a value of $\theta$ slightly greater than $\theta =1$, showing us where the edge of the steady-state region is. This allows us to draw bifurcation diagrams, demonstrating where a stable steady state or free fall will occur. We do not plot the value of the steady state $a^{*}$, when it occurs, but it will decrease as distance from the bifurcation interface increases. In figure 12 we show bifurcation plots of $Pe$ against $St$ for different $\lambda$ and for different $T_b$, and in figure 13 we show plots of $Pe$ against $T_b$ for different $\gamma$ and for different $\lambda$. These plots have been obtained using Mathematica, with a value of $\theta =1+10^{-6}$, which is sufficiently close to $\theta =1$.

Figure 12. Plots showing the interface between where the solutions tend to steady states $a^{*}$ and $b^{*}$ and where the fluid undergoes free fall, with $\gamma =8$ in both cases. Plots were generated in Mathematica by letting $\theta =1+10^{-6}$. In (a) we vary $\lambda$, with values shown $\lambda =10^{-6},0.1,0.2,0.3,0.4$ and $T_b=0$. In (b) we vary $T_b$, with values shown $T_b=10^{-1},10^{-2},10^{-3},10^{-4},10^{-5},0$ and $\lambda =0.1$.

Figure 13. Plots showing the interface between where the solutions tend to steady states $a^{*}$ and $b^{*}$ and where the fluid undergoes free fall, with $St=1$ in both cases. Plots were generated in Mathematica by letting $\theta =1+10^{-6}$. In (a) we vary $\gamma$, with values shown $\gamma =10,20,30,40,50$ and $\lambda =0.1$. In (b) we vary $\lambda$, with values shown $\lambda =10^{-6},1,2,3,4,5$ and $\gamma =30$.

We can also consider how the thickness of the crust is influenced by the dimensionless parameters in the model by plotting the result of parametric sweeps of the Péclet number against $\gamma , St, T_b$, and $\lambda$ in figure 14. For each curve we found the steady-state value by determining the root of (4.15) for 100 values of $Pe$ in the range plotted. For some parameters steady states cease to exist when the Péclet number is above a certain value, as shown in figures 12 and 13. In these cases we have only drawn the curve up to that value. In each case we see that the crust thickness increases with $Pe$. As seen in figures 8, 10 and 11, $a^{*}$ increases with the Péclet number, since material falls at a faster rate and so the crust will form lower down. This could explain why the crust increases in size, since $b^{*}$ would not increase by more than $a^{*}$. Conversely, as $\gamma$ increases the thickness of the crust decreases, since higher heating from beneath the liquid causes $a^{*}$ to decrease. A higher value of the Stefan number also reduces the crust thickness, since there will be less latent heat at the interface $z=-a^{*}$ and so the radiative heating will cause more material to be consumed and the interface to be driven higher, narrowing the crust. We might expect that a higher value of $T_b$ would lead to a thinner crust, since the crust exists over a narrower temperature range, due to condition (4.1d). However, the opposite is true, as shown in figure 14. This is due to condition (4.1e), whereby a larger crust is required to balance the additional weight that comes from the upper region $z>-b^{*}$. Similarly, larger $\lambda$ leads to a smaller steady-state crust, since the crust becomes more viscous and does not need to be as thick to balance the weight from above.

Figure 14. Plots of the crust thickness $a^{*}-b^{*}$ against $Pe$ using parameter sweeps in (a) $\gamma$, (b) $St$, (c) $T_b$, (d) $\lambda$. When not stated in the plot, fixed parameter values are $St=1, \gamma =10, \lambda =1$ and $T_b=0.5$.