2.1. Boundary and initial conditions

The equations of conservation of mass in the solid and gas phases (2.6a)–(2.6b) are both hyperbolic. We impose the inflowing concentration of solid at $x=x_0$, and a flux of gas $f_*$ through the free boundary, so that

(2.8a)\begin{equation} C=C^{in} \quad \text{at }x=x_0, \quad f=f_* \quad \text{at }x=s(t). \end{equation}

We prescribe a cold temperature of the inflowing solid material

(2.8b)\begin{equation} T_s=T_s^{in} \quad \text{at }x=x_0. \end{equation}

At the free boundary, we impose

(2.8c)\begin{equation} T_g=T_s=\frac{R}{s} \quad \text{at }x=s(t),\end{equation}

which is a crude model for a point radiative heat source at $x=0$, with constant intensity $R$, incident on this boundary. This boundary condition is motivated by the electric arc which forms in the (roughly cylindrical) crater region within the material bed of a silicon submerged arc furnace (Schei et al. Reference Schei, Tuset and Tveit1998). For other applications, different forms for the boundary condition (2.8c) may be more appropriate.

The position of the free boundary $x=s(t)$ is determined by the Stefan-like condition on the concentration of solid material

(2.8d)\begin{equation} C=C_* \quad \text{at }x=s(t), \end{equation}

where $C_*< C^{in}$. Thus the free boundary of the domain is the position where sufficient solid material has been reacted away. For silicon furnaces, this condition describes the position of the edge of the gas-filled crater that forms deep within the furnace: $C_*$ is the critical concentration at which the solid structure loses structural integrity, and any remaining solid material collapses into the crater.

Finally, we require an initial configuration for the system as initial conditions on $C$, $f$, $T_s$, $T_g$ and the position $s$ of the free boundary. These initial conditions must be consistent with the boundary conditions described above.

2.2. Non-dimensionalisation

We non-dimensionalise the model (2.6) and (2.8) by scaling $x=x_0 x'$, $s=x_0 s',$ $t=x_0t'/u_s,$ $C=C^{in}C',$ $f=[f]f',$ $T_s=[T_s]T_s'$ and $T_g=[T_g]T_g'$, where the prime notation denotes dimensionless variables. The temperature scalings are given by the boundary condition at $x=s(t)$, since this is the major source of heat to the system, and so we fix $[T_s]=[T_g]=R/[s],$ where $[s]$, to be determined, is the expected size of $s$. We define $\delta$ to be the inverse Péclet number

(2.9)\begin{equation} \delta=Pe^{{-}1}=\frac{k_s}{x_0u_sMc_{p,s}C^{in}}. \end{equation}

Throughout this paper we assume that $\delta$ is small, so that the effective thermal conductivity of the porous solid structure is small over the relatively large length scale of the domain. Thus the heat provided to the solid material at the boundary $x=s(t)$ cannot penetrate far into the domain: we anticipate a conductive boundary layer at this boundary, of width $O(\delta )$.

Assuming that the majority of the chemical reaction occurs in this hot $O(\delta )$ boundary layer at $x=s(t)$, the length scale $[s]$ is given by a balance in the solid mass equation (2.6a),

(2.10)\begin{equation} \frac{u_s}{\delta x_0}=k_{Q}\exp\left(\frac{-a_Q[s]}{R}\right), \end{equation}

which, along with the definition of $\delta$, gives $[s]=(R/a_Q)\log (k_{Q}k_s/(Mc_{p,s}u_s^{2}C^{in})).$ The temperature scalings are therefore $[T_s]=[T_g]=a_Q/\log (k_Qk_s/(Mc_{p,s}u_s^{2} C^{in})).$ We note that for these temperature scalings to be positive, the velocity $u_s$ of the solid material must be sufficiently small. Specifically, we require

(2.11)\begin{equation} u_s<\sqrt{\frac{k_Qk_s}{Mc_{p,s}C^{in}}}. \end{equation}

For larger $u_s$ than this, the (bounded) Arrhenius chemical reaction rate cannot match the rate at which material flows in, and the model ceases to be valid. Finally, balancing the gas mass equation (2.6b) requires $[f]=u_sC^{in}$.

With these scalings, and dropping the prime notation, the equations become

(2.12a)$$\begin{gather} \delta\left(\frac{\partial{C}}{\partial{t}}-\frac{\partial{C}}{\partial{x}}\right)={-}Q, \end{gather}$$

(2.12b)$$\begin{gather}\delta\frac{\partial{f}}{\partial{x}}=Q, \end{gather}$$

(2.12c)$$\begin{gather}C\left(\frac{\partial{T_s}}{\partial{t}}-\frac{\partial{T_s}}{\partial{x}}\right)=\delta\frac{\partial{{}^{2}T_s}}{\partial{x^{2}}}+\mu(T_g-T_s)-\frac{\gamma}{\delta} Q, \end{gather}$$

(2.12d)$$\begin{gather}K\frac{\partial}{\partial x}(fT_g)=\frac{1}{\delta}T_sQ-\mu(T_g-T_s), \end{gather}$$

on the domain $x\in (s(t),1)$, where the dimensionless chemical reaction rate is

(2.12e)\begin{equation} Q(C,T_s)=C\exp\left(A\left(1-\frac{1}{T_s}\right)\right), \end{equation}

and with boundary conditions

(2.12f)$$\begin{gather} C=C_*,\quad f=f_*, \quad T_s=T_g=\frac{\rho}{s} \quad \text{at }x=s(t), \end{gather}$$

(2.12g)$$\begin{gather}C=1,\quad T_s=T^{in} \quad \text{at }x=1, \end{gather}$$

with $f_*,C_*$ and $T^{in}$ appropriately rescaled.

Here, we have introduced the dimensionless parameters

(2.13a–e)\begin{equation} A=\frac{a_Q[s]}{R}, \quad \mu=\frac{\lambda x_0}{u_sMc_{p,s}C^{in}}, \quad K=\frac{c_{p,g}}{c_{p,s}}, \quad \gamma=\frac{[s]\Delta H}{Mc_{p,s}R}, \quad \rho=\frac{[s]}{x_0}. \end{equation}

In addition to (2.12), we require initial conditions for $C$, $f$, $T_g$, $T_s$ and $s$. The system (2.12) is a fifth-order system of differential equations in space, with an additional degree of freedom in the position $s$ of the free boundary. We therefore expect the model to be correctly specified since we impose six boundary conditions.

Throughout, we assume that the dimensionless parameters $K,$ $\gamma,$ $A,$ $\rho$, $C_*$ and $f_*$ in the model (2.12) are of order unity (see Luckins et al. (Reference Luckins, Oliver, Please, Sloman and Van Gorder2022) for justification in the case of a silicon furnace). The inverse Péclet number, $\delta$, is expected to be of the order of $10^{-3}-10^{-2}$ (Luckins et al. Reference Luckins, Oliver, Please, Sloman and Van Gorder2022), and is the smallest parameter in the system. The asymptotic analysis in this paper is in the limit of $\delta \rightarrow 0$. We also assume that the input temperature of the solid, $T^{in}$, is small. Specifically, we require

(2.14)\begin{equation} \delta T^{in}\gg \exp\left(A\left(1-\frac{1}{T^{in}}\right)\right), \end{equation}

so that the rate of reaction at this temperature is much less than the rate at which material moves through the whole domain by advection. Thus, if there were no heating, we would expect the solid material not to appreciably react away over the time it spends moving through the domain. We further assume that $\delta \ll T^{in}$, meaning that, although $T^{in}$ is small in the sense of (2.14), it is not as small as $\delta$. This further assumption that $\delta \ll T^{in}$ does not impact the dominant behaviour of the system, and is only needed when analysing the correction terms to the leading-order behaviour in certain parameter regimes. The final dimensionless parameter is $\mu$, which encodes the rate of convective heat transfer between the phases. We will explore various distinguished limits relating $\mu$ to $\delta$.

4.1. Region AII

We begin by considering the behaviour over an $O(1)$ length scale, which we denote region AII in figure 4. We recall that $T^{in}\ll 1$ is small, in the specific sense (2.14). At the boundary $x=1$ we have $T_s= T^{in}$, and we therefore look for a solution expansion of the form

(4.2a,b)\begin{equation} T_s=T^{in}\left(\hat{T}^{0}_s+o\left(T^{in}\right)\right), \quad C=\hat{C}^{0}+o\left(T^{in}\right), \end{equation}

and also expand the gas variables in region AII by

(4.3a,b)\begin{equation} T_g=\hat{T}^{0}_g+o\left(T^{in}\right), \quad f=\hat{f}^{0}+o\left(T^{in}\right). \end{equation}

(Here, we use hat notation to differentiate from the variables used in region AI below.) Since by (2.14) we have that $Q(T^{in})\ll \delta T^{in}$, at leading order in (2.12a) and (2.12c) we obtain

(4.4a,b)\begin{equation} \frac{\partial{\hat{C}^{0}}}{\partial{t}}-\frac{\partial{\hat{C}^{0}}}{\partial{x}}=0, \quad \hat{C}^{0}\left(\frac{\partial{\hat{T}_s^{0}}}{\partial{t}}-\frac{\partial{\hat{T}^{0}_s}}{\partial{x}}\right)=0. \end{equation}

Solving, using the boundary conditions at $x=1$, we must have, for all $x$

(4.5a,b)\begin{equation} \hat{T}_s^{0}=1, \quad \hat{C}^{0}=1. \end{equation}

Since $\hat {T}_s$ and $\hat {C}$ are constant to leading order, the gas variables in this region, given by (4.1), are also uniform. The error in these leading-order solutions is $O(\delta ^{-1}Q(T^{in}))=o(T^{in})$, since $Q\ll \delta T^{in}$.

4.2. Region AI

We note that we cannot impose the boundary conditions at the free boundary $x=s(t)$ on the constant solution found in the previous section. We therefore look for a boundary layer at $x=s(t)$, denoted region AI, where $T_s=O(1)\gg T^{in}$. Since $\delta \ll 1$, we change variables, setting $x=s(t)+\delta X$. Since $T^{in}\gg \delta$, we look for an expansion in powers of $T^{in}$ in region AI, setting

(4.6a–c)\begin{align} T_s=T_s^{0}(X,t)+O\left(T^{in}\right), \quad C=C^{0}(X,t)+O\left(T^{in}\right), \quad s=s^{0}(t)+O\left(T^{in}\right). \end{align}

We also expand the gas variables in powers of $T^{in}$, taking

(4.7a,b)\begin{equation} T_g=T_g^{0}(X,t)+O\left(T^{in}\right), \quad f=f^{0}(X,t)+O\left(T^{in}\right). \end{equation}

Substituting the expansions (4.6a–c) into (2.12a) and (2.12c), at leading order we obtain

(4.8a)$$\begin{gather} -\left(\dot{s}_0+1\right)\frac{\partial{C^{0}}}{\partial{X}}={-}C^{0}\exp\left(A\left(1-\frac{1}{T_s^{0}}\right)\right), \end{gather}$$

(4.8b)$$\begin{gather}-\left(\dot{s}_0+1\right)C^{0}\frac{\partial{T_s^{0}}}{\partial{X}}=\frac{\partial{{}^{2}T_s^{0}}}{\partial{X^{2}}}-\gamma C^{0}\exp\left(A\left(1-\frac{1}{T_s^{0}}\right)\right), \end{gather}$$

where the dot notation denotes a time derivative, $\dot {s}_0=\mathrm {d} s_0/\mathrm {d}t$. From (2.12 f), we have the boundary conditions

(4.8c)\begin{equation} C^{0}=C_*, \quad T_s^{0}=\frac{\rho}{s_0}, \quad \text{at }X=0.\end{equation}

As $X\rightarrow \infty$ we match with the outer region AII, and so require

(4.8d)\begin{equation} C^{0}\rightarrow 1, \quad T_s^{0}\rightarrow 0 \quad \text{as}\ X\rightarrow\infty.\end{equation}

The far-field temperature must approach zero in (4.8d) so that the chemical reaction rate $Q$ approaches zero in this limit. This is necessary for the problem to be correctly specified. Region AI is the reaction region, with the energy for the reaction provided by the heat radiation on the boundary $x=s(t)$. This reduced system (4.8) in region AI determines $T_s^{0}$, $C^{0}$ and also gives a relation between $\dot {s}_0$ and $s_0$. We solve (4.8) numerically in § 10 below. We note that the leading-order system (4.8) is quasi-steady, in that the only time derivatives are $\dot {s}_0$, with no time derivatives of $T_s$ or $C_s$ at leading order. Since the region-AII problem has a constant, time-independent solution, we see that the entire system is quasi-steady in case A. The relation between $\dot {s}_0$ and $s_0$ found numerically as part of the solution of (4.8) in § 10 therefore describes fully the quasi-steady time-dependent behaviour of the system. In particular we may find steady state(s) of the system for which $\dot {s}=0$.

In fact, we obtain a quasi-steady boundary-layer problem, similar to (4.8), in each of the non-zero $\mu$ cases explored below. In cases with larger $\mu$, however, we will find fully time-dependent behaviour outside of region I, so that, even though the region-I behaviour is still quasi-steady, the overall system is not quasi-steady. In the analysis of many of these cases we find and discuss the time-dependent solutions, but when computing solutions numerically in §§ 10–11 we will restrict ourselves to the steady-state case, for simplicity.

Given the solution of this solid problem (4.8), from (4.1) with the change of variables $x=s(t)+\delta X$, we see that the leading-order gas variables in region AI are

(4.9a)$$\begin{gather} f^{0}(X,t)=f_*+(1+\dot{s}_0)(C^{0}-C_*), \end{gather}$$

(4.9b)$$\begin{gather}T_g^{0}(X,t)=\frac{\dfrac{f_*\rho}{s_0}+K^{{-}1}\int_{\bar{X}=0}^{X} T_s^{0}C^{0}\exp\left(A\left(1-\dfrac{1}{T^{0}_s}\right)\right)\,\mathrm{d}\bar{X}}{f_*+(1+\dot{s}_0)(C^{0}-C_*)}. \end{gather}$$

In region AII both $\hat {T}_g^{0}$ and $\hat {f}^{0}$ are uniform. By matching with region AI at $x=s_0$, we find

(4.10a)$$\begin{gather} \hat{f}^{0}(t)=F_s(t):=f_*+(1+\dot{s}_0)(1-C_*), \end{gather}$$

(4.10b)$$\begin{gather}\hat{T}_g^{0}(t)=G_s(t):=\frac{\dfrac{f_*\rho}{s_0}+K^{{-}1}\int_{\bar{X}=0}^{\infty} T_s^{0}C^{0}\exp\left(A\left(1-\dfrac{1}{T^{0}_s}\right)\right)\,\mathrm{d}\bar{X}}{f_*+(1+\dot{s}_0)(1-C_*)}. \end{gather}$$

4.3. Composite solution

We construct composite leading-order solutions for case A as follows. Since $C$, $f$ and $T_g$ are $O(1)$ in all regions of the domain, and uniform to leading order in region AII, the leading-order composite solutions in case A (which we denote by superscript A) are simply

(4.11a–c)\begin{equation} C^{A}=C^{0}\left(\frac{x-s_0(t)}{\delta},t\right), \quad f^{A}=f^{0}\left(\frac{x-s_0(t)}{\delta},t\right), \quad T_g^{A}=T_g^{0}\left(\frac{x-s_0(t)}{\delta},t\right). \end{equation}

However, $T_s$ is $O(1)$ in region AI, but $O(T^{in})$ in region AII. We note that this small temperature in region AII must induce a $O(T^{in})$ switchback correction (Hinch Reference Hinch1991) to the leading-order solution in the boundary layer AI. We therefore have some freedom in how to choose the leading-order composite expansion (Hinch Reference Hinch1991). For simplicity we choose to set the $O(T^{in})$ correction in region AI to be uniform, and equal to $T^{in}$. Thus we form the composite expansion

(4.12)\begin{equation} T_s^{A}=T_s^{0}\left(\frac{x-s_0(t)}{\delta}\right)+T^{in}, \end{equation}

which has $O(T^{in})$ error in region AI, and $O(\delta ^{-1}Q(T^{in}))$ error in region AII.

In case A, we have seen that all the material consumption (chemical reaction) occurs within the boundary layer AI at $x=s(t)$, and that all variables have constant, uniform profiles in the outer region AII. We also observed that the leading-order behaviour is quasi-steady throughout the domain, with the time-derivatives of $s$ in region AI the only time derivatives in the reduced problem.

5.1. Region BI

As in case A, we look for a reaction region at the free boundary $x=s(t)$. This is an $O(\delta )$ boundary layer at $x=s(t)$, denoted region BI in figure 5. We look for an expansion solution in powers of $\eta _B$ in region BI, namely

(5.3a)$$\begin{gather} T_s=T_s^{0}(X,t)+\eta_BT_s^{1}(X,t)+O\left(\eta_B^{2}\right), \quad C=C^{0}(X,t)+\eta_BC^{1}(X,t)+O\left(\eta_B^{2}\right), \end{gather}$$

(5.3b)$$\begin{gather}T_g=T_g^{0}(X,t)+\eta_BT_g^{1}(X,t)+O\left(\eta_B^{2}\right), \quad f=f^{0}(X,t)+\eta_Bf^{1}(X,t)+O\left(\eta_B^{2}\right), \end{gather}$$

(5.3c)$$\begin{gather}s=s_0(t)+\eta_Bs_1(t)+O\left(\eta_B^{2}\right). \end{gather}$$

Substituting into (2.12), (2.12b) and (2.12d) for the gas variables decouple on this length scale, and so we obtain exactly the same system at leading order for just the solid variables $T_s^{0}$, $C^{0}$ and $s_0$ as in case A, namely (4.8). We will find below that $C=1+O(\eta _B)$ everywhere outside of region BI, and so the matching conditions as $X\rightarrow \infty$ are also identical to case A, given by (4.8d).

5.2. Region BII

Since the convective heat-transfer rate is non-zero, unlike case A, the equations for the gas variables do not decouple in (2.12) in the entire domain. Instead, outside the boundary layer BI, we must have heat transfer between the solid and gas phases, and expect the solid temperature scaling $T_s\approx \eta _B$. We therefore make the rescaling

(5.4)\begin{equation} T_s=\eta_B\left(1+\eta_B\hat{T}_s\right), \end{equation}

in the $O(1)$-length scale outer region, denoted region BII, and for clarity of notation we also write $C=\hat {C},$ $f=\hat {f}$ and $T_g=\hat {T}_g$. The chemical reaction rate here is therefore

(5.5)\begin{equation} \hat{Q}(\hat{C},\hat{T}_s)=\hat{C}\exp\left(A\left(1-\frac{1}{T_s}\right)\right)=\hat{C}e^{A}\delta\eta_B\left(e^{A\hat{T}_s}+O(\eta_B)\right), \end{equation}

and thus the heat consumed by the chemical reaction at this temperature $\eta _B$ is the same order ($\delta \eta _B$) as the heat transferred from the gas to the solid. We expand all dependent variables in powers of $\eta _B\ll 1$, setting

(5.6)\begin{equation} \left. \begin{aligned} \hat{T}_s & =\hat{T}_s^{0}(x,t)+O(\eta_B), \quad \hat{C}=\hat{C}^{0}(x,t)+O(\eta_B),\\ \hat{T}_g & =\hat{T}_g^{0}(x,t)+O(\eta_B), \quad \hat{f}=\hat{f}^{0}(x,t)+O(\eta_B). \end{aligned} \right\} \end{equation}

Since the reaction rate is small, at leading order we find that $\hat {f}^{0}$ and $\hat {T}_g^{0}$ are uniform in $x$. Matching at $x=s(t)$ with region BI, using the same notation as in (4.10), we find

(5.7a,b)\begin{equation} \hat{f}^{0}=F_s, \quad \hat{T}_g^{0}=G_s. \end{equation}

From (2.12a) we find that $\hat {C}^{0}$ is conserved along characteristics travelling right-to-left through region BII. As we will see below, $C$ is uniform to leading order over all regions to the right of this, and so by using the boundary condition at $x=1$, and matching, we will find that $\hat {C}^{0}=1$ is constant throughout region BII. The solid heat equation (2.12c), at leading order (order $\eta _B$), reduces to

(5.8)\begin{equation} M\hat{T}_g^{0}=\gamma \hat{C}^{0}\exp\left(A(1+\hat{T}_s^{0})\right), \end{equation}

so that $\hat {T}_s^{0}$ is also uniform, and given by

(5.9)\begin{equation} \hat{T}_s^{0}=\frac{1}{A}\log\left(\frac{MG_s}{\gamma}\right)-1. \end{equation}

This completes the leading-order analysis in region BII. Physically, in this region all the heat transferred from the gas to the solid is used for the chemical reaction, so that the solid temperature remains uniform. Region BII is a reaction region (like region BI), although here the heat for the chemical reaction is provided by the heat transfer between phases, rather than the radiation onto the free boundary.

5.3. Region BIV

For this size of the convective heat-transfer coefficient $\mu =O(\eta _B)$, there is also a second possible dominant balance on an $O(1)$ length scale, namely the balance of advection (term 1 in (3.1)) with the heat transfer from the gas (term 3 in (3.1)), when the chemical reaction rate is small $Q\ll \delta \eta _B$. We see that this balance occurs at the right-hand side of the domain (denoted region BIV) since, near the boundary $x=1$, $T_s\approx T^{in}$, and we note that $\eta _B>T^{in}$, by the assumption (2.14). In region BIV, the solid material is therefore heated up to the critical reacting temperature $\eta _B$ by the heat transfer from the gas. To study region BIV we make the change of variables

(5.10)\begin{equation} T_s=\eta_B\bar{T}_s, \end{equation}

for $\bar {T}_s<1$, and we write the boundary value $T^{in}=\eta _B\bar {T}^{in}$, where by (2.14) $\bar {T}^{in}<1$. For clarity of notation we also write $C=\bar {C}$, $f=\bar {f}$ and $T_g=\bar {T}_g$ in this region BIV. Then the chemical reaction rate in BIV is

(5.11)\begin{equation} Q=\bar{C}\exp\left(A\left(1-\frac{1}{\eta_B\bar{T}_s}\right)\right)=\bar{C}e^{A}(\delta\eta_B)^{1/\bar{T}_s}\ll\delta\eta_B. \end{equation}

Again we expand in powers of $\eta _B$, setting

(5.12)\begin{equation} \left. \begin{aligned} \bar{T}_s & =\bar{T}_s^{0}(x,t)+O(\eta_B), \quad \bar{C}=\bar{C}^{0}(x,t)+O(\eta_B), \\ \bar{T}_g & =\bar{T}_g^{0}(x,t)+O(\eta_B), \quad \bar{f}=\bar{f}^{0}(x,t)+O(\eta_B). \end{aligned} \right\} \end{equation}

At leading order in (2.12a), (2.12b) and (2.12d), we see that $\bar {C}^{0}, \bar {f}^{0}$ and $\bar {T}_g^{0}$ are all uniform through region BIV. From (2.12c) we find that $\bar {T}_s^{0}$ satisfies

(5.13)\begin{equation} \bar{C}^{0}\left(\frac{\partial{\bar{T}^{0}_s}}{\partial{t}}-\frac{\partial{\bar{T}^{0}_s}}{\partial{x}}\right)=M\bar{T}_g^{0}. \end{equation}

By the method of characteristics, applying the boundary conditions at $x=1$, we find

(5.14)\begin{equation} \bar{T}_s^{0}=\bar{T}^{in}+\int_{\xi=0}^{1-x} \frac{M\bar{T}_g^{0}(t+x-1+\xi)}{\bar{C}^{0}(t+x-1+\xi)}\,\mathrm{d}\xi.\end{equation}

The integrand in (5.14) is positive, and so for any fixed $t$, $\bar {T}_s^{0}$ increases monotonically as $x$ decreases from $1$ into the domain. Thus in region BIV no chemical reaction occurs since the solid material is below the critical temperature $\eta _B$, and the transfer of heat from the gas to the solid has the effect of heating up the solid material.

5.4. Region BIII

The solution in region BIV is only valid so long as $\bar {T}_s^{0}<1$. As $\bar {T}_s^{0}$ approaches 1 (so that $T_s$ approaches the critical temperature $\eta _B$), there must be a transition region, denoted region BIII, over which the chemical reaction rate (term 5 in (3.1)) enters the dominant balance. Physically, region BIII is the transition region where the heat transferred from gas to solid ceases to be used for the chemical reaction, and is instead used to heat the solid up from its input temperature $T^{in}$. To analyse this region we use the change of variables

(5.15)\begin{equation} T_s=\eta_B\left(1+\eta_B\tilde{T}_s\right). \end{equation}

The other variables are denoted $C=\tilde {C},\ f=\tilde {f}$ and $T_g=\tilde {T}_g$ in region BIII for clarity. Since we require all of the heat loss to the chemical reaction (term 5 in (3.1)), heat transfer from gas to solid (term 3), and advection of heat within the solid (term 1) to balance, we require an $O(\eta _B)$ length scale, and so make the change of variables $x=v(t)+\eta _By$, where $y=O(1)$, and $v(t)$, to be determined, is the position of the transition layer BIII. As in regions BII and BIV, we see from (2.12a), (2.12b) and (2.12d) that to leading order in $\eta _B$ all of $\tilde {C},\tilde {f}$ and $\tilde {T}_g$ are uniform through region BIII. Matching to both left (region BII) and right (region BIV), we see that throughout regions BII–BIV

(5.16a–c)\begin{equation} \hat{C}^{0}=\bar{C}^{0}=1, \quad \hat{f}^{0}=\bar{f}^{0}=F_s(t), \quad \hat{T}_g^{0}=\bar{T}_g^{0}=G_s(t). \end{equation}

With these changes of variable, and expanding in powers of $\eta _B$, at leading order in (2.12c) we obtain the ordinary differential equation (ODE)

(5.17)\begin{equation} -(\dot{v}+1)\frac{\partial{\tilde{T}_s^{0}}}{\partial{y}}=MG_s-\gamma \exp \left(A\left(1+\tilde{T}_s\right)\right), \end{equation}

for $\tilde {T}_s^{0}$, which has solution

(5.18)\begin{equation} \tilde{T}_s^{0}={-}1-\frac{1}{A}\log\left(\frac{\gamma}{MG_s}\left[1+B\exp\left(\frac{AMG_s}{\dot{v}+1}y\right)\right]\right). \end{equation}

The position, $v(t)$, of the transition layer, and the constant of integration, $B(t)$, must be determined by matching with the neighbouring regions. We note that $v$ is fixed by the leading-order matching

(5.19)\begin{equation} 1=\bar{T}_s^{0}\Big|_{x=v}=\bar{T}^{in}+M\int_{\xi=0}^{1-v} G_s(t+v-1+\xi)\,\mathrm{d}\xi. \end{equation}

Since the integrand $G_s>0$ is positive, the integral is decreasing with $v$, and there is therefore a unique solution, $v(t)<1$, of (5.19), for every $t$. The value of $B$ in (5.18) may be fixed by matching at higher order, which requires the $O(\eta _B)$ correction in all regions of the domain BI-BIV.

We note that the position $v(t)$ given by (5.19) requires no knowledge of the form $\tilde {T}_s^{0}$ of the solution in region BIII. The function of region BIII is therefore to simply smooth the transition from uniform $T_s$ in region BII, to $\bar {T}_s^{0}$ given by (5.14) in region BIV.

We also notice that the position, $v(t)$, of the transition layer lies within the domain so long as $v(t)>s_0(t)$. From (5.19) for $v(t)$, we see that if any of $M$, $G_s$ or $1-s_0$ are too small, then $v(t)< s_0(t)$. In this case there is no region BII or BIII: region BIV fills the entire domain $[s_0(t),1]$, except for the boundary layer at $x=s_0$. We see that our solution approaches the solution for case A in § 4 as $M$ becomes small, since in this case $v(t)< s_0(t)$, and so everywhere in $[s_0(t),1]$ we must have

(5.20a–d)\begin{equation} C\sim 1, \quad f\sim F_s, \quad T_g\sim G_s,\quad T_s=\eta_B\bar{T}_s^{0}=\eta_B\left(\bar{T}^{in}+O(M)\right)\sim T^{in}. \end{equation}

5.5. Composite solution

Since $C$, $f$ and $T_g$ are uniform at leading order outside region BI, the leading-order composite solutions are simply

(5.21a–c)\begin{equation} C^{B}=C^{0}\left(\frac{x-s_0(t)}{\delta}\right), \quad f^{B}=f^{0}\left(\frac{x-s_0(t)}{\delta}\right), \quad T_g^{B}=T_g^{0}\left(\frac{x-s_0(t)}{\delta}\right). \end{equation}

As in case A, we note that $T_s=O(1)$ in region BI, but is $O(\eta _B)$ in all other regions. We construct a leading-order composite solution by assuming that the $O(\eta _B)$ correction in region BI is simply uniform, and equal to $\eta _B$, so that the composite $T_s$ solution is accurate with $O(\eta _B)$ error in region BI, and with $O(\eta _B^{2})$ error in regions BII–BIV. The solution is defined piecewise on either side of the transition layer BIII, centred at the point $x=v(t)$, where the solution is continuous by construction. We therefore use the composite solution

(5.22)\begin{equation} T_s^{B}=\begin{cases} T_s^{0}\left(\dfrac{x-s_0(t)}{\delta}\right)+\eta_B, & \text{if }x< v(t),\\ \eta_B\bar{T}_s^{0}(x,t), & \text{if }x\geq v(t), \end{cases} \end{equation}

where $\bar {T}_s^{0}$ is given by (5.14).

Thus in case B, while $\mu =O(\eta _B)$ is non-zero and so there is some heat transfer between phases, only an $O(\eta _B)$ amount of chemical reaction occurs outside region BI, and this only if $v(t)>s_0(t)$. The primary role of the heat transferred from the gas to the solid is to heat up the solid material. While this does not affect the leading-order problem in region BI, which is the same as that in AI in case A, it reduces its accuracy: there is now an $O(\eta _B)$ ‘switchback’ correction term in region BI, greater than the $O(T^{in})$ correction in case A.

6.1. Region CI

As in the previous cases, we find an $O(\delta )$ reaction region (region CI) at $x=s(t)$, where the heat for the reaction is provided by the radiation onto the free boundary. As before, we make the change of variables $x=s(t)+\delta X$, and expand the dependent variables in powers of $\eta _C$, by

(6.2a)$$\begin{gather} T_s=T_s^{0}(X,t)+\eta_CT_s^{1}(X,t)+O\left(\eta_C^{2}\right),\quad C=C^{0}(X,t)+\eta_CC^{1}(X,t)+O\left(\eta_C^{2}\right), \end{gather}$$

(6.2b)$$\begin{gather}T_g=T_g^{0}(X,t)+\eta_CT_g^{1}(X,t)+O\left(\eta_C^{2}\right), \quad f=f^{0}(X,t)+\eta_Cf^{1}(X,t)+O\left(\eta_C^{2}\right), \end{gather}$$

(6.2c)$$\begin{gather}s=s_0(t)+\eta_Cs_1(t)+O\left(\eta_C^{2}\right). \end{gather}$$

We obtain the same leading-order system of (4.8a)–(4.8c) in region CI as in cases A and B. As $X\rightarrow \infty$, we still require $T_s^{0}\rightarrow 0$ in order for the problem to be correctly specified. However, since the chemical reaction rate is $O(1)$ in the outer region, we will no longer have uniform concentration at leading order in the outer region. We must therefore impose the far-field conditions

(6.3)\begin{equation} T_s^{0}\rightarrow 0,\quad C^{0}\rightarrow C^{\infty}(t) \quad \text{ as }X\rightarrow\infty,\end{equation}

where $C^{\infty }(t)$ is yet to be determined by matching with the outer region. Thus we can no longer solve the region CI problem independently of the rest of the domain: it is now intrinsically coupled with the outer problem, since $C$ varies at leading order over the outer region. The leading-order gas variables within region CI are given by (4.9), but the far-field values as $X\rightarrow \infty$ are now

(6.4a)$$\begin{gather} f^{0}\rightarrow F_{s}(t):=f_*+(1+\dot{s}_0)(C^{\infty}-C_*), \end{gather}$$

(6.4b)$$\begin{gather}T_g^{0}\rightarrow G_{s}(t):=\frac{\dfrac{f_*\rho}{s_0}+K^{{-}1}\int_{\bar{X}=0}^{\infty} T_s^{0}C^{0}\exp\left(A\left(1-\dfrac{1}{T^{0}_s}\right)\right)\,\mathrm{d}\bar{X}}{f_*+(1+\dot{s}_0)(C^{\infty}-C_*)}. \end{gather}$$

6.2. Region CII

Outside region CI, over an $O(1)$ length scale, the heat transferred from the gas to the solid is sufficient that the chemical reaction is a dominant effect, and the solid temperature remains close to uniform. We refer to this region as region CII. Like region BII in case B, region CII is a reaction region within which the energy for the reaction is provided by the heat transferred between phases. Unlike region BII, however, the heat transfer is now fast enough that the gas temperature decreases through this region, rather than remaining uniform. In region CII, we make the rescaling

(6.5)\begin{equation} T_s=\eta_C\left(1+\eta_C\hat{T}_s\right), \end{equation}

and we also use the notation $C=\hat {C}$, $f=\hat {f}$, and $T_g=\hat {T}_g$ for clarity. We expand the dependent variables in powers of $\eta _C$, setting

(6.6)\begin{equation} \left. \begin{aligned} \hat{T}_s & =\hat{T}_s^{0}(x,t)+O(\eta_C), \quad \hat{C}=\hat{C}^{0}(x,t)+O(\eta_C), \\ \hat{T}_g & =\hat{T}_g^{0}(x,t)+O(\eta_C), \quad \hat{f}=\hat{f}^{0}(x,t)+O(\eta_C). \end{aligned} \right\} \end{equation}

At leading order, the system (2.12a)–(2.12d) becomes

(6.7a)$$\begin{gather} \frac{\partial{\hat{C}^{0}}}{\partial{t}}-\frac{\partial{\hat{C}^{0}}}{\partial{x}}={-}\hat{C}^{0}\exp\left(A\left(1+\hat{T}_s^{0}\right)\right), \end{gather}$$

(6.7b)$$\begin{gather}\frac{\partial{\hat{f}^{0}}}{\partial{x}}=\hat{C}^{0}\exp\left(A\left(1+\hat{T}_s^{0}\right)\right), \end{gather}$$

(6.7c)$$\begin{gather}\mu\hat{T}_g^{0}=\gamma\hat{C}^{0}\exp\left(A\left(1+\hat{T}_s^{0}\right)\right), \end{gather}$$

(6.7d)$$\begin{gather}K\frac{\partial}{\partial x}\left(\hat{f}^{0}\hat{T}^{0}_g\right)={-}\mu\hat{T}_g^{0}. \end{gather}$$

Combining (6.7b)–(6.7d) we obtain the reduced system

(6.8a,b)\begin{equation} \frac{\partial{\hat{f}^{0}}}{\partial{x}}=\frac{\mu}{\gamma}\hat{T}_g^{0}, \quad K\frac{\partial}{\partial x}\left(\hat{f}^{0}\hat{T}^{0}_g\right)={-}\mu\hat{T}_g^{0}, \end{equation}

for  $\hat {f}^{0}$ and $\hat {T}_g^{0}$. Matching with region CI, we must have $\hat {f}^{0}=F_{s}(t)$ and $\hat {T}_g^{0}=G_{s}(t)$ at $x=s_0$. Combining (6.8a,b), and using the boundary conditions at $x=s_0$, we find the first integral

(6.9)\begin{equation} \hat{f}^{0}\left(\gamma+K\hat{T}_g^{0}\right)=F_s(\gamma+KG_s).\end{equation}

Using (6.9) to eliminate $\hat {f}^{0}$ from (6.8a), we obtain a separable equation for $\hat {T}_g^{0}$, namely

(6.10)\begin{equation} \frac{1}{\hat{T}_g^{0}(\gamma+K\hat{T}_g^{0})^{2}}\frac{\partial{\hat{T}_g^{0}}}{\partial{x}}={-}\frac{\mu}{\gamma KF_s(\gamma+KG_s)}, \end{equation}

with solution

(6.11a)\begin{equation} \hat{T}_g^{0}(x,t)={-}\frac{\gamma W(a_{C})}{K(1+W(a_{C}))},\end{equation}

where $W$ is the principal branch of the Lambert W function, and $a_{C}=a_{C}(x,t)$ is given by

(6.11b)\begin{equation} a_{C}(x,t)={-}\frac{KG_{s}}{\gamma+KG_{s}}\exp\left(-\frac{1}{\gamma+KG_{s}}\left(\frac{\gamma\mu}{KF_{s}}(x-s_0)+KG_{s}\right)\right). \end{equation}

From (6.9), we therefore find

(6.11c)\begin{equation} \hat{f}^{0}=\frac{F_{s}(\gamma+KG_{s})}{\gamma+K\hat{T}_g^{0}}. \end{equation}

Given $\hat {T}_g^{0}$ from (6.11), the leading-order concentration $\hat {C}^{0}$ satisfies (6.7a), or

(6.12)\begin{equation} \frac{\partial{\hat{C}^{0}}}{\partial{t}}-\frac{\partial{\hat{C}^{0}}}{\partial{x}}={-}\frac{\mu}{\gamma}\hat{T}_g^{0}. \end{equation}

Matching with region CI at $x=s_0(t)$, we require $\hat {C}^{0}=C^{\infty }(t)$ at $x=s_0(t)$, and we will see below that the correct matching condition for $\hat {C}^{0}$ at the right side of region CII is $\hat {C}^{0}=1$ at $x=1$. By the method of characteristics, (6.12) with the boundary condition at $x=1$ has solution

(6.13)\begin{equation} \hat{C}^{0}(x,t)=1-\frac{\mu}{\gamma}\int_{\xi=0}^{1-x}\hat{T}_g^{0}(1-\xi,\xi+t+x-1)\,\mathrm{d}\xi,\end{equation}

for $\hat {T}_g^{0}(x,t)$ given by (6.11a). Evaluating at $x=s_0(t)$, the matching condition for $\hat {C}^{0}$ gives

(6.14)\begin{equation} C^{\infty}(t)=1-\frac{\mu}{\gamma}\int_{\xi=0}^{1-s_0(t)}\hat{T}_g^{0}(1-\xi,\xi+t+s_0(t)-1)\,\mathrm{d}\xi.\end{equation}

Since $\hat {T}_g^{0}$ depends on $s_0$, $F_s$ and $G_s$, we see that $C^{\infty }$ depends on these as well. Physically, the heat and flux of gas leaving region CI determine the rate of solid material consumption in region CII, which affects the position of the free boundary $s_0$ (found as part of the solution of the region CI problem), and therefore the gas temperature. The processes in regions CI and CII are thus fully coupled. However, since we have solved the region CII problem to find $C^{\infty }(t)$, at leading order we are left with the region CI problem (4.8a)–(4.8c) with far-field condition (6.3), with $C^{\infty }(t)$ given by (6.14).

Concluding the analysis in region CII, we find from (6.7c) that

(6.15)\begin{equation} \hat{T}_s^{0}=\frac{1}{A}\log\left(\frac{\mu\hat{T}_g^{0}}{\gamma \hat{C}^{0}}\right)-1, \end{equation}

given the forms (6.11a) of $\hat {T}_g^{0}$ and (6.13) of $\hat {C}^{0}$.

6.3. Region CIV

In region CII, we have $T_s=\eta _C+O(\eta _C^{2})$, and so we cannot impose the boundary condition $T_s=T^{in}$ at $x=1$ (as discussed for case B, the condition (2.14) ensures $T^{in}<\eta _C$). There must therefore be a second boundary layer at $x=1$, within which the solid temperature is heated from its input temperature $T^{in}$ up to the critical temperature $\eta _C$. This heating region is analogous to region BIV in case B, and we denote it region CIV. In region CIV we change variables, setting

(6.16)\begin{equation} T_s=\eta_C\bar{T}_s, \end{equation}

with $\bar {T}_s<1$, so that

(6.17)\begin{equation} \bar{Q}:=Q(C,\eta_C\bar{T}_s)=O\left(\delta^{1/\bar{T}_s}\right)\ll\delta, \end{equation}

and thus we cannot expect a dominant balance in the solid heat equation (2.12c) between the heating from the gas and heat lost to chemical reaction. Instead, we find a balance between the heating from the gas and the advection term, which fixes the length scale of region CIV to be $O(\eta _C)$. Within region CIV we therefore make the change of variables $x=1-\eta _C Y$, where $Y=O(1)$ is positive. For clarity of notation, we set $C=\bar {C}$, $f=\bar {f}$ and $T_g=\bar {T}_g$ in region CIV, and we also rewrite the input temperature $T^{in}=\eta _C\bar {T}^{in}$, so that the boundary condition becomes

(6.18)\begin{equation} \bar{T}_s=\bar{T}^{in} \quad \text{ at }Y=0. \end{equation}

Again, we expand in powers of $\eta _C$, setting

(6.19)\begin{equation} \left. \begin{aligned} \bar{T}_s & =\bar{T}_s^{0}(Y,t)+O(\eta_C), \quad \bar{C}=\bar{C}^{0}(Y,t)+O(\eta_C), \\ \bar{T}_g & =\bar{T}_g^{0}(Y,t)+O(\eta_C), \quad \bar{f}=\bar{f}^{0}(Y,t)+O(\eta_C). \end{aligned} \right\} \end{equation}

Due to the short length scale, at leading order in (2.12a), (2.12b) and (2.12d), we see that all of $\bar {C}^{0}$, $\bar {f}^{0}$ and $\bar {T}_g^{0}$ are uniform across region CIV. The boundary condition at $Y=0$ (or $x=1$) fixes $\bar {C}^{0}=1$, while $\bar {f}^{0}$ and $\bar {T}_g^{0}$ must be fixed by matching to the left. At leading order in (2.12c), we find that the solid temperature $\bar {T}_s^{0}$ satisfies

(6.20)\begin{equation} \frac{\partial{\bar{T}_s^{0}}}{\partial{Y}}=\mu\bar{T}_g^{0}, \end{equation}

and thus, using the boundary condition (6.18), we find that $\bar {T}_s^{0}$ is linear in $Y$

(6.21)\begin{equation} \bar{T}_s^{0}=\bar{T}^{in}+\mu\bar{T}_g^{0}Y. \end{equation}

6.4. Region CIII

Exactly as in case B, the solution (6.21) in region CIV is valid only while $\bar {T}_s^{0}<1$, and we cannot directly match the linear profile (6.21) with the constant $T_s\sim \eta _C$ in region CII. Thus, as in case B, we must introduce a transition layer, region CIII, between regions CII and CIV, in which the heat transferred from the gas is no longer used in the chemical reaction, but is instead used to heat the solid material from its input temperature $T^{in}$. In this transition layer, we must have $T_s\approx \eta _C$, and we expect a balance between all three of heat advection (term 1), heat transfer from the gas (term 3) and heat loss to the chemical reaction (term 5) in (3.1). The length scale of such a region must therefore be $O(\eta _C^{2})$. We may show, as in case B, that all of $f$, $T_g$ and $C$ are uniform across region CIII, so that matching across CII–CIV, we find

(6.22a–c)\begin{equation} \hat{C}^{0}|_{x=1} =\bar{C}^{0}=1, \quad F_1(t):=\hat{f}^{0}|_{x=1} =\bar{f}^{0}, \quad G_1(t):=\hat{T}_g^{0}|_{x=1} =\bar{T}_g^{0}. \end{equation}

(This fixes the boundary condition at $x=1$ for $\hat {C}^{0}$, in region CII.) As in case B, we find that the ODE for the solid temperature variation simply smooths the transition from region CII to region CIV. The changes of variable for this analysis are summarised in table 3. The position $w(t)$ of the transition layer is found by matching $T_s$ across regions CIII–CIV, which, at leading order, requires $1=\bar {T}^{in}+\mu G_1 w,$ so that

(6.23)\begin{equation} w(t)=\frac{1}{\mu G_1(t)}\left(1-\bar{T}^{in}\right).\end{equation}

6.5. Composite solution

As in previous cases, we may construct a leading-order composite solution. Since $C$, $f$ and $T_g$ are $O(1)$ in all regions, and uniform in regions CIII and CIV, the composite solutions are

(6.24)$$\begin{gather} C^{C}=C^{0}\left(\frac{x-s_0(t)}{\delta}\right)+\hat{C}^{0}(x,t)-C^{\infty}(t), \end{gather}$$

(6.25)$$\begin{gather}f^{C}=f^{0}\left(\frac{x-s_0(t)}{\delta}\right)+\hat{f}^{0}(x,t)-F_s(t), \end{gather}$$

(6.26)$$\begin{gather}T_g^{C}=T_g^{0}\left(\frac{x-s_0(t)}{\delta}\right)+\hat{T}_g^{0}(x,t)-G_s(t). \end{gather}$$

The solid temperature is $O(1)$ in region CI but $O(\eta _C)$ elsewhere. We therefore construct a composite solution by setting the $O(\eta _C)$ correction in region CI to be simply $\eta _C$, so that the composite solution is accurate at leading order in all regions. As in case B, the composite solution for $T_s$ is defined piecewise about the centre, $1-\eta _Cw(t)$, of the transition layer CIII. Thus we use the composite solution

(6.27)\begin{equation} T_s^{C}=\begin{cases} T_s^{0}\left(\dfrac{x-s_0(t)}{\delta}\right)+\eta_C, & \text{if }x<1-\eta_C w(t),\\ \eta_C\bar{T}_s^{0}\left(\dfrac{1-x}{\eta_C}\right), & \text{if }x\geq 1-\eta_C w(t). \end{cases} \end{equation}

We have seen that case C is very similar to case B. In each, the solid is heated up from its input temperature $T^{in}$ by heat transferred from the gas in region IV. Once it has reached the critical temperature $\eta$ the chemical reaction then starts to take place, and the heat transferred from the gas to the solid in this region (II) is used primarily for this reaction, with no change in the solid temperature at leading order until it reaches the hot boundary layer (I) at $x=s$. The differences in case C are that, since the convective heat transfer from the gas is now greater, the solid material heats up much faster (over the $O(\eta _C)$ length scale of region CIV, rather than the $O(1)$ length scale of BIV), and that the reaction rate in region CII is now high enough to affect all of the solid concentration $C$, gas mass flux $f$ and gas temperature $T_g$ at leading order. As we increase $\mu$ still further, the gas temperature must lose so much heat to the solid that it ceases to be $O(1)$. This is the situation considered in case D below.

7.1. Behaviour common to cases Da and Db

There are several regions of the domain near to the free boundary $x=s$ where the behaviour is the same for both cases Da and Db.

7.1.1. Region DI

As usual we find a reaction boundary layer, region DI, of width $O(\delta )$ at $x=s$, within which the heat for the reaction is provided by the radiation onto the free boundary. Expanding the dependent variables in powers of $\eta _D$, by setting

(7.1a)$$\begin{gather} T_s=T_s^{0}(X)+\eta_DT_s^{1}(X)+O\left(\eta_D^{2}\right), \quad C=C^{0}(X)+\eta_DC^{1}(X)+O\left(\eta_D^{2}\right), \end{gather}$$

(7.1b)$$\begin{gather}T_g=T_g^{0}(X)+\eta_DT_g^{1}(X)+O\left(\eta_D^{2}\right), \quad f=f^{0}(X)+\eta_Df^{1}(X)+O\left(\eta_D^{2}\right), \end{gather}$$

(7.1c)$$\begin{gather}s=s_0+\eta_Ds_1+O\left(\eta_D^{2}\right), \end{gather}$$

we obtain the same equations (4.8a)–(4.8c) at leading order since $\mu \ll \delta ^{-1}$ (although with $\dot {s}_0=0$, since we restrict ourselves to steady-state solutions). The far-field solid concentration must satisfy (6.3) as in case C, where $C^{\infty }$ is yet to be determined by matching with the outer problem.

7.1.2. Region DII

Outside region DI, we require $T_s\sim \eta _D$, in order to balance the terms 3 and 5 in the solid heat equation (3.1). However, since then $Q=O(\mu )\gg 1$, this region, which we denote region DII, must have a shorter length scale than $O(1)$ in order to obtain a balance in the gas temperature equation. In this region the heat for the reaction is provided by the heat transfer between phases. Region DII is therefore a second, wider boundary layer, width $O(\mu ^{-1})$, about $x=s_0$, and so we make the change of variables $x=s_0+y/\mu$. We look for a solution

(7.2)\begin{equation} T_s=\eta_D\left(1+\eta_D\hat{T}_s\right), \end{equation}

and we also use the notation $C=\hat {C}$, $f=\hat {f}$, and $T_g=\hat {T}_g$ in region DII for clarity. We expand in powers of $\eta _D\ll 1$, setting

(7.3)\begin{equation} \left. \begin{aligned} \hat{T}_s & =\hat{T}_s^{0}(y)+O(\eta_D), \quad \hat{C}=\hat{C}^{0}(y)+O(\eta_D),\\ \hat{T}_g & =\hat{T}_g^{0}(y)+O(\eta_D), \quad \hat{f}=\hat{f}^{0}(y)+O(\eta_D). \end{aligned} \right\} \end{equation}

Substituting into (2.12a)–(2.12d), we find that at leading order

(7.4a)$$\begin{gather} -\frac{\mathrm{d}{\hat{C}^{0}}}{\mathrm{d}{y}}={-}\hat{C}^{0}\exp\left(A\left(1+\hat{T}_s^{0}\right)\right), \quad \hat{T}_g^{0}=\gamma\hat{C}^{0}\exp\left(A\left(1+\hat{T}_s^{0}\right)\right), \end{gather}$$

(7.4b)$$\begin{gather}\frac{\mathrm{d}{\hat{f}^{0}}}{\mathrm{d}{y}}=\hat{C}^{0}\exp\left(A\left(1+\hat{T}_s^{0}\right)\right), \quad K\frac{\mathrm{d}}{\mathrm{d} y}(\hat{f}^{0}\hat{T}_g^{0})={-}\hat{T}_g^{0}. \end{gather}$$

Precisely as in (6.11) for region CII of case C, we find the solutions

(7.5a,b)\begin{equation} \hat{T}_g^{0}={-}\frac{\gamma W(a_{D})}{K(1+W(a_{D}))}, \quad \hat{f}^{0}=\frac{F_s(\gamma+KG_s)}{\gamma +K\hat{T}_g^{0}}, \end{equation}

of (7.4b), where

(7.6)\begin{equation} a_{D}={-}\frac{KG_s}{\gamma+KG_s}\exp\left(-\frac{1}{\gamma+KG_s}\left(\frac{\gamma }{KF_s}y+KG_s\right)\right), \end{equation}

and the matching constants $F_s$ and $G_s$ are as defined in (6.4) (with $\dot {s}_0=0$, since we are in steady state). The leading-order solid concentration therefore satisfies

(7.7)\begin{equation} -\frac{\mathrm{d}{\hat{C}^{0}}}{\mathrm{d}{y}}=\frac{W(a_{D})}{K(1+W(a_{D}))},\end{equation}

along with $\hat {C}^{0}=C^{\infty }$ at $y=0$, and a matching condition to the right, as $y$ becomes large. We will find below that $C$ is uniform to leading order in all regions to the right of region DII, in both cases Da and Db. Thus the matching condition is $\hat {C}^{0}\rightarrow 1$ as $y\rightarrow \infty.$ From (7.5a) we see that $\hat {T}_g^{0}\rightarrow 0$ as $y\rightarrow \infty$, and so the solution of (7.7) is

(7.8)\begin{equation} \hat{C}^{0}=1-\frac{KF_s(\gamma+KG_s)\hat{T}_g^{0}}{\gamma(\gamma+K\hat{T}_g^{0})}. \end{equation}

Applying the boundary conditions at $y=0$, where $\hat {T}_g^{0}=G_s$, we therefore fix the value of $C^{\infty }$ to be

(7.9)\begin{equation} C^{\infty}=1-\frac{KF_sG_s}{\gamma}. \end{equation}

Having found $\hat {T}_g^{0}$ and $\hat {C}^{0}$, the solid temperature variation, given by (7.4a) (right), is then

(7.10)\begin{align} \hat{T}_s^{0}&=\frac{1}{A}\log\left(\frac{\hat{T}_g^{0}}{\gamma \hat{C}^{0}}\right)-1\nonumber\\ &={-}1-\frac{1}{A}\log\left(K\left(1+W(a_{D})\right)\left(\frac{1}{-W(a_{D})}-\frac{F_s(\gamma+KG_s)}{\gamma}\right)\right). \end{align}

We note that this region DII is exactly analogous to region CII in case C with $\mu =O(1)$, except that, since the convective heat-transfer coefficient $\mu$ is now larger, the decay length of $\hat {T}_g^{0}$ is shorter. We see that (7.4) hold until $\hat {T}_g^{0}$ becomes order $\eta _D$, at which point term 4 must enter the dominant balance of (3.1). Since for small $a_{D}$, $W(a_{D})\sim a_{D}$, we see that $\hat {T}_g^{0}\sim \eta _D$ when $-\gamma a_{D}/K=\eta _D$, which occurs at the point $y=y_*$, defined by $\gamma a_{D}(y_*)/K=\eta _D$. Rearranging, this gives

(7.11)\begin{equation} y_*=\frac{1}{\gamma}KF_s(\gamma+KG_s)\left(-\log(\eta_D)+\log\left(\frac{\gamma G_s}{\gamma +KG_s}\right)-\frac{KG_s}{\gamma+KG_s}\right).\end{equation}

We see that point $y_*=O(-\log (\eta _D))$ is large. As $y$ becomes larger than $O(1)$, we strictly leave region DII, and enter region DII’, an extension of DII, within which the length scale is $O(-\log (\eta _D)/\mu )$. The dominant physical processes in region DII’ (heat transfer from gas to solid, which provides energy for the chemical reaction) are the same as in region DII, except that $f$ is uniform at leading order, and thus $T_g$ is exponentially decaying.

In the original variables, the point $y_*$ is at $x_1:=s_0+y_*/\mu$. As mentioned at the start of this section, we notice that if $1\ll \mu \ll -\log (\eta _D)$ then $x_1-s_0\gg 1$ and so we revert to a structure like case C, since $\hat {T}_g$ remains $O(1)$ throughout the domain. Conversely, in the regime $-\log (\eta _D)\ll \mu$ studied in this section, we have $x_1-s_0\ll 1$ and so $\hat {T}_g^{0}$ becomes small within the domain. In the distinguished limit $\mu =O(-\log (\eta _D))$, we see that $x_1-s_0=O(1)$, and we may be in either of these cases, or instead have $x_1$ close to $1$, for which a number of additional subcases arise.

7.1.3. Region DIII

With $\mu \gg -\log (\eta _D)$, the point $x_1$ is an $O(-\log (\eta _D)/\mu )\ll 1$ distance from $s_0$ (and so is definitely contained within the domain). Around this point $x_1$, there must be an $O(\mu ^{-1})$ transition layer, denoted region DIII, in which term 4 enters the dominant balance in (3.1). In this region we use the change of spatial variables $x=x_1+z/\mu$. Since now $T_g=O(\eta _D)$, the balance of terms 3 and 5 requires that term 5, the chemical reaction rate, be $O(\delta \mu \eta _D)$. Heat transfer between the two phases is slower than in regions DII-DII’ since the difference between the two temperatures is now smaller, and so the chemical reaction driven by this temperature difference is slower in this region than in region DII and DII’. Thus, the solid temperature must be of the form

(7.12)\begin{equation} T_s=\eta_D\left(1+\frac{\eta_D\log(\eta_D)}{A}+\eta_D\mathcal{T}_s\right).\end{equation}

We also use the change of variables $T_g=\eta _D\mathcal {T}_g$ in region DIII, and set $f=\mathcal {F}$, $C=\mathcal {C}$ for clarity of notation, and we expand in powers of $\eta _D$, setting

(7.13)\begin{equation} \left. \begin{aligned} \mathcal{T}_s & =\mathcal{T}_s^{0}(z)+O(\eta_D), \quad \mathcal{C}=\mathcal{C}^{0}(z)+O(\eta_D), \\ \mathcal{T}_g & =\mathcal{T}_g^{0}(z)+O(\eta_D), \quad \mathcal{F}=\mathcal{F}^{0}(z)+O(\eta_D). \end{aligned} \right\} \end{equation}

At leading order in (2.12a) and (2.12b), we find that both $\mathcal {C}^{0}$ and $\mathcal {F}^{0}$ are uniform in region DIII. This is because the chemical reaction rate is now $Q=O(\delta \eta _D\mu )$, and so there is only an $O(\eta _D)$ mass consumption over region DIII. Indeed, we will find that in all regions to the right of region DIII (in either case Da or Db), $C$ and $f$ are uniform up to $O(\eta _D)$. Thus by matching to the right we will find that $\mathcal {C}^{0}=1$, and matching at leading order with region DII, we have

(7.14)\begin{equation} \mathcal{F}^{0}=F_*:=\frac{F_s(\gamma+KG_s)}{\gamma}=\hat{f}^{0}|_{y=y_*}+O(\eta_D). \end{equation}

At leading order, the heat equations (2.12c) and (2.12d) for $\mathcal {T}_s^{0}$ and $\mathcal {T}_g^{0}$ become

(7.15)$$\begin{gather} \mathcal{T}_g^{0}-\left(1+\eta_D\frac{\log(\eta_D)}{A}\right)=\gamma \exp\left(A\left(1+\mathcal{T}_s^{0}\right)\right), \end{gather}$$

(7.16)$$\begin{gather}KF_*\frac{\mathrm{d}{\mathcal{T}_g^{0}}}{\mathrm{d}{z}}={-}\mathcal{T}_g^{0}+\left(1+\eta_D\frac{\log(\eta_D)}{A}\right). \end{gather}$$

Here, we have included the $O(-\eta _D\log (\eta _D))$ terms, since these are larger than the next order $\eta _D$, and we can find an explicit solution taking both $O(1)$ and $O(-\eta _D\log (\eta _D))$ terms into account together. The ODE (7.16) for $\mathcal {T}_g^{0}$ has the general solution

(7.17)\begin{equation} \mathcal{T}_g^{0}=1+\eta_D\frac{\log(\eta_D)}{A}+\exp\left(-\frac{1}{KF_*}z\right), \end{equation}

where we have fixed the constant of integration by matching with region DII. Finally, from (7.15) we compute

(7.18)\begin{equation} \mathcal{T}_s^{0}={-}\frac{1}{A}\left(\log(\gamma)+\frac{1}{KF_*}z\right)-1, \end{equation}

completing the analysis in region DIII.

7.1.4. Regions DaV and DbVI

Regions DI–DIII as described above are common to both cases Da and Db; these are all reaction regions, with the chemical reaction playing a dominant role. We now pause in our analysis moving left to right through the domain, and instead investigate the behaviour near the $x=1$ boundary, where there is negligible chemical reaction. We will find different behaviour for the cases Da and Db, depending on whether the overall heat flux $Kf-1=KF_*-1+O(\eta _D)$ is negative or positive, respectively.

On the right side of the domain, the solid temperature must decrease from close to $\eta _D$ at the left of the domain to the input value $T^{in}=\eta _D\bar {T}^{in}$. (As in regions BIV and CIV, here the heat transferred from the gas to the solid is used to heat the solid material from its input temperature, $T^{in}$). We therefore rescale

(7.19a,b)\begin{equation} T_s=\eta_D\bar{T}_s \quad \text{and} \quad T_g=\eta_D\bar{T}_g, \end{equation}

noting that $T_g$ must be $O(\eta _D)$ in this region too because it is $O(\eta _D)$ in region DIII, and can only decrease through the rest of the domain. With $\bar {T}_s<1+\eta _D\log (\eta _D)/A$, we notice that the chemical reaction rate $\bar {Q}\ll \delta \mu \eta _D$ is small. In order to balance the heat advection with heat transfer between phases in the gas heat equation, we must scale the space variable to an $O(1/\mu )$ domain. Thus the region over which $T_s$ changes must be a boundary layer. We make the change of variables $x=x_c+Z/\mu$, where $Z=O(1)$, and $x_c$ (the position of this boundary layer) is an $O(1)$ point in the domain, to be determined. As well as (7.19a,b), we set $f=\bar {f}$ and $C=\bar {C}$ in this region for clarity of notation, but note that since $\bar {Q}\ll \delta \mu \eta _D$ is small, $\bar {C}=1+O(\eta _D)$ and $\bar {f}^{0}=F_*+O(\eta _D)$ are uniform over this region at leading order, fixed by matching with regions to the right and left, respectively. We expand the temperatures in powers of $\eta _D$, by setting

(7.20a,b)\begin{equation} \bar{T}_s=\bar{T}_s^{0}(Z)+O(\eta_D), \quad \bar{T}_g=\bar{T}_g^{0}(Z)+O(\eta_D). \end{equation}

At leading order in (2.12c) and (2.12d), we therefore find

(7.21a,b)\begin{equation} -\frac{\mathrm{d}{\bar{T}_s^{0}}}{\mathrm{d}{Z}}=\bar{T}_g^{0}-\bar{T}_s^{0}, \quad KF_*\frac{\mathrm{d}{\bar{T}_g^{0}}}{\mathrm{d}{Z}}={-}\bar{T}_g^{0}+\bar{T}_s^{0}, \end{equation}

which admit a general solution

(7.22a,b)\begin{equation} \bar{T}_s^{0}=b_1+KF_*b_2\exp\left(\frac{KF_*-1}{KF_*} Z\right), \quad \bar{T}_g^{0}=b_1+b_2\exp\left(\frac{KF_*-1}{KF_*}Z\right), \end{equation}

where $b_1$ and $b_2$ are constants to be determined. We notice that the position of the boundary layer, $x_c$, must depend on the sign of $KF_*-1$, and we consider the cases Da ($KF_*<1$) and Db ($KF_*>1$) separately below. In case Da with $KF_*<1$, we find that this boundary layer must be at the left side of the domain (region DaV), whereas in case Db with $KF_*>1$, we must take $x_c=1$ so that the boundary layer is region DbVI, on the right side of the domain. From the non-dimensionalisation, we see that

(7.23)\begin{equation} KF_*<1 \quad \text{if}\quad c_{p,g}f< c_{p,s}u_sC_s, \end{equation}

and vice versa. Thus, since both phases have the same temperature at leading order on an $O(1)$ length scale, we are comparing the heat flux in the gas (left-to-right) with that in the solid (right-to-left).

7.2. Case Da: $KF_*<1$

If $KF_*<1$, then the heat flux in the solid is greater than that in the gas, and so the overall heat transport is right-to-left.

7.2.1. Region DaV

From the form of the solution (7.22a,b), we require that $x_c$ is on the left side of the domain, so that the solution remains bounded as $Z\rightarrow \infty$ and we leave the boundary layer. This boundary layer is denoted region DaV in figure 7.

As $Z\rightarrow \infty$, we have $\bar {T}_s^{0},\bar {T}_g^{0}\rightarrow b_1$, so in this case we must take $b_1=\bar {T}^{in},$ in order to satisfy the boundary condition at $x=1$. The value of $b_2$ is fixed by matching at $Z=0$, or $x=x_c$. We notice that the difference in solutions at $Z=0$ is

(7.24)\begin{equation} (\bar{T}_g^{0}-\bar{T}_s^{0})\Big|_{Z=0}=(1-KF_*)b_2=O(1), \end{equation}

thus $T_g-T_s=O(\eta _D)$ at $Z=0$. We must therefore have $x_c$ within the $O(1/\mu )$ transition layer DIII about $x_1$, in which $T_g-T_s=O(\eta _D)$. Since $T_s=\eta _D(1+\eta _D\log (\eta _D)/A)+O(\eta _D^{2})$ (by (7.12)) throughout region DIII, we fix $b_2$ to be

(7.25)\begin{equation} b_2=\frac{1}{KF_*}\left(1+\eta_D\frac{\log(\eta_D)}{A}-\bar{T}^{in}\right). \end{equation}

Thus the solutions (7.22a,b) in this boundary-layer region DaV (to the right of $x=x_c$) are

(7.26a)$$\begin{gather} \bar{T}_s^{0}=\bar{T}^{in}+\left(1+\eta_D\frac{\log(\eta_D)}{A} -\bar{T}^{in}\right)\exp\left(-\frac{1-KF_*}{KF_*} Z\right), \end{gather}$$

(7.26b)$$\begin{gather}\bar{T}_g^{0}=\bar{T}^{in}+\frac{1}{KF_*}\left(1+\eta_D\frac{\log(\eta_D)}{A} -\bar{T}^{in}\right)\exp\left(-\frac{1-KF_*}{KF_*} Z\right). \end{gather}$$

At $x=x_c$, or $Z=0$, we have $\bar {T}_g^{0}=(1+\eta _D\log (\eta _D)/A)/(KF_*)$. Imposing continuity of $T_g$ at $x=x_c$ determines the position of $x_c$ within region DIII, i.e. $x_c$ is the point at which

(7.27)\begin{equation} \mathcal{T}_g^{0}=\frac{1+\eta_D\log(\eta_D)/A}{KF_*}. \end{equation}

Using the form of $\mathcal {T}_g^{0}$ given by (7.17), we find that

(7.28)\begin{equation} x_c=x_1+\frac{1}{\mu}KF_*\log\left(\frac{1}{(1-KF_*)(1+\eta_D\log(\eta_D)/A)}\right).\end{equation}

7.2.3. Composite solution

We can form a leading-order, composite, steady-state solution as follows. Since $C=1$ at leading order in regions DIII-DaV’, the leading-order composite solution is simply

(7.29)\begin{equation} C^{Da}=C^{0}\left(\frac{x-s_0}{\delta}\right)+\hat{C}^{0}(\mu(x-s_0))-C^{\infty}. \end{equation}

The gas temperature is $O(1)$ in regions DI and DII, but $O(\eta _D)$ in regions DIII–DaV’. Since we have not computed the $O(\eta _D)$ correction in regions DI–DII, we assume that this correction is $\eta _D(1+\eta _D\log (\eta _D)/A)$ so that the composite solution is correct to $O(\eta _D)$ in regions DIII–DaV’. The composite solution is defined piecewise on either side of the transition layer DaIV. The centre of this transition layer is $x_c$, given by (7.28), which is defined so that there is a continuous transition between $T_s=\eta _D\mathcal {T}_g^{0}$ in region DIII and $T_s=\eta _D\bar {T}_g$ in region DaV. For the composite solution, we choose to define a slightly different changeover point $x_c^{*}=x_c+O(\eta _D/\mu )$, so that the composite solution is continuous at this point. We therefore use

(7.30)\begin{equation} T_g^{Da}=\begin{cases} T_g^{0}\left(\dfrac{x-s_0}{\delta}\right)+\hat{T}_g^{0}(\mu(x-s_0))-G_s & \\ \quad +\eta_D\mathcal{T}_g^{0}(\mu(x-x_1))-\exp\left(-\dfrac{\mu}{KF_*}(x-x_1)\right), & \text{if }x< x_c^{*},\\ \eta_D\bar{T}_g^{0}(\mu(x-x_c)), & \text{if }x>x_c^{*}, \end{cases} \end{equation}

which, from the form of $\mathcal {T}_g^{0}$ in (7.17), simplifies to

(7.31)\begin{equation} T_g^{Da}=\begin{cases} T_g^{0}\left(\dfrac{x-s_0}{\delta}\right)+\hat{T}_g^{0}(\mu(x-s_0))-G_s+\eta_D\left(1+\eta_D\dfrac{\log(\eta_D)}{A}\right), & \text{if }x< x_c^{*},\\ \eta_D\bar{T}_g^{0}(\mu(x-x_c)), & \text{if }x>x_c^{*}. \end{cases} \end{equation}

The solid temperature is similarly $O(1)$ in region DI, but $O(\eta _D)$ everywhere else. As in previous cases, we set the $O(\eta _D)$ correction in region DI to be $\eta _D(1+\eta _D\log (\eta _D)/A)$ so that the composite solution is accurate to $O(\eta _D)$ in regions DII–DaV’. We neglect the $O(\eta _D^{2})$ corrections in regions DII and DIII for this composite solution. We form the leading-order composite solution

(7.32)\begin{equation} T_s^{Da}= \begin{cases} T_s^{0}\left(\dfrac{x-s_0}{\delta}\right)+\eta_D\left(1+\eta_D\dfrac{\log(\eta_D)}{A}\right), & \text{if }x< x_c,\\ \eta_D\bar{T}_s^{0}(\mu(x-x_c)), & \text{if }x>x_c. \end{cases} \end{equation}

7.3. Case Db: $KF_*>1$

In the case $KF_*>1$, the heat flux in the gas is greater than in the solid. We therefore require the position of the boundary layer over which $T_s$ decays to $T_s^{in}$ to be at the far right of the domain, at $x_c=1$. This is so that the variable $Z<0$, and the solution (7.22a,b) is bounded as we leave the boundary layer ($Z\rightarrow -\infty$). This boundary layer is therefore region DbVI in figure 7.

In this case, we notice that as $Z\rightarrow -\infty$, both $\bar {T}_s^{0},\bar {T}_g^{0}\rightarrow b_1$, so that far from $x=1$, we have $T_g-T_s=O(\eta _D^{2})$. In region DIII, we have $T_g-T_s=O(\eta _D)$. We therefore cannot match directly with region DIII: there must be another region within which $T_s-T_g=O(\eta _D^{2})$. We denote this region DbIV.

7.3.1. Region DbIV

We note that the analysis of region DIII holds so long as $T_g-T_s\gg \eta _D^{2}$. As $T_g$ is decaying through region DIII, $T_g-T_s$ is decreasing. At the point where $T_g-T_s=O(\eta _D^{2})$, the advection (term 1 in (3.1)) must also enter be included in the dominant balance. From (7.17), we see that $T_g-T_s=\eta _D^{2}$ when

(7.33)\begin{equation} \exp\left(-\frac{1}{KF_*}z\right)=\eta_D, \end{equation}

i.e. at the point

(7.34)\begin{equation} z_*={-}KF_*\log(\eta_D). \end{equation}

Like $y_*$ given by (7.11), we note that $z_*=O(-\log (\eta _D))$ is large. Strictly speaking, there is therefore another region DbIII’ with length scale $O(-\log (\eta _D)/\mu )$ to the right of region DIII, but since this is a subregion of region DIII (the same equations hold, simply over a longer length scale), we do not analyse this region separately. Although $z_*$ is large, in the original variables, this is the point $x_2:=s_0+(y_*+z_*)/\mu \ll 1$, since $\mu \gg -\log (\eta _D)$. Thus $x_2$ remains within the $O(-\log (\eta _D)/\mu )$ boundary layer around $x=s_0$.

Region DbIV is the transition layer around the point $x_2$, in which all the terms 1, 3, 4 and 5 balance in the solid heat equation (3.1). Since $T_g-T_s=O(\eta _D^{2})$ in region DbIV, we make the change of variables

(7.35a,b)\begin{align} T_s=\eta_D\left(1+2\eta_D\frac{\log(\eta_D)}{A}+\eta_D\frac{\tau_s}{A}\right), \quad T_g=\eta_D\left(1+2\eta_D\frac{\log(\eta_D)}{A}+\eta_D\frac{\tau_g}{A}\right), \end{align}

so that at leading order both $T_s$ and $T_g$ are uniform through region DbIV. The heat exchange terms (3 and 4) now combine to be $O(\mu \eta _D^{2})$, and the heat lost to the chemical reaction (term 5) is likewise $O(\mu \eta _D^{2})$, since

(7.36)\begin{equation} Q(C,T_s)=\delta \mu\left(\eta_D^{2}Ce^{A+\tau_s}+O\left(\eta_D^{3}\right)\right). \end{equation}

As in regions DII and DIII, we require an $O(1/\mu )$ length scale to balance advection and heat transfer terms in the gas heat equation, and so we make the change of variables $x=x_2+Z/\mu$. As in region DIII, we note that since the reaction rate is small, $C=1+O(\eta _D)$, $f=F_*+O(\eta _D)$ are uniform at leading order throughout this region. Expanding in powers of $\eta _D$, taking

(7.37a,b)\begin{equation} \tau_s=\tau_s^{0}(Z)+O(\eta_D), \quad \tau_g=\tau_g^{0}(Z)+O(\eta_D), \end{equation}

we find that at leading order (2.12c) and (2.12d) become

(7.38a,b)\begin{equation} -\frac{\mathrm{d}{\tau_s^{0}}}{\mathrm{d}{Z}}=(\tau_g^{0}-\tau_s^{0})-\gamma A e^{A+\tau_s^{0}}, \quad KF_*\frac{\mathrm{d}{\tau_g^{0}}}{\mathrm{d}{Z}}={-}(\tau_g^{0}-\tau_s^{0}). \end{equation}

This transition region DbIV plays the same role as regions BIII, CIII and DaIV, in cases B, C and Da, respectively. It is in this region that term 1 enters the dominant balance of (3.1), while term 5 leaves the dominant balance. However, while $T_g$ is uniform at leading order over the transition region DbIV, we must solve for the variation, $\tau _g^{0}$, simultaneously with the solid temperature variation, $\tau _s^{0}$, via the coupled system (7.38a,b), rather than solving a single ODE for the solid temperature variation, as we did in regions BIII, CIII and DaIV.

While we cannot solve the system (7.38a,b) explicitly, we can study solutions using a phase-plane analysis. We first rescale $\bar {Z}=KF_*Z$ so that the system (7.38a,b) becomes

(7.39a,b)\begin{equation} \frac{\mathrm{d}{\tau_s^{0}}}{\mathrm{d}{\bar{Z}}}={-}c_1(\tau_g^{0}-\tau_s^{0})+c_2 e^{\tau_s^{0}}, \quad \frac{\mathrm{d}{\tau_g^{0}}}{\mathrm{d}{\bar{Z}}}={-}(\tau_g^{0}-\tau_s^{0}), \end{equation}

in terms of the two parameters $c_1=KF_*$ and $c_2=KF_*\gamma A e^{A}$. Solutions of this system (7.39a,b) are shown as trajectories in the phase-plane diagram in figure 8. The black lines shown in figure 8 are the nullclines

(7.40a,b)\begin{equation} \tau_g^{0}=\tau_s^{0} \quad \text{and} \quad \tau_g^{0}=\tau_s^{0}+\frac{c_2}{c_1}e^{\tau_s^{0}}, \end{equation}

on which $\mathrm {d}\tau _g/\mathrm {d}\bar {Z}=0$ and $\mathrm {d}\tau _s/\mathrm {d}\bar {Z}=0$, respectively. As $\bar {Z}\rightarrow -\infty$, both $\tau _s^{0}$ and $\tau _g^{0}$ must grow exponentially, in order to match with region DIII. From the phase-plane diagram in figure 8, we see that all such trajectories follow just above the $\mathrm {d}\tau _s^{0}/\mathrm {d}\bar {Z}=0$ nullcline as $\bar {Z}\rightarrow -\infty$. This is to be expected, since the advection term 1 must become small as $\bar {Z}\rightarrow -\infty$ in order to obtain the dominant balance between terms 3 and 5 in DIII. We observe in figure 8 that these trajectories stay just above the $\mathrm {d}\tau _s^{0}/\mathrm {d}\bar {Z}=0$ nullcline, before either peeling away to the left, or turning the corner and crossing both nullclines, eventually approaching curves with $\tau _g^{0}$ equal to a constant. In order to match to the right of region DbIV, the solution trajectory must not cross the nullclines, so that $\tau _g^{0}>\tau _s^{0}$ everywhere, and both $\tau _g^{0},\tau _s^{0}\rightarrow -\infty$ as $\bar {Z}\rightarrow \infty$.

Figure 8.

Phase plane for the steady-state system (7.39a,b). Nullclines are in black and trajectories in blue. The green dashed curve is the line $\tau _g^{0}=\tau _g^{limit}(\tau _s^{0})$ given by (7.42). We take physically relevant values $\gamma =A=1$ and $KF_*=1$, so that $c_1=2$ and $c_2=5.44$.

The trajectories do not cross the nullclines for any $\bar {Z}$, provided that the inequality

(7.41)\begin{equation} \left.\frac{\mathrm{d}{\tau_g^{0}}}{\mathrm{d}{\bar{Z}}}\right/\frac{\mathrm{d}{\tau_s^{0}}}{\mathrm{d}{\bar{Z}}}=\frac{\mathrm{d}{\tau_g^{0}}}{\mathrm{d}{\tau_s^{0}}}<\frac{\mathrm{d}}{\mathrm{d}\tau_s^{0}}\left(\tau_s^{0}+\frac{c_2}{c_1}e^{\tau_s^{0}}\right) \end{equation}

holds at all points on the solution trajectory. Substituting in the derivatives from (7.39a,b), we see that (7.41) reduces to the inequality

(7.42)\begin{equation} \tau_g^{0}>\frac{(c_1\tau_s^{0}+c_2e^{\tau_s^{0}})\left(1+\frac{c_2}{c_1}e^{\tau_s^{0}}\right)-\tau_s^{0}}{c_1+c_2e^{\tau_s^{0}}-1}=:\tau_g^{limit}. \end{equation}

The trajectory $\tau _g^{0}=\tau _g^{limit}(\tau _s^{0})$ is plotted as the dashed green curve in figure 8, and is seen to divide the trajectories crossing the nullclines from those remaining above them. As $\tau _s^{0}\rightarrow \infty$, we have

(7.43)\begin{equation} \tau_g^{limit}\sim\tau_s^{0}+\frac{c_2}{c_1} e^{\tau_s^{0}}+\frac{1}{c_1}+\cdots. \end{equation}

Therefore, as long as $\tau _g^{0}>\tau _s^{0}+\gamma \exp (A(1+\tau _s^{0}))+1/c_1$ in the limit as $\bar {Z}\rightarrow -\infty$ (which must be true from the matching with region DIII), the trajectories never cross the nullclines.

We notice that, on these trajectories, $\tau _s^{0}$ becomes negative as $\bar {Z}\rightarrow \infty$, and so the exponential term in (7.39a,b), representing the heat consumed by the chemical reaction, becomes exponentially small in this limit. As $\bar {Z}\rightarrow \infty$, we therefore see by combining (7.39a,b) that $\tau _s^{0}-c_1\tau _g^{0}$ approaches a constant, so that the trajectories approach straight lines, as can be seen in figure 8. We also notice that some such trajectories follow the curve $\tau _g^{limit}(\tau _s^{0})$ for a long time before peeling off to the $\tau _s^{0}-c_1\tau _g^{0}$ behaviour. The curve $\tau _g^{limit}(\tau _s^{0})$ approaches the nullclines as both $\tau _s^{0}$ and $\tau _g^{0}$ decrease below zero, and the nullclines become closer together. Thus on these trajectories remaining close to $\tau _g^{limit}(\tau _s^{0})$, both the exponential chemical reaction term in (7.39a,b) and the heat-transfer terms $\tau _g^{0}-\tau _s^{0}$ are very small, so that both $\mathrm {d}\tau _g^{0}/\mathrm {d}\bar {Z}$ and $\mathrm {d}\tau _s^{0}/\mathrm {d}\bar {Z}$ are small. The solution we are concerned with must lie on one of these very slowly moving trajectories, since $\tau _g^{0}-\tau _s^{0}$ cannot grow until we reach the boundary-layer region DbVI. Thus the transition layer region DbIV must actually extend fully into region DbV, which may equally be considered a subregion of either DbIV or DbVI.

7.3.2. Regions DbV and DbVI

To the right of region DbIV, the heat lost to the chemical reaction, term 5 in (3.1), is not included in the dominant balance, so the system is purely advective, until the heat transfer between phases becomes important again in region DbVI. Thus in region DbV all dependent variables are constant at leading order, namely

(7.44)\begin{equation} T_s=T_g=\eta_D\left(1+2\eta_D\frac{\log(\eta_D)}{A}\right), \quad f=F_*, \quad C=1, \end{equation}

where we have fixed the constants $f=F_*$ by matching with region DbIV, and $C=1$ by matching with region DbVI.

Using this behaviour in region DbV, we may complete our solution in region DbVI. By matching between regions DbV and DbVI, we find $b_1=1+2\eta _D\log (\eta _D)/A,$ and applying the boundary condition $\bar {T}_s^{0}=\bar {T}^{in}$ at $Z=0$, we find

(7.45)\begin{equation} b_2=\frac{1}{KF_*}(\bar{T}^{in}-a_1)=\frac{1}{KF_*}\left(\bar{T}^{in}-1-2\eta_D\frac{\log(\eta_D)}{A}\right). \end{equation}

Thus, in region DbVI

(7.46)$$\begin{gather} \bar{T}_s^{0}=1+2\eta_D\frac{\log(\eta_D)}{A}+\left(\bar{T}^{in}-1-2\eta_D\frac{\log(\eta_D)}{A}\right)\exp\left(\frac{KF_*-1}{KF_*}Z\right), \end{gather}$$

(7.47)$$\begin{gather}\bar{T}_g^{0}=1+2\eta_D\frac{\log(\eta_D)}{A}+\frac{1}{KF_*}\left(\bar{T}^{in}-1-2\eta_D\frac{\log(\eta_D)}{A}\right)\exp\left(\frac{KF_*-1}{KF_*}Z\right), \end{gather}$$

completing the analysis of case Db.

7.3.3. Composite solution

We can therefore construct a leading-order, steady-state composite solution for case Db. Since $C$ and $f$ are everywhere $O(1)$, and are constant to leading order in regions DIII–DbVI, the composite expansions are simply

(7.48)$$\begin{gather} C^{Db}=C^{0}\left(\frac{x-s_0}{\delta}\right)+\hat{C}^{0}\left(\mu(x-s_0)\right)-C^{\infty}, \end{gather}$$

(7.49)$$\begin{gather}f^{Db}=f^{0}\left(\frac{x-s_0}{\delta}\right)+\hat{f}^{0}\left(\mu(x-s_0)\right)-F_s. \end{gather}$$

Both $T_s$ and $T_g$ are $O(1)$ at $x=s$, but become $O(\eta _D)$ as we move further into the domain. As in the previous cases there is freedom to choose a composite solution. Since $T_s\sim \eta _D$ is uniform through regions DII–DbV, we use the composite solid temperature

(7.50)\begin{equation} T_s^{Db}=T_s^{0}\left(\frac{x-s_0}{\delta}\right)+\eta_D\bar{T}_s^{0}\left(\mu(x-1)\right), \end{equation}

so that $T_s$ is accurate with error $O(\eta _D)$ in region DI, with error $O(\eta _D^{2})$ in regions DbIV–DbVI, but with $O(\eta _D\log (\eta _D))$ error in regions DII–DbIII. Similarly, we choose the composite gas temperature

(7.51)\begin{equation} T_g^{Db}=T_g^{0}\left(\frac{x-s_0}{\delta}\right)+\hat{T}_g^{0}\left(\mu(x-s_0)\right)-G_s+\eta_D\bar{T}_g^{0}\left(\mu(x-1)\right), \end{equation}

accurate to $O(\eta _D)$ in regions DI–DII, to $O(\eta _D\log (\eta _D))$ in region DIII and to $O(\eta _D^{2})$ in regions DbIV–DbVI.

8.1. Region EI

As usual we look for a boundary layer, region EI, of width $O(\delta )$ at $x=s$, setting $x=s+\delta X$. Unlike the previous cases, we now have heat transfer between the solid and gas phases on the $O(\delta )$ length scale, and so the solid and gas equations in region EI do not decouple. At leading order in $\delta$, the equations in region EI are quasi-steady, so that we have $f=f_*+(1+\dot {s})(C-C_*)+O(\delta )$ here. To leading order in $\delta$, the equations for $C$, $T_s$ and $T_g$ are

(8.1a)$$\begin{gather} -(1+\dot{s})\frac{\partial{C}}{\partial{X}}={-}C\exp\left(A\left(1-\frac{1}{T_s}\right)\right), \end{gather}$$

(8.1b)$$\begin{gather}-(1+\dot{s})C\frac{\partial{T_s}}{\partial{X}}=\frac{\partial{{}^{2}T_s}}{\partial{X^{2}}}+M(T_g-T_s)-\gamma C\exp\left(A\left(1-\frac{1}{T_s}\right)\right), \end{gather}$$

(8.1c)$$\begin{gather}K\frac{\partial}{\partial X}(fT_g)=T_sC\exp\left(A\left(1-\frac{1}{T_s}\right)\right)-M(T_g-T_s), \end{gather}$$

with boundary conditions

(8.1d)\begin{equation} C=C_*, \quad T_s=T_g=\frac{\rho}{s} \quad \text{at }X=0. \end{equation}

As in cases A–D, we will see (in § 8.2 below) that the main source of error in our approximation in this region EI is the switchback error due to the non-zero temperatures in region EII. Thus, although we have only neglected terms of order $\delta$ in (8.1), the error in the solution is in fact greater than $O(\delta )$. We note that all the physical processes come into the dominant balance in this boundary layer EI, except for the time derivatives. Like the region-I problem in cases A–D, (8.1) is quasi-steady.

The system (8.1) is closed by applying appropriate matching conditions as $X\rightarrow \infty$, but it is not immediately clear what the correct conditions are. In the far field, as $X\rightarrow \infty$, we require that the chemical reaction rate becomes small in order that the rate of change of $C$ becomes small. This means we expect $T_s$ to be small as we leave the boundary layer. Using the far-field conditions

(8.2)\begin{equation} T_s\rightarrow 0, \quad C\rightarrow C^{\infty}(t) \quad \text{as }X\rightarrow\infty, \end{equation}

as in cases C and D, where $C^{\infty }(t)$ is to be determined, the system (8.1) may be shown to be correctly specified so long as

(8.3)\begin{equation} K\left(f_*+(1+\dot{s})(C^{\infty}-C_*)\right)<(1+\dot{s})C^{\infty}, \end{equation}

where, as always, $s$ is found as part of the solution of the region-I problem. Since $\dot {s}$ and $C^{\infty }$ are functions of time, we note that (8.3) may hold for certain times, and not others. This condition (8.3) is precisely the requirement that, in the far field of the moving boundary layer EI, the heat flux left-to-right in the gas is less than the heat flux right-to-left in the solid. Thus the far-field conditions (8.2) can only be imposed if the overall direction of heat flow is right-to-left within the moving boundary layer, i.e. into the boundary layer. In this case, we note that the conditions (8.2) ensure that $T_g\rightarrow 0$ as $X\rightarrow 0$ also, and so $T_s-T_g\rightarrow 0$ in the far-field limit.

We restrict ourselves to the case $K<1$, i.e. $c_{p,g}< c_{p,s}$, which we expect to be physically relevant in general (and is certainly relevant for the silicon furnace (Luckins et al. Reference Luckins, Oliver, Please, Sloman and Van Gorder2022)). In this case we note that

(8.4)\begin{equation} C^{\infty}-K(C^{\infty}-C_*)=C^{\infty}(1-K)+KC_*>0 \end{equation}

is positive. The inequality (8.3) may therefore be rearranged to

(8.5)\begin{equation} \dot{s}>{-}1+\frac{Kf_*}{C^{\infty}-K(C^{\infty}-C_*)}=:U(C^{\infty}). \end{equation}

Thus the matching conditions (8.2) are appropriate if the position of the free boundary does not decrease too rapidly.

We note that (regardless of the size of the convective heat-transfer coefficient $\mu$) for physically relevant solutions we require $\dot {s}>-1$, since the fastest the free boundary can move inwards is the speed of the solid material. Since $K>1$, we see from (8.5) that $U>-1$. Thus we expect both the situations when (8.5) does and does not hold to be possible, and physically relevant. If (8.5) does not hold, so that $\dot {s}< U$, the far-field conditions (8.2) cannot be imposed. The overall direction of heat transfer in the far field of the boundary layer is now left-to-right, and so we must understand the decay behaviour of the region EI problem, in order to choose the correct far-field conditions. Since all physical effects are included in the region EI equations (8.1), in the far field we will simply have some of these effects dropping out of the dominant balance as the length scale becomes large. We will find the appropriate far-field behaviour of the region EI model by examining the system on greater length scales.

8.2. The $O(1)$ domain

We now investigate the behaviour in the majority of the domain for either of the cases $\dot {s}>U$ or $\dot {s}< U$. On an $O(1)$ length scale, we note from the mass equations,(2.12a)–(2.12b), that the chemical reaction rate must be small, at most $O(\delta )$. Since the chemical reaction rate is small, $T_s$ must be small. We suppose that $T_s=O(\theta )$ for some small $\theta \ll 1$ to be determined. From the (dominant) $M(T_g-T_s)$ term in either energy equation (2.12c) or (2.12d), we must therefore have $T_g=O(\theta )$, as well. By balancing the advection or chemical reaction terms in the solid heat equation, we see that the difference $T_g-T_s$ must be at most order

(8.6)\begin{equation} \epsilon:=\max\left(\delta\theta,\exp\left(-\frac{A}{\theta}\right)\right)\ll\theta,\end{equation}

so that the two temperatures are equal at leading order in $\theta$. Thus the model (2.12) reduces to a single-temperature system over $O(1)$ length scales.

Specifically, we write $T_s=\theta T+\epsilon T_s^{1}$ and $T_g=\theta T+\epsilon T_g^{1}$. Adding the two heat equations (2.12c) and (2.12d) together, the convective heat transfer between the phases cancel, and we obtain a single equation for the leading-order temperature $T$, namely

(8.7)\begin{equation} \delta\theta\left(\frac{\partial}{\partial t}(CT)+\frac{\partial}{\partial x}\left((Kf-C)T\right)\right)={-}\gamma C\exp\left(A\left(1-\frac{1}{\theta T}\right)\right). \end{equation}

Here, we have absorbed the term for the heat transfer between phases due to the mass transfer into the advection term, by using the solid mass equation (2.12a). We have neglected the heat diffusion term, which is always smaller than the advection term on length scales greater than $O(\delta )$. Only the two remaining terms in (8.7) – heat advection and heat consumed by chemical reaction – can possibly enter the dominant balance, and we must choose the appropriate $\theta$ to ensure the correct balance. We will see that the $O(\epsilon )$ correction is smaller than the corrections due to other effects, and so the single-temperature model (8.7) is indeed valid. We note that (8.7) is a first-order, hyperbolic equation for $T$, and so our entire system is hyperbolic on the $O(1)$ length scale: the $C$ characteristic curves travel right-to-left, the $f$ characteristic curves travel left-to-right and the $T$ characteristic curves from (8.7) travel in direction $Kf-C$. The correct choice of scaling $\theta$, and so the balance in (8.7), depends on the direction of the $T$ characteristics, and hence on the sign of $Kf-C$.

We notice that the chemical reaction rate must be at most order $\delta \theta$ over this $O(1)$ length scale in order to be balanced by the advection term in (8.7). From the equations of conservation of mass, (2.12a)–(2.12b), this means that both $C$ and $f$ are uniform at leading order in $\theta$ over $O(1)$ domains. We will later see that there are no regions of $O(1)$ variation in $C$ or $f$ outside of region EI, and thus both

(8.8a,b)\begin{equation} C=1+\epsilon_C\hat{C} \quad \text{and} \quad f=f^{\infty}(t)+\epsilon_C\hat{f} \end{equation}

are uniform at leading order over the entire $O(1)$ domain $x>s(t)$. Here, the scaling $\epsilon _C\ll 1$ of the correction terms is yet to be determined. This depends on which boundary the $T$ characteristic curves originate: we will find $\epsilon _C\sim \delta ^{-1}Q(T^{in})=o(T^{in})$ for when these originate at the $x=1$ boundary, and $\epsilon _C=\theta _*^{2}$ when these originate at $x=s$ (in either case $\epsilon _C\gg \epsilon$, the error in the single-temperature approximation, defined by (8.6)). In (8.8a,b) we have fixed $C^{\infty }=1$ at leading order using the boundary condition at $x=1$, and

(8.9)\begin{equation} f^{\infty} (t)=f_*+(1+\dot{s}(t))(1-C_*), \end{equation}

by matching with the boundary layer EI as $X\rightarrow \infty$. We note that the leading-order gas flux $f^{\infty }$ is spatially uniform but not constant in time. Thus the direction of the characteristic curves for $T$ are uniform over the entire $O(1)$ domain, and are given by

(8.10)\begin{equation} Kf-C=Kf^{\infty}(t)-1 +O(\epsilon_C). \end{equation}

Since the characteristic curve slopes depend only on time and do not vary spatially, we do not expect characteristic curves to cross anywhere in the domain.

We will have different behaviour (and different $\theta$ scalings) in regions where the characteristics of (8.7) originate at different boundaries of the $(x,t)$ domain. We note that these $T$-characteristic curves, with gradient $Kf^{\infty }-1$, can originate at $x=1$ only if $Kf^{\infty }<1$, or equivalently only if

(8.11)\begin{equation} \dot{s}< V:=\frac{1-K(f_*+(1-C_*))}{K(1-C_*)}. \end{equation}

Also, characteristic curves can originate at $x=s(t)$ only if $Kf^{\infty }-1>\dot {s}$, or equivalently (since we are assuming throughout that $K<1$) only if

(8.12)\begin{equation} \dot{s}< U,\end{equation}

where (in a slight abuse of notation) we write $U:=U|_{C^{\infty }=1}$ for the function $U(C^{\infty })$ defined as in (8.5), evaluated at $C^{\infty }=1$. Thus characteristic curves in the outer domain originate at $x=s$ if and only if the heat transfer is left-to-right in the far field of the (moving) boundary layer EI. We note that $U$ and $V$ have opposite signs: if (as in case Ea) $K(f_*+1-C_*)<1$ then $U<0< V$, whereas if (as in case Eb) $K(f_*+1-C_*)>1$ then $U>0>V$. For either sign of $K(f_*+1-C_*)-1$, there are a variety of solution structures depending on which, if any, of (8.11) and (8.12) hold, since (8.11) and (8.12) are neither mutually exclusive nor mutually inclusive. In particular, the steady-state solutions in the cases Ea and Eb, are illustrated in figure 9(a,b). The characteristic diagrams for the $T$ problem (8.7), showing the different structures possible (in the general, non-steady case), are presented in figure 10. As $\dot {s}$ changes continuously in time, there are other structures possible than those shown here.

Figure 10.

Characteristic diagrams in the cases $V>0>U$ and $V<0< U$, for different sizes of $\dot {s}$. The red characteristics mark the boundaries between regions with the characteristics originating at different boundaries of the domain, or at the initial data. In general, we expect to have different dominant behaviour on either side of these bounding characteristics.

We consider regions with characteristic curves originating at $x=1$ and $x=s$ separately, as the temperature scaling $\theta$, and so too the solution behaviour, is different in these two cases. To fully determine the non-steady solutions, we should also investigate the case when characteristic curves originate at the $t=0$ axis in the diagrams of figure 10. The solution in this case is determined by the initial conditions, in particular on the magnitude of the temperature, but for simplicity we do not discuss this here.

8.2.2. The $T$ characteristics originate at $x=s(t)$

In a region where the characteristics originate at $x=s$, then the temperature scaling $\theta =\theta _*$ must be determined by the temperature leaving the boundary layer EI. We suppose that the temperature is approximately constant, and make the changes of variables $T=\theta _*(1+\theta _*\hat {T})$ and $C=1+\theta _*^{2}\hat {C}$ and $f=f^{\infty }+\theta _*^{2}\hat {f}$. In (8.7), we see that $\theta _*$ must be determined by a balance between the chemical reaction term and the heat advection, and we therefore require that $\theta _*$ is the smallest positive solution of

(8.13)\begin{equation} \delta\theta_*^{2}e^{A/\theta_*}=1,\end{equation}

namely,

(8.14)\begin{equation} \theta_*= \frac{-A}{2W_{{-}1}\left(-\frac{A}{2}\sqrt{\delta}\right)} =\frac{A}{\log\left(\dfrac{1}{\delta}\right)}\left(1+\frac{\log(\log(\delta)^{2})}{\log(\delta)}+\cdots\right), \end{equation}

and we note that $\theta _*\rightarrow 0$ as $\delta \rightarrow 0$. As an aside, we also note that the temperature through the majority of the domain in case Db is

(8.15)\begin{equation} \eta_D\left(1+2\eta_D\frac{\log(\eta_D)}{A}\right)= \frac{A}{\log\left(\dfrac{1}{\delta\mu}\right)}\left(1+\frac{\log\left(\log(\delta\mu)^{2}\right)}{\log(\delta\mu)}+\cdots\right), \end{equation}

which is the same as $\theta _*$ to leading order, since we had $\mu \ll \delta ^{-1}$ in case Db. Indeed all of $\eta _B,\eta _C,\eta _D\sim A/(-\log (\delta ))$ at leading order, although the correction terms differ, and, being logarithmic, are not particularly small.

At leading order in $\theta _*$, (2.12a)–(2.12b) and (8.7) become

(8.16a,b)\begin{gather} \frac{\partial{\hat{C}}}{\partial{t}}-\frac{\partial{\hat{C}}}{\partial{x}}={-}\exp\left(A\left(1+\hat{T}\right)\right), \quad \frac{\partial{\hat{f}}}{\partial{x}}=\exp\left(A\left(1+\hat{T}\right)\right), \end{gather}

(8.17)\begin{gather} \frac{\partial{\hat{T}}}{\partial{t}}+(Kf^{\infty}-1)\frac{\partial\hat{T}}{\partial x}={-}\gamma \exp\left(A\left(1+\hat{T}\right)\right).\end{gather}

Solving (8.17) by the method of characteristics, we find the characteristic curves, parameterised by $\xi$ and the time $t_0$ at which $x=s_0$,

(8.18)$$\begin{gather} t(\xi,t_0)=\xi+t_0, \end{gather}$$

(8.19)$$\begin{gather}x(\xi,t_0)=s(t_0)+\int_{\hat{\xi}=0}^{\xi}\left[K\left(f_*+(1+\dot{s}(\hat{\xi}+t_0))(1-C_*)\right)-1\right]\,\mathrm{d}\hat{\xi}. \end{gather}$$

Along these characteristic curves,

(8.20)\begin{equation} \hat{T}(\xi)={-}\frac{1}{A}\log\left(Ae^{A}\gamma \xi\right). \end{equation}

Here, we have fixed the constant of integration so that $\hat {T}$ blows up as $\xi \rightarrow 0$, i.e. as the characteristic curve approaches $x=s_0$. We may then find $\hat {C}$ and $\hat {f}$ from (8.16a,b), with the requirements that $\hat {C}=0$ at $x=1$ and $\hat {f}=0$ at $x=s(t)$.

In particular, in steady state, the full solution is

(8.21a,b)\begin{equation} T=\theta_*-\frac{\theta_*^{2}}{A}\log\left(\frac{A\gamma e^{A}}{KF^{\infty}-1}(x-s)\right), \quad C=1+\theta_*^{2}\frac{Kf^{\infty}-1}{A\gamma}\log\left(\frac{x-s}{1-s}\right), \end{equation}

and, as always in steady state, $f=f_*+C-C_*$. Specifically, this is the region EbII solution in the steady-state case Eb, illustrated in figure 9(b), since in case Eb, characteristic curves originate at $x=s$ and end at $x=1$.

Given this steady-state solution, we now understand the far-field behaviour, and the appropriate boundary conditions for the region-EbI problem. Since (8.1a) is first order in $C$, we prescribe the leading-order far-field behaviour

(8.22)\begin{equation} C\rightarrow 1 \quad \text{ as }X\rightarrow\infty, \end{equation}

while for the $T_s$ boundary condition, we can impose the correct far-field balance

(8.23)\begin{equation} (Kf^{\infty}-1)\frac{\mathrm{d}{T_s}}{\mathrm{d}{X}}+\gamma\exp\left(A\left(1-\frac{1}{T_s}\right)\right)\rightarrow 0 \quad \text{ as }X\rightarrow\infty. \end{equation}

We note that when characteristics of the outer problem move into region EI, the leading-order problem (8.1) is accurate to order $\theta \ll 1$, because the temperature (which is order $\theta$) outside the boundary layer induces a switchback correction within the boundary layer. In particular, in the steady-state case Ea, we found $\theta =T^{in}$, and thus the switchback error in case Ea is $O(T^{in})$. Conversely, when $\dot {s}< U$ and so the characteristics of the outer problem originate at $x=s$, the boundary conditions are given by (8.22)–(8.23), and we do not have a switchback correction due to the outer region. Instead the main source of error is the $O(\theta _*^{2})$ error in the far-field boundary condition (8.22) for $C$.

8.3. Region EbIII

In situations where characteristics start at $x=s$ and end at $x=1$ (such as in the steady-state case Eb in figure 9b), we must introduce a boundary layer at $x=1$ of width $O(\delta )$, in order to impose the boundary condition $T=T^{in}$ at $x=1$. (In steady state, we denote this boundary-layer region EbIII.) This is the region, analogous to regions BIV and CIV, in which the heat transferred from the gas to the solid heats the solid from its input temperature, $T^{in}$. We make the change of variables $x=1-\delta Y$, and set

(8.24a,b)\begin{equation} T_s=\theta_*\bar{T}_s, \quad T_g=\theta_*\bar{T}_g, \end{equation}

for $\bar {T}_s<1$. The chemical reaction rate is then

(8.25)\begin{equation} \bar{Q}\sim \exp\left(-\frac{A}{\theta_*\bar{T}_s}\right)\ll e^{{-}A/\theta_*}=\delta\theta_*^{2}, \end{equation}

by definition of $\theta _*$. Thus as the chemical reaction rate is small, so that $C=1+O(\theta _*^{2})$ and $f=f^{\infty }+O(\theta _*^{2})$ are uniform to leading order in the boundary layer. We expand $\bar {T}_s=\bar {T}_s^{0}(Y,t)+O(\theta _*)$ and $\bar {T}_g=\bar {T}_g^{0}(Y,t)+O(\theta _*)$, and at leading order in the heat equations (2.12c) and (2.12d) we find the quasi-steady equations

(8.26a)\begin{equation} \frac{\partial{\bar{T}_s^{0}}}{\partial{Y}}=\frac{\partial{{}^{2}\bar{T}_s^{0}}}{\partial{Y^{2}}}+M(\bar{T}_g^{0}-\bar{T}_s^{0}), \quad Kf^{\infty}\frac{\partial{\bar{T}^{0}_g}}{\partial{Y}}=M(\bar{T}_g^{0}-\bar{T}_s^{0}), \end{equation}

(the time dependence is through $f^{\infty }(t)$) with boundary conditions

(8.26b)\begin{equation} \bar{T}_s^{0}=\frac{T^{in}}{\theta_*} \quad \text{ at }Y=0, \quad \bar{T}_s^{0},\bar{T}_g^{0}\rightarrow 1 \quad \text{ as }Y\rightarrow\infty. \end{equation}

The solution of (8.26) is

(8.27a,b)\begin{equation} \bar{T}_s^{0}=1-\left(1-\frac{T^{in}}{\theta_*}\right)\exp\left(\chi(t) Y\right), \quad \bar{T}_g^{0}=1+\frac{\chi(t)-1}{Kf^{\infty}}\left(1-\frac{T^{in}}{\theta_*}\right)\exp\left(\chi(t) Y\right), \end{equation}

where the function of time $\chi (t)$ is given by

(8.28)\begin{equation} \chi(t)=\frac{1}{2}\left(\frac{M}{Kf^{\infty}}+1\right)-\frac{1}{2}\sqrt{\left(\frac{M}{Kf^{\infty}}-1\right)^{2}+4M}. \end{equation}

8.4. Composite solutions for cases Ea and Eb

In case Ea, the characteristic curves originate at $x=1$, and so all dependent variables are uniform to leading order outside region EI. In this case the steady-state composite solution is

(8.29a–c)\begin{equation} C^{Ea}=C^{0}\left(\frac{x-s_0}{\delta}\right), \quad T_s^{Ea}=T_s^{0}\left(\frac{x-s_0}{\delta}\right)+T^{in},\quad T_g^{Ea}=T_g^{0}\left(\frac{x-s_0}{\delta}\right)+T^{in}. \end{equation}

The steady-state composite solution in case Eb, when characteristics originate at $x=s$, is

(8.30a)$$\begin{gather} C^{Eb}=C^{0}\left(\frac{x-s_0}{\delta}\right), \end{gather}$$

(8.30b)$$\begin{gather}T_s^{Eb}=T_s^{0}\left(\frac{x-s_0}{\delta}\right)+\theta_*\bar{T}_s^{0}\left(\frac{1-x}{\delta}\right)-\theta_*, \end{gather}$$

(8.30c)$$\begin{gather}T_g^{Eb}=T_g^{0}\left(\frac{x-s_0}{\delta}\right)+\theta_*\bar{T}_g^{0}\left(\frac{1-x}{\delta}\right)-\theta_*. \end{gather}$$

8.5. Additional structures in some time-dependent cases

Some time-dependent solutions will have the same structure as the steady-state cases Ea and Eb. Additional structures are also possible in situations where characteristic curves originate at both $x=1$ and $x=s$. The boundary between these regions is the characteristic curve (8.19) which comes from $x=s$ at the point where $\dot {s}=\min (U,V)$ and subsequently $\dot {s}<\min (U,V)$ (shown in red in figure 10). We denote this dividing characteristic curve by $x_0(t)$. Around this characteristic curve, we expect there to be a thin layer, of width $\sqrt {\delta }$ to balance the time-derivative and diffusion terms, in which diffusion smooths out the discontinuity between the two regions. A schematic of the solution structure in this case is shown in figure 11.