1. Introduction
This paper is motivated by applications in which flow constriction can arise due to wall-mounted ice formation, such as aero-engine icing or flow blockages in large-scale plumbing systems (Gent, Dart & Cansdale Reference Gent, Dart and Cansdale2000; Elliot & Smith Reference Elliot and Smith2017). In these settings, flow constriction is undesirable and must either be prevented or treated subsequently. Total flow blockage due to ice formation, for instance, can lead to a pipe fracture due to the slight volumetric expansion accompanying the transition of water to ice. Theoretical modelling informed by these applications is therefore of potential practical value, especially when it yields physical estimates for quantities such as typical freezing times, the spatial extent of ice formation, and the pressure and shear levels encountered in flowing water.
Some other flow configurations that are significantly affected by ice growth include ship motion (Norde, van der Weide & Hoeijmakers Reference Norde, van der Weide and Hoeijmakers2022), aero-space glaciated icing (Trontin & Villedieu Reference Trontin and Villedieu2018), ice accumulation on moving vehicles with a flexible surface (Li et al. Reference Li, Jin, Chen, Han and Cong2020) and experimental work on ice crystal growth on axisymmetric bodies (Bucknell et al. Reference Bucknell, McGilvray, Gillespie, Jones, Reed and Collier2018). The freezing of water droplets also receives much attention. For example, Jingru et al. (Reference Jingru, Jingbin, Zhongwei, Kang and Haojun2024) describe experiments concerned with the freezing process of water droplets for air jet technology, while theoretical analyses and accompanying computations are described by Elliot & Smith (Reference Elliot and Smith2017), and Gorin, Bonn & Kellay (Reference Gorin, Bonn and Kellay2022) examine droplet impact onto a cold surface. Experiments in the context of droplet impact were performed by Hammond et al. (Reference Hammond, Quero, Ivey, Miller, Purvis, McGregor and Tan2012).
Numerous studies have investigated flow through channels or pipes with fixed constrictions (i.e. not dependent on the flow behaviour), primarily motivated by applications in physiological systems, as well as mechanical and engine dynamics. Several examples especially concerned with flow structure are due to Smith (Reference Smith1979) on flow separation, Smith & Duck (Reference Smith and Duck1980) who modelled large-scale separating flow in a channel, Sobey (Reference Sobey2000) on interactive boundary layers in channels and Dennis & Smith (Reference Dennis and Smith1980) on channel flow simulations. Beyond rigid geometries, other research has addressed the effects of elastic walls, where flexible constrictions arise in channel flow (Gajjar & Sibanda Reference Gajjar and Sibanda1996; Davies & Carpenter Reference Davies and Carpenter1997; Guneratne & Pedley Reference Guneratne and Pedley2006; Pruessner & Smith Reference Pruessner and Smith2015). Related phenomena have also been explored in the context of phase change: for instance, Moore, Mughal & Papageorgiou (Reference Moore, Mughal and Papageorgiou2017) considered ice formation beneath a thin liquid film subject to oncoming airflow in a boundary layer, though the resulting ice does not necessarily act as a constriction. More directly coupling fluid flow and phase change, Jolley, Dang & Smith (Reference Jolley, Dang and Smith2024) developed and analysed a model which describes the melting of a wall-mounted ice hump exposed to warm incoming flow, deep within a boundary layer.
Of particular relevance to this paper is the work of Weigand & Beer (Reference Weigand and Beer1992), who examined wall-mounted freezing of fluid in laminar flow within a parallel plate channel. Their analysis revealed that the axial position of the ice boundary converges towards the channel centreline and scales with the square root of time, exhibiting characteristics similar to the classical Stefan problem (Fowler Reference Fowler1997). In their formulation, however, the position of the ice boundary is decoupled from the oncoming channel flow via reductive techniques, allowing for a more tractable analysis of the freezing dynamics. Related studies that explore the effect of ice formation on fluid flow in confined geometries include the axi-symmetric pipe model of Ozisik & Mulligan (Reference Ozisik and Mulligan1969), the investigation by Kikuchi et al. (Reference Kikuchi, Shigemasa, A. and Ogata1986) into the influence of a warm-fluid region upstream of ice formation and the work of Sampson & Gibson (Reference Sampson and Gibson1981) which examines nozzle blockage arising from ice growth. In contrast to the aforementioned studies, there appear to be relatively few, if any, contributions in the open literature that account for the influence of an oncoming channel flow on the evolution of wall-mounted ice further downstream. This hence motivates the present work.
In this work, we develop a continuum model describing the interaction between fluid flow and two constricting, wall-mounted ice boundaries in a semi-infinite, two-dimensional channel. Focusing on high-Reynolds-number flow, the fluid dynamics is governed by the boundary layer equations, which are coupled to a thermal equation that describes the evolution of fluid temperature throughout the channel. The positions of the unknown constricting ice boundaries are governed by Stefan conditions. Via direct numerical simulations of the full model, we demonstrate that ice forming on the channel walls can subsequently grow inwards, eventually leading to severe constriction of the channel flow. Conversely, in other parameter regimes, there is no ice formation. Using both a numerical and asymptotic approach, we identify regions of parameter space in which ice formation does or does not occur.
The remainder of the paper is constructed as follows. In § 2, we state the equations governing the fluid flow, associated temperature and ice boundaries, as well as necessary boundary and inlet conditions. Following this, exemplar numerical solutions of the model are presented in § 3, which illustrate the presence of constricting ice boundaries suggesting that the channel gap width reduces to zero in a finite time (
$t=t_c$
) and at a fixed point longitudinally. An asymptotic analysis of the model about this critical point is then presented in § 4, whereby the model variables are scaled and expanded in constant powers of the small quantity
$t_c-t\ll 1$
; interestingly, these exponents are shown to depend on the incoming far-field fluid and thermal dynamics. Furthermore, the values of these exponents are used to indicate when severe channel constriction due to ice growth is or is not expected, in terms of the model parameters. Lastly, we discuss the model behaviour in § 5 and outline some potential avenues of further research.
2. Model development
In this section, we develop a dimensionless model which describes the interaction between fluid flow and a constricting ice boundary. The independent variable
$t$
denotes dimensionless time, whereas
$x^*\geq 0$
and
$0\leq y^* \leq f\leq 1$
describe the Cartesian longitudinal and lateral coordinates of a two-dimensional gap within a channel, where
$y^*=f$
represents the upper gap boundary.
The upper channel wall is fixed at
$y^* = 1$
so that
$f=1$
corresponds to the absence of ice formation on the wall. When
$f\lt 1,$
then
$y^*=f$
denotes the position of the ice boundary. We assume symmetry of the model solutions about
$y^*=0,$
such that the lower gap boundary and lower channel wall are at
$y^*=-f$
and
$y^*=-1,$
respectively. The contained channel fluid flow has dimensionless pressure
$p^*(x^*,\,y^*,\,t)$
and
$x^*\hbox{-}$
and
$y^*\hbox{-}$
directional dimensionless velocities
$u^*(x^*,\,y^*,\,t)$
and
$v^*(x^*,\,y^*,\,t).$
The dimensionless fluid temperature is denoted by
$\theta ^*(x^*,\,y^*,\,t)$
and the fluid freezing temperature is normalised to
$\theta ^*=-1.$
As seen later, the characteristic channel length scale is sufficiently large to allow the ice temperature to be related linearly with
$y^*.$
Given that the fluid freezing temperature is
$ \theta ^* = -1$
and ice forms exclusively on a wall maintained at this temperature, the ice temperature is consequently fixed at
$-1$
. Throughout this paper, we assume that ice has the same density as the fluid. We assume that any ice forming on the channel walls is pure and that phase change at the ice–water interface is suitably governed by a classical Stefan condition. The fluid flow is taken to be laminar and quasi-steady, with the dominant time scale arising from the comparatively slow evolution of the moving ice–water boundary.
2.1. Governing equations
We assume that the channel fluid is incompressible and Newtonian, and adopt the dimensionless steady Navier–Stokes equations with incompressibility condition
\begin{align} u^*\frac {\partial u^*}{\partial x^*}+v^*\frac {\partial u^*}{\partial y^*}&=-\frac {\partial p^*}{\partial x^*}+{\textit{Re}}^{-1}\bigg (\frac {\partial ^2 u^*}{\partial x^{*2}}+\frac {\partial ^2 u^*}{\partial y^{*2}}\bigg ), \end{align}
\begin{align} u^*\frac {\partial v^*}{\partial x^*}+v^*\frac {\partial v^*}{\partial y^*}&=-\frac {\partial p^*}{\partial y^*}+{\textit{Re}}^{-1}\bigg (\frac {\partial ^2 v^*}{\partial x^{*2}}+\frac {\partial ^2 v^*}{\partial y^{*2}}\bigg ), \end{align}
\begin{align} \frac {\partial u^*}{\partial x^*}+\frac {\partial v^*}{\partial y^*} & = 0, \end{align}
where
${\textit{Re}}$
is the Reynolds number. We assume that
$\theta ^*$
satisfies the dimensionless steady heat transfer equation
\begin{equation} u^*\frac {\partial \theta ^*}{\partial x^*}+v^*\frac {\partial \theta ^*}{\partial y^*} = (\sigma {\textit{Re}})^{-1}\bigg (\frac {\partial ^2 \theta ^*}{\partial x^{*2}}+\frac {\partial ^2 \theta ^*}{\partial y^{*2}}\bigg ), \end{equation}
where
$\sigma$
is the Prandtl number. As discussed shortly, we assume that the freezing time scale is significantly larger than that of the fluid flow, which, to repeat, allows motion to be treated as quasi-steady as per (2.1) and (2.2).
A schematic illustrating the model geometry, where
$y=f(x,\,t)$
denotes the position of the constricting ice boundary and
$x=\delta$
represents the point at which the wall temperature drops from
$\theta =0$
to
$\theta =-1,$
as per (2.5a
).
In practical terms, taking a characteristic velocity of
$U^* = 1 \,\mathrm{m\,s^{-1}}$
, which is typical for household pipes (Żabnieńska–Góra & Dudkiewicz Reference Żabnieńska–Góra and Dudkiewicz2018), a length scale
$L^* = 10^{-3} {-} 10^{-2}\,\mathrm{m}$
and a representative kinematic viscosity
$\nu ^* = 10^{-6}\,\mathrm{m^2\,s^{-1}}$
for water of
$0^\circ$
C (Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot2002) gives an approximate Reynolds number
${\textit{Re}} = U^* L^*/\nu ^* \approx 10^3 {-} 10^4$
. A typical value of the Prandtl number for water is
$\sigma =8$
, and we henceforth focus on the asymptotically high-Reynolds-number case, whilst taking
$\sigma =\mathcal{O}(1)$
. The time scale of the water flow,
$L^*/U^*,$
can feasibly vary in the range of
$ 10^{-3}{-}10^{-2}\, \text{s}$
in view of the values mentioned previously. The time scale of phase change which depends on the Stefan number (defined as
$\mathrm{St} = c_p\Delta \theta /L_f$
, where
$c_p, \Delta \theta$
and
$L_f$
are the specific fluid heat, characteristic temperature difference and latent heat) varies widely. Following Jolley et al. (Reference Jolley, Dang and Smith2024), we thus take the Stefan number to be such that the phase-change time scale is significantly larger than that of the fluid flow. Lastly, we follow Moore et al. (Reference Moore, Mughal and Papageorgiou2017) and Batley & Smith (Reference Batley and Smith2025) in neglecting the effects of viscous dissipation, which is justifiable given that the Eckert number is given by
${Ec} = U^{*2}/c_p\Delta \theta \approx 2.38 \times 10^{-4}$
with
$c_p = 4.2\times 10^3 \,\mathrm{J\,kg^{-1}\, K^{-1}}$
and
$\Delta \theta = 1\, \mathrm{K}.$
As discussed previously, we pay particular interest to case
${\textit{Re}}\gg 1,$
so that the flow is driven by an inertial-viscous balance. In view of this, we introduce the variables
$x^*={\textit{Re}}x, y^*=y$
and expansions
\begin{equation} [u^*,\,v^*,\,p^*,\,\theta ^*] \sim [u,\,{\textit{Re}}^{-1}v,\,p,\,\theta ](X,\,Y,\,t), \end{equation}
so that (2.1) and (2.2) provide
\begin{align} uu_x+vu_y &= -p_x + u_{yy} , \end{align}
\begin{align} u_x + v_y &= 0, \end{align}
\begin{align} u\theta _x+v\theta _y &= \sigma ^{-1}\theta _{yy}, \\[12pt] \nonumber \end{align}
to leading order in
${\textit{Re}}^{-1}\ll 1.$
We also have
$p:=p(x,\,t)$
, as the lateral momentum balance reduces to
$\partial p/\partial y=0$
to leading order in
${\textit{Re}}^{-1}.$
2.2. Boundary and initial conditions
We now present the boundary conditions which couple the governing equations from (2.4). The ice–fluid interface is fixed at the freezing temperature of
$\theta =-1$
and we specify a sudden drop in temperature from
$\theta =0$
to
$\theta =-1$
on the upper wall at
$x=\delta ,$
so that
\begin{equation} \theta \Big \rvert _{y=f} = \begin{cases} -\mathrm{H}(x-\delta ) \quad\,\, \text{if} \quad f=1,\\ -1 \qquad \,\qquad \,\, \text{if} \quad f\lt 1, \end{cases} \end{equation}
where
$\delta \gt 0$
is an
$\mathcal{O}(1)$
constant and
${H}$
is the Heaviside function. On the ice boundary and channel wall, we impose the no-slip and no-penetration conditions
\begin{equation} u=v=0 \quad \text{at} \, y=f. \end{equation}
We impose symmetry conditions at
$y=0$
so that
\begin{equation} \frac {\partial u}{\partial y}=v=\frac {\partial \theta }{\partial y}=0 \quad \text{at} \, y=0. \end{equation}
We also prescribe the steady inlet conditions
\begin{align} u&=u_0(y) \quad \text{at} \, x=0, \end{align}
\begin{align} \theta &=\theta _0(y) \quad\, \text{at} \, x=0. \end{align}
Finally, we prescribe the Stefan condition which determines the position of the ice boundary
\begin{align} \frac {\partial f}{\partial t}&=-\theta _y \quad \text{at} \, y=f. \end{align}
Given that the channel wall is rigid and permanent (as with metal channels for instance), this Stefan condition only applies for all
$x$
beyond that where
$f\lt 1,$
the formation of ice not occurring when
$f\geq 1$
is obtained. As such, we have
$f=1$
when (2.5f
) provides
$f\gt 1$
. We note that, given
$\theta =-1$
on both the channel wall and ice–fluid interface
$y=f$
, there is no thermal flux term contributed by the ice temperature in (2.5f
). For simplicity, we assume that there is no initial ice profile at
$t=0,$
so that
$f(x,\,0)=1$
.
3. Numerical results
To illustrate the behaviour of the partial differential equation (PDE) system from (2.4) and (2.5), we now present and discuss some exemplar numerical solutions. We also present parameter regions in which ice formation does and does not occur. We assume parabolic inlet conditions of the form
\begin{equation} (u_0,\,\theta _0) = (\lambda ,\,\phi ) (1-y^2), \end{equation}
where
$\lambda \gt 0$
and
$\phi \lt 0$
are constant. For simplicity, we fix
$\delta =0.1$
. The methods used to obtain numerical solutions are presented in Appendix A.
In figure 2, we present
$\theta ,\,u$
and
$v$
at
$t= 0.18$
, as well as
$f$
and
$p$
at fixed times in
$t\in [0,\,0.72],$
for
$\phi =-10/4$
and
$\lambda =15/4$
. Figure 2(a) illustrates the formation and growth of an ice prong, which results from the wall temperature reducing to the freezing temperature at
$x=0.1$
. The leading edge of the growing ice prong, located at
$x\approx 0.11$
for all
$t$
, is intriguingly positioned downstream of the wall temperature drop at
$x=0.1$
. This occurs as the temperature near the channel walls remains above
$ \theta = -1$
for
$ 0.1\leq x \leq 0.11$
(see figure 2
f); consequently, the right-hand side of the Stefan condition from (2.5f
) within this
$x$
-range is positive (leading to
$f\gt 1$
), thereby inhibiting ice growth. The asymptotic behaviour of
$\theta$
and
$f$
close to the wall temperature drop is investigated in Appendix B.
Numerical solutions of the model from (2.4) and (2.5) with
$\sigma =1,$
using the inlet conditions from (3.1) where
$\phi =-10/4$
and
$\lambda =15/4$
. Panels (b,c) and (d) show solutions for
${\theta },\,u$
and
$v$
at
$t=0.72$
, respectively. Panels (a) and (e) illustrate
$f$
and
$p$
for uniformly distributed, fixed values of
$t\in [0,\,0.72],$
where the arrows point in the direction of increasing
$t.$
Panel (f) shows
$\theta$
at
$t=0.4$
for uniformly distributed, fixed values of
$x\in [0.09,\,0.115]$
.
Figure 2 shows the ice prong growing inwards with
$t$
towards
$y=0$
(the symmetry line of a full channel) and forming a narrow gap about
$x \approx 0.16$
, where
$u,\,v$
and
$p$
are significantly larger than elsewhere in the channel. The fluid temperature in this gap is a small negative perturbation of the freezing temperature, and the gap width about
$x\approx 0.16$
approaches zero in a finite time. The features described here are derived analytically and discussed further in § 4, where the asymptotic behaviour of the model solutions are examined as the gap width between
$y=0$
and
$y=f$
approaches zero. As the fluid passes through the narrow gap and travels downstream, it is gradually warmed towards the freezing temperature,
$\theta =-1,$
by the ice boundaries, so that
$\theta \rightarrow -1^{-}$
and hence
$\theta _y\rightarrow 0$
as
$x\rightarrow \infty .$
As such, the ice prong width reduces as
$x\rightarrow \infty$
, as suggested by figure 2.
In contrast to the results shown in figure 2, numerical simulations suggest that certain parameter regimes do not give rise to ice formation, i.e.
$ f = 1$
for all
$t$
. To illustrate this, figure 3 displays regions in
$(\lambda ,\,\sigma )$
-space in which ice formation either does or does not occur, for four different values of
$\phi$
, obtained via numerical solutions of (2.4) and (2.5). The value of
$\sigma$
and
$\lambda$
are seen to play an important role in determining the absence or formation of ice. Specifically, for larger
$\sigma ,$
the minimal value of
$\lambda$
required for ice formation decreases. This aligns with expectation, as increasing
$\sigma$
corresponds to either increasing or decreasing the momentum or thermal diffusivity, respectively, either of which inhibits the transfer of heat from the warm-wall region, thereby promoting ice formation. Furthermore, the minimal value of
$\lambda$
required for ice formation increases as
$|\phi |$
decreases. This suggests that the warmer the supercooled inlet fluid, the faster it must be transported to the cool-wall region
$(0.1\leq x)$
, thereby minimising its exposure to warming in the warm-wall region
$(x\lt 0.1).$
Regions in
$(\lambda ,\,\sigma )$
-space in which ice formation is or is not expected, for
$|\phi |=6,\,10,\,14,\,18.$
The solid black lines which separate the two regions are obtained numerically from (2.4) and (2.5). The solid black arrows point to the regions in which ice is or is not expected and the dashed arrow points in the increasing direction of
$|\phi |$
.
As seen in figure 3, the demarcation curves separating regions of
$(\lambda ,\,\sigma )$
parameter space in which ice does or does not form are straight lines in log–log space with gradient of negative unity, and hence take the power-law form
$\sigma = C(\phi )/\lambda$
for a fixed value of
$\phi .$
Here, a large range of
$\sigma$
is shown for completeness; however, as noted in § 2,
$\sigma \approx 8$
is appropriate for water. The power-law relationship can, in fact, be deduced analytically: in the absence of ice
$(f=1),$
the boundary layer equations from (2.4) admit the analytical solutions
$v=0$
and
$u=\lambda (1-y^2)$
, so that (2.4c
) reveals that the thermal dynamics depends solely on the compound parameter
$\lambda \sigma .$
For a fixed value of
$\phi$
, we expect that there is a unique value of this compound parameter such that the right-hand side of the Stefan condition from (2.5f
) is zero, thus identifying the demarcation curve. Direct numerical simulations of (2.4) and (2.5) strongly suggest that when ice forms, the gap width between
$y=0$
and
$y=f$
always approaches zero at some finite value of
$t.$
In this case, the regions in figure 3 also correspond to regions where the channel does or does not close in finite time. A more detailed analysis of the regions in which channel closure does or does not occur is undertaken in the following section, using an alternative approach.
4. Asymptotic behaviour at the onset of channel closure
We now examine the asymptotic behaviour of the model from (2.4) and (2.5) at the onset of channel closure due to the constricting ice boundary, as seen in figure 2. Interestingly, this analysis provides insight into when channel closure due to ice growth is or is not expected. In the following analysis, we exploit the small quantities
$0\lt t_c-t\ll 1$
and
$|x-x_c|\ll 1$
, where
$t_c$
and
$x_c$
represent the time and
$x$
-coordinate at which the ice-prong impacts
$y=0,$
so that
$f(x_c,\,t_c)=0.$
Via a dominant balance argument, the transformation and expansions of the variables which provide a full balance when
$0\lt t_c-t\ll 1$
and
$|x-x_c|\ll 1$
are
\begin{equation} (x-x_c,y) = (t_c-t)^{\frac {m+1}{2}}(X,Y) \end{equation}
and
\begin{align} (u,v) \sim (t_c-t)^{-\frac {m+1}{2}}&\big [U(X,\,Y),V(X,\,Y)\big ],\, \quad p \sim (t_c-t)^{-m-1}P(X), \nonumber\\&\qquad f \sim (t_c-t)^{\frac {m+1}{2}}F(X), \end{align}
The Stefan condition from (2.5f ) suggests that the temperature expands as
\begin{equation} \theta \sim -1 + (t_c-t)^m\varTheta (X,\,Y). \end{equation}
The value of
$m\gt 0$
is obtained in § 4.1 and its inclusion here is required to allow a thermal variation when
$-X\gg 1$
. For instance, whilst taking
$m=0$
provides the same full balance as with
$m\gt 0,$
it necessitates that
$\theta = -1$
when
$-X\gg 1$
; this uniform temperature then forces
$F$
to remain linear in the scaled streamwise coordinate
$X$
, a feature which is only associated with
$F$
away from
$x=x_c$
. This constraint then motivates the above-mentioned scalings, which consider an asymptotically smaller spatial region than with
$m=0.$
As seen shortly in § 4.1, the quantity
$m$
can also indicate whether channel closure due to ice growth occurs.
In view of the above-mentioned scalings and expansions from (4.1b ), the governing equations from (2.4) provide the full balance
\begin{align} U_X+V_Y&=0, \end{align}
\begin{align} UU_X+VU_Y &= -P_X+U_{YY}, \end{align}
\begin{align} U\varTheta _X+V\varTheta _Y &= \sigma ^{-1}\varTheta _{YY}, \end{align}
to leading order in
$0\lt t_c-t\ll 1,$
where
$-\infty \lt X\lt \infty$
and
$0\leq Y\leq F(X)$
. As first discussed in § 2.1, we assume that the time derivative terms absent from (2.1) and (2.2) remain sufficiently small such that they can again be neglected in (4.2). We note that, however, new balances will eventually arise when unsteady contributions overtake the quasi-steady solutions, given the large velocity expansions in (4.1b
), for example. Using (2.5), the system from (4.2) is subject to the boundary conditions
\begin{align} U&=V=0 \quad \text{and} \quad \varTheta =0 \quad \,\,\,\, \text{at} \, Y=F(X), \end{align}
\begin{align} \frac {\partial U}{\partial Y}&=V=0 \quad \text{and} \quad \frac {\partial \varTheta }{\partial Y}=0 \quad \text{at} \, Y=0, \end{align}
as well as the Stefan condition
\begin{equation} (m+1)\big (F-XF_X\big )=2\varTheta _Y \quad \text{at} \, Y=F(X). \end{equation}
The appropriate initial conditions which are imposed at some sufficiently large value of
$-X$
, as well as the value of
$m,$
are obtained via an analysis of (4.2) as
$X\rightarrow -\infty$
, as follows in § 4.1.
4.1. Asymptotic analysis as
$X\rightarrow -\infty$
To analyse the system from (4.2) as
$X\rightarrow -\infty$
, it is first instructive to introduce the transformations
\begin{equation} Y=\zeta F, \quad U=\frac {\bar {U}}{F}, \quad V=\bar {V}+\frac {\zeta F_X\bar {U}}{F}, \end{equation}
which fixes the new lateral coordinate onto the unit domain as
$0\leq \zeta \leq 1.$
In view of this transformation, (4.2) provides
\begin{align} \bar {U}_X+\bar {V}_\zeta &=0, \end{align}
\begin{align} \bar {U}\bar {U}_X+\bar {V}\bar {U}_\zeta -\frac {F_X \bar {U}^2}{F} &= -F^2P_X+\frac {\bar {U}_{\zeta \zeta }}{F}, \end{align}
\begin{align} \bar {U}\varTheta _X+\bar {V}\varTheta _\zeta &= \sigma ^{-1}\frac {\varTheta _{\zeta \zeta }}{F}, \end{align}
subject to the boundary conditions
\begin{align} \bar {U}&=\bar {V}=0 \quad \text{and} \quad \varTheta =0 \quad \quad \text{at} \, \zeta =1, \end{align}
\begin{align} \frac {\partial \bar {U}}{\partial \zeta }&=\bar {V}=0 \quad \text{and} \quad \frac {\partial \varTheta }{\partial \zeta }=0 \quad\, \text{at} \, \zeta =0, \end{align}
and the Stefan condition
\begin{equation} (m+1)F(F-XF_X)=2\varTheta _\zeta \quad \text{at} \, \zeta =1. \end{equation}
This system is numerically integrated forwards in
$X$
, starting from a sufficiently large value of
$-X$
such that the solution converges. We use the
$-X\gg 1$
expansions of (4.4) presented shortly in (4.5) as initial conditions; the value of
$m$
is also obtained via this asymptotic analysis. Numerical solutions and quantities arising from (4.4) are presented and discussed in Appendix C.
The appropriate expansions when
$-X\gg 1$
for the new scaled variables are obtained by requiring that the original variables contribute an
$\mathcal{O}(1)$
effect when
$x=x_c-\mathcal{O}(1)$
. Writing
$-X=\xi ,$
for
$\xi \gg 1$
, we have
\begin{equation} \bar {U} \sim F_0R(\zeta ), \qquad \bar {V}\rightarrow 0, \qquad P \sim Q\xi ^{-2}, \qquad \varTheta \sim \xi ^{b+1}T(\zeta ), \qquad F \sim F_0\xi +a\xi ^b, \end{equation}
where
$Q\lt 0$
,
$a$
is an arbitrary amplitude and
$F_0\gt 0$
describes the ice shape far upstream. The inclusion of the unknown constant
$b$
is required to prevent
$T=0$
when
$\xi \gg 1$
. Furthermore, the presence of the power
$b+1$
and
$b$
in the expansions for
$\varTheta$
and
$F,$
respectively, is required to obtain a non-trivial balance in the Stefan condition from (4.4f
). To relate
$b$
to
$m$
, we require the correction term from (4.1c
) when
$\xi \gg 1$
,
\begin{equation} (t_c-t)^m\Big (-(x-x_c)(t_c-t)^{-\frac {m+1}{2}}\Big )^{b+1}T(\zeta ), \end{equation}
to contribute an
$\mathcal{O}(1)$
effect when
$x=x_c-\mathcal{O}(1)$
, so that
\begin{equation} m = \frac {1+b}{1-b}. \end{equation}
In view of (4.5), the system from (4.4) provides the ordinary differential equations
\begin{align} R''&=F_0^2(2Q+R^2), \end{align}
\begin{align} T''&=-F_0^2\sigma (b+1)RT, \end{align}
to leading order when
$\xi \gg 1,$
with the corresponding boundary conditions
\begin{align} R&=T=0 \quad \, \text{at} \, \zeta =1, \end{align}
\begin{align} R'&=T'=0\quad \text{at} \, \zeta =0, \end{align}
\begin{align} T'&=F_0a \vphantom{=0} \quad\quad\kern2pt \text{at} \, \zeta =1. \end{align}
We note that the last condition here has been simplified by using the value of
$m$
provided in (4.6b
). Using the first from (4.5), conservation of mass also requires that
\begin{equation} \int _0^1 R\,\mathrm{d}\zeta =\frac {\mathcal{C}}{F_0}, \end{equation}
where
$\mathcal{C}$
is the known mass flux; see Jeffery (Reference Jeffery1915) and Hamel (Reference Hamel1917). Given that the shape parameter
$F_0$
and mass flux
$\mathcal{C}$
are known and determined from the flow and thermal behaviour when
$x=x_c-\mathcal{O}(1)$
,
$R$
and
$Q$
are determined first via (4.7a
), (4.7f
) and attendant boundary conditions. Subsequently, the eigenvalue
$b$
and function
$T$
are determined via (4.7b
)–(4.7e
). Given that
$R$
is non-negative, the sign on the right-hand-side of (4.7b
) indicates that
$b+1\gt 0$
to allow a solution. In addition,
$b$
must satisfy
$b\lt 1$
for the expansion in
$F$
from (4.5) to be asymptotically valid, so that
$-1\lt b\lt 1$
. This indicates that
$0\lt m\lt \infty$
, thus confirming via (4.1b
) the finite-time blow up of solutions as
$t\rightarrow t_c$
; possible regularising effects are proposed in § 5. Nevertheless, the value of
$b$
provides an interesting indication of whether channel closure due to ice growth occurs, as is now discussed.
In the following analysis, we restrict our attention to the cases where (4.7) admits zero or one values of
$-1\lt b\lt 1$
. The transition between these cases represents the bifurcation between the non-existence and existence of solutions to the system from (4.4), in view of the
$-X\gg 1$
behaviour from (4.5). Given that, in practice, (4.7) can admit values of
$b$
larger than unity (although such values are not physically relevant), the transition between zero and one admissible values of
$b$
is therefore captured by the smallest admissible value of
$b,$
say
$ b=:b_{\text{s}}$
, as it crosses the demarcation curve
$ b_{\text{s}} = 1$
. Physically, the existence of solutions
$(b_{\text{s}}\lt 1)$
to (4.4) corresponds to channel closure due to ice growth. We now obtain the value of
$b_{\text{s}}$
for some extreme flow rates.
4.2. Analysis of the eigenvalue problem from (4.7) for extreme flow rates
Here, we exploit the two extreme cases of
$R$
being relatively small and large to obtain the value of
$Q$
and eigenvalue
$b_{\text{s}}$
explicitly from (4.7), in terms of the shape and mass flux parameters
$F_0$
and
$\mathcal{C}.$
We then compare these results with those obtained numerically from the full eigenvalue problem from (4.7).
4.2.1. The lubrication limit,
$\mathcal{C}F_0\ll 1$
We first address the case of
$R$
being relatively small, in which (4.7a
) suggests that a lubrication-type balance between the first and second terms (i.e. a viscous-pressure balance) occurs when
$\mathcal{C}F_0\ll 1$
. This condition is the necessary requirement to neglect the nonlinear term. As such, (4.7a
) and (4.7f
) provide
\begin{equation} R \sim \frac {3\mathcal{C}}{2F_0}(1-\zeta ^2), \qquad Q \sim -\frac {3\mathcal{C}}{2F_0^3} , \end{equation}
to leading order in
$\mathcal{C}F_0\ll 1$
. In view of this, solving (4.7b
) with
$T'(0)=0$
gives
\begin{equation} T(\zeta ;\, \beta ) \sim \frac {C_1}{\sqrt {\zeta }}\bigg [ {\mathrm{W}}_{\frac {\sqrt {\beta }}{4},\frac {1}{4}}\left (\sqrt {\beta }\,\zeta ^2\right )-\gamma {\mathrm{M}}_{\frac {\sqrt {\beta }}{4},\frac {1}{4}}\left (\sqrt {\beta }\,\zeta ^2\right ) \bigg ], \end{equation}
such that
\begin{equation} \gamma = \frac {\sqrt {\pi }\,\,\left [4\,\varGamma \left (\frac {1}{4}-\frac {\sqrt {\beta }}{4}\right )+\varGamma \left (-\frac {\sqrt {\beta }}{4}-\frac {3}{4}\right )+\sqrt {\beta }\,\varGamma \left (-\frac {\sqrt {\beta }}{4}-\frac {3}{4}\right )\right ]}{\varGamma \left (\frac {1}{4}-\frac {\sqrt {\beta }}{4}\right )\varGamma \left (-\frac {\sqrt {\beta }}{4}-\frac {3}{4}\right )}, \end{equation}
where
${W}$
and
${M}$
are the Whittaker-
${W}$
and -
${M}$
functions,
$\varGamma$
is the Gamma function,
$C_1$
is an arbitrary constant and
$\beta (b_{\text{s}}) = 3\mathcal{C}F_0\sigma (b_{\text{s}}+1)/2$
is an eigenvalue.
The value of
$\beta$
which leads to
$b_{\text{s}}$
is
$\beta =\beta _{\text{s}} :\approx 2.829$
, so that
\begin{equation} b_{\text{s}} \sim \frac {2\beta _{\text{s}}}{3\mathcal{C}F_0\sigma }-1. \end{equation}
Given that
$\beta _{\text{s}}$
is known,
$C_1$
can now in theory be determined numerically for a given value of
$a$
via (4.7e
). In view of (4.10), the condition
$b_{\text{s}}\lt 1$
which represents channel closure due to ice growth is only realisable when
$\mathcal{C}F_0\ll 1$
if
$\sigma = \mathcal{O}(1/\mathcal{C}F_0)$
i.e. a large Prandtl number is required for closure in the presence of lubrication-type flow. In this instance, the demarcation curve
$b_{\text{s}}=1$
, which separates regions in which ice closure is expected
$(b_{\text{s}}\lt 1)$
and not expected
$(b_{\text{s}}\gt 1)$
, is hence parametrised by
\begin{equation} \sigma \sim \frac {\beta _{\text{s}}}{3\mathcal{C}F_0}. \end{equation}
A discussion of this result as well as a comparison with the full numerical solutions of (4.7) is presented in § 3.
4.2.2. The plug-flow limit,
$\mathcal{C}F_0\gg 1$
In the case of
$R$
being relatively large, (4.7a
) suggests that a plug flow-type balance between the second and third terms (i.e. an inertial-pressure balance) occurs when
$|Q|\gg 1$
, which is valid away from a boundary layer at
$\zeta =1,$
so that
$R\sim \sqrt {2|Q|}.$
The effects of viscosity on the flow are non-negligible in the boundary layer, in which the scalings
$R \sim \sqrt {2|Q|}\bar {R}(\eta )$
and
$\zeta =1-\kappa \eta$
hold, where
$\kappa = 2^{1/4}F_0^{-1}|Q|^{-1/4}\ll 1$
and
$\eta ,\,\bar {R}=\mathcal{O}(1)$
. As such, (4.7a
) provides
\begin{equation} \bar {R}''=2(\bar {R}^2-1), \end{equation}
where
$\bar {R}(0)=0$
and
$\lim _{\eta \rightarrow \infty }\bar {R}'=0,$
the solution of which is (Goldstein Reference Goldstein1938)
\begin{equation} \bar {R} = -2+3\tanh ^2(\eta +\ell ), \quad \ell =\text{arctanh}\bigg (\sqrt {\frac {2}{3}}\bigg ). \end{equation}
This solution facilitates the automatic matching of the outer and boundary layer expansions since
$\lim _{\eta \rightarrow \infty }\bar {R}=1,$
so that the boundary layer expansion represents a fully composite solution for all
$\zeta .$
In view of this, (4.7f
) provides the transcendental equation for
$Q$
,
\begin{equation} \sqrt {2|Q|}\Big [1+3\kappa \ell -3\kappa \tanh (\ell +\kappa ^{-1})\Big ] \sim \frac {\mathcal{C}}{F_0}. \end{equation}
Given that
$\kappa ^{-1}\gg 1,$
we have that
$\tanh (\ell +\kappa ^{-1})\sim 1$
(neglecting only exponentially small terms here) and hence finally
\begin{equation} Q \sim -2\Bigg (\frac {\mathcal{C}}{3(1-\ell )+\sqrt {2F_0\mathcal{C}+9(1-\ell )^2}}\Bigg )^{4}. \end{equation}
This result, together with (4.7f
) and the above-mentioned outer solution, imply that
$\mathcal{C}F_0\gg 1$
is a sufficient condition to neglect the first term from (4.7a
).
As just discussed, the viscous effects on the flow are non-negligible in the boundary layer at
$\zeta =1$
; however, this layer is passive in the thermal response as
$T$
from (4.7b
) can satisfy
$T(1)=T'(0)=0$
using the solution away from the layer, namely
$R\sim \sqrt {2|Q|}$
. In view of this, solving (4.7b
) and using
$T'(0)=0$
provides
\begin{equation} T(\zeta ;\,\lambda ) \sim C_2 \cos (\sqrt {\lambda }\zeta ), \end{equation}
to leading order, where
$C_2$
is an arbitrary constant which can be determined by (4.7e
) and
$\lambda (b_{\text{s}}) = \sqrt {2|Q|}F_0^2\sigma (b_{\text{s}}+1).$
The value of
$\lambda$
which satisfies
$T(1)=0$
and gives rise to
$b_{\text{s}}$
is
$\lambda =\pi ^2/4,$
so that
\begin{equation} b_{\text{s}} \sim \frac {\pi ^2}{4\sqrt {2|Q|}F_0^2\sigma }-1. \end{equation}
Given that
$\sqrt {|Q|}\gg 1,$
we therefore require
$\sigma = \mathcal{O}(|Q|^{-1/2})$
to allow
$b_s=1$
to be realisable. In this instance, the demarcation curve
$b_{\text{s}}=1$
is parametrised by
\begin{equation} \sigma \sim \frac {\pi ^2}{8\sqrt {2|Q|}F_0^2}. \end{equation}
To draw comparison between the demarcation curves obtained via the lubrication and plug-flow limit approximations from (4.11) and (4.16), we note that as
$CF_0\rightarrow \infty$
, we have from (4.13b
) that
$Q\rightarrow -\mathcal{C}^2/(2F_0^2)$
, so that
$b_{\text{s}}\rightarrow \pi ^2/(4\mathcal{C}F_0\sigma )-1.$
As such, the demarcation curve as
$CF_0\rightarrow \infty$
is parametrised by
$\sigma =\pi ^2/(8\mathcal{C}F_0).$
Remarkably, the demarcation curve
$b_s=1$
from the lubrication
$(\mathcal{C}F_0\ll 1)$
and plug-flow
$(\mathcal{C}F_0\gg 1)$
limits are thus described by a similar power law, only differing by a multiplicative constant. This finding will be illustrated in the following subsection.
4.2.3. Comparison between asymptotic and numerical results
We now compare the asymptotic results obtained for the two above-mentioned extreme flow rates with numerical solutions of the eigenvalue problem from (4.7). First,
$R$
and
$Q$
are determined via (4.7a
), (4.7f
) and attendant boundary conditions. Then, by fixing
$T'(0)=0$
, the value of
$b$
is selected by solving for
$T$
using a shooting method and requiring that
$T(1)=0$
is satisfied for an arbitrary value of
$T(0)$
. This value can then be subsequently determined via (4.7e
) for any prescribed value of
$a,$
thereby fully determining
$T.$
Solid black lines represent the quantities (a)
$|b_{\text{s}}|$
and (b)
$Q$
obtained numerically from (4.7) for
$\sigma =1$
and
$F_0=0.1,\,1,\,10$
. The dashed green line represents the lubrication limit approximation of
$Q$
and
$b_{\text{s}}$
from (4.8) and (4.10), whereas the the dotted pink line represents the plug-flow limit approximation of
$Q$
and
$b_{\text{s}}$
from (4.13b
) and (4.15). The arrow points in the increasing direction of
$F_0$
.
Demarcation curve
$b_{\text{s}}=1$
for
$F_0=0.1,\,1,\,10$
which separates regions in which channel closure due to ice growth is expected
$(b_{\text{s}}\gt 1)$
and not expected
$(b_{\text{s}}\lt 1)$
. The black curve represents
$b_{\text{s}}=1$
obtained numerically from (4.7). The dashed green and dotted pink line represent the lubrication limit approximation from (4.11) and the plug-flow limit approximation from (4.16). For each demarcation curve, the upward and downward pointing arrows indicate regions in parameter space such that
$b_{\text{s}}\gt 1$
and
$b_{\text{s}}\lt 1,$
respectively.
In figure 4, we illustrate the quantities
$|b_{\text{s}}|$
(panel a) and
$Q$
(panel b) as functions of
$\mathcal{C}$
for various values of
$F_0$
, obtained both numerically and via the above-mentioned asymptotic analyses. Strikingly, there is an excellent agreement between the asymptotic and numerical solutions. Of perhaps more interest, the demarcation curve
$b_{\text{s}}=1,$
plotted
$(\mathcal{C},\,\sigma )$
-space for various values of
$F_0$
, is presented in figure 5, again illustrating an excellent agreement between the numerical and asymptotic approximations. As demonstrated in § 4.2.2, approximations of the demarcation curve via the lubrication and plug-flow limit both relate
$\sigma$
proportionally to
$\mathcal{C}^{-1},$
with the proportionality constant differing only by approximately one half. Indeed, this close similarity is reflected well in figure 5, where the inset highlights the minor transition between the two asymptotic approximations. In this specific example (and others not shown here), each
$b_{\text{s}}=1$
curve can be seen to transition monotonically between the two asymptotic approximations, which thus tightly bound the demarcation curve for all values of
$\mathcal{C}.$
5. Discussion
In this work, we develop and analyse a model which describes the interaction between high-Reynolds-number channel flow and wall-mounted ice that evolves over a long axial length scale. The fluid flow and temperature dynamics in the channel are governed by the boundary layer and heat transfer equations, which are treated as steady due to the comparatively slowly evolving ice boundary (Moore et al. Reference Moore, Mughal and Papageorgiou2017; Jolley et al. Reference Jolley, Dang and Smith2024). The location of this ice boundary (i.e. ice–water interface), which is maintained at the ice melt temperature, is related to the thermal flux of the fluid there via a classical Stefan condition (Fowler, Reference Fowler1997). The channel walls on which the ice forms are also fixed at the constant melt temperature and, together with the aforementioned long axial length scale assumption, the ice temperature is approximated by this constant. This justifies the absence of the ice temperature flux in the Stefan condition.
Direct numerical simulations of the model described here demonstrate that ice can form on the channel walls and subsequently grow inwards towards the channel centreline, eventually resulting in significant flow constriction there. These simulations also strongly suggested that when ice forms, then constriction always occurs (i.e. there are no steady-state solutions). The emergence of these constricting ice boundaries suggest that the channel gap width may reduce to zero in scaled finite time
$(t=t_c)$
and at a fixed point longitudinally. This observation motivates an asymptotic analysis of the model local to the critical point, wherein the model variables and scaled and expanded in constant powers of
$t_c-t\ll 1.$
The associated exponents were all found to be related to an eigenvalue,
$b,$
which must in general be calculated numerically from a system of nonlinear ordinary differential equations which describe the fluid and thermal dynamics far upstream from the critical point. The asymptotic analysis confirms the finite-time blow-up of solutions as the channel gap width reduces to zero, although such behaviour is anticipated in the absence of regularising mechanisms such as viscous dissipation or pressure-induced melting (regelation). Viscous dissipation in particular will become increasingly significant as the ice constriction narrows, due to the corresponding increase in fluid velocity there.
Conversely, direct numerical simulations indicate that in certain regions of parameter space, ice does not form on the channel walls. This observation motivated further analysis to identify the conditions under which ice growth is expected. First, analysis of the full model reveals that ice formation occurs when the Prandtl number,
$\sigma$
, satisfies
$\sigma \gt C/\lambda$
, where
$\lambda$
scales linearly with the prescribed longitudinal inlet fluid velocity and
$C$
is a constant that depends on the prescribed inlet temperature. Second, the asymptotic analysis near the critical point (as previously described) offers supplementary insight into ice growth parameter regimes. In particular, when the eigenvalue
$b$
(which can be expressed in terms of the model parameters) satisfies
$b\gt 1$
, the asymptotic expansions far upstream of the critical point cease to be valid, indicating that channel constriction due to ice growth (and thus ice formation) does not occur.
As mentioned, a natural extension of the present work would be to incorporate additional physical effects, such as viscous dissipation and pressure-induced melting (regelation), which may regularise the model by preventing the channel gap width from reducing to zero as
$t\rightarrow t_c$
. Another possible extension would be to consider an axisymmetric cylindrical geometry in place of the parallel channel configuration, thereby improving the model’s relevance to wall-mounted ice formation in pipe or tube flows, as encountered in plumbing applications. Furthermore, in applications related to mechanical systems, it may be more appropriate to consider a fixed fluid flow power in place of prescribing an inlet condition, as was done in this paper.
Acknowledgements
Thanks are due to UCL for the support of R.B. and to the anonymous reviewers for helpful suggestions.
Funding
We thank the Leverhulme Trust; the work was supported by a Leverhulme Trust Research Project Grant RPG-2023-233 supporting J.J.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical methods for the full problem from (2.4) and (2.5)
We now discuss the numerical methods used to solve the governing system of equations from (2.4) and (2.5). Similarly to the Prandtl transposition, we introduce straight-channel geometry via the new coordinate system
$[\bar {X}(x,y),\bar {Y}(x,y)]$
, where
$\hat {X}=x$
and
$\hat {Y}=y/f$
such that
$\bar {Y}(x,f(x)) = 1.$
By defining
$\hat {\theta }=\theta$
,
$\hat {p}=p,\,{ \hat {u}=u}$
and
\begin{equation} {\hat {v} = \frac {v-\hat {Y}\hat {u}f_{\hat {X}}}{f}}, \end{equation}
the boundary layer and thermal equations from (2.5) provide
\begin{align} \hat {u}_{\hat {X}} + \hat {v}_{\hat {Y}} &= -\frac {f_{\hat {X}}}{f}\hat {u}, \end{align}
\begin{align} \hat {u}\hat {u}_{\hat {X}}+\hat {v}\hat {u}_{\hat {Y}} &= -\hat {p}_{\hat {X}}+\frac {1}{f^2}\hat {u}_{\hat {Y}\hat {Y}}, \end{align}
\begin{align} \hat {u}\hat {\theta }_{\hat {X}}+\hat {v}\hat {\theta }_{\hat {Y}} &= \frac {1}{f^2}\hat {\theta }_{\hat {Y}\hat {Y}}. \end{align}
The attendant boundary condition is simply given by (2.5) but with
$x$
replaced with
$\hat {X}$
and so forth, and
$\hat {Y}=1$
replacing variable evaluations at
$y=f.$
We introduce the discretisations
$\hat {X}_i=i\Delta x$
and
$\hat {Y}_j=j\Delta y$
, such that
$i=1,2,\ldots ,M$
and
$j=1,2,\ldots ,N$
, and
$\Delta x=\max \{\hat {X}\}/(M-1)$
and
$\Delta y=1/(N-1)$
are the step sizes. We denote
$\hat {U}(\hat {X}_i,\,\hat {Y}_j) \approx u_{i,\,j}$
etc., and implement the following finite difference scheme to solve (A2):
\begin{align} \frac {u_{i,j}-u_{i-1,j}}{\Delta x}+\frac {v_{i-1,j+1}-v_{i-1,j-1}}{2\Delta y} & = -\frac {f_{i}-f_{i-1}}{f_i\Delta x}u_{i-1,j}, \\[-28pt] \nonumber \end{align}
\begin{align} u_{i-1,j}\frac {u_{i,j}-u_{i-1,j}}{\Delta x} & + v_{i-1,j}\frac {u_{i-1,j+1}-u_{i-1,j-1}}{2\Delta y} = -\frac {p_i-p_{i-1}}{\Delta x} \nonumber\\& + \frac {1}{f_i^2}\frac {u_{i-1,j+1}-2u_{i-1,j}+u_{i-1,j-1}}{\Delta y^2}, \\[-28pt] \nonumber \end{align}
\begin{align} u_{i-1,j}\frac {\theta _{i,j}-\theta _{i-1,j}}{\Delta x}+v_{i-1,j}\frac {\theta _{i-1,j+1}-\theta _{i-1,j-1}}{2\Delta y} & = \frac {1}{f_i^2}\frac {\theta _{i-1,j+1}-2\theta _{i-1,j}+\theta _{i-1,j-1}}{\Delta y^2}, \end{align}
Given inlet conditions provided on
$u_{1,\,j}$
and
$\theta _{1,\,j},$
we first obtain the value of
$p_1$
by initialising a value, say
$p_1^*$
. This value is then used to obtain an initial value of
$v_{1,\,j}^*$
for all
$j,$
and the discrepancy of the value
$v_{1,\,N}^*$
with respect to the no-slip condition
$v_{1,\,N}^*=0$
is assessed. If this discrepancy is greater than a fixed tolerance, then
$p_1 \neq p_1^*$
and another value of
$p_1^*$
is chosen. This process is continued until the aforementioned discrepancy is below a fixed tolerance. Due to the linear dependence on pressure in the momentum equation here, a single iteration will obtain
$p_1.$
This whole process then continues until
$i=M,$
which produces the solution at the first time step. The solution of
$f_i$
is then obtained by discretising the Stefan condition via a forward Euler method. In figure 6, we provide supplementary evidence of the convergence of the numerical solution of
$f$
for increasing values of
$N$
, with
$M=N.$
Numerical solution of
$f$
from (2.4) and (2.5) for seven values of grid point number
$N$
uniformly distributed between
$300$
and
$900$
, with
$M=N.$
The arrow points in the direction of increasing
$N.$
The inset highlights the minimum of
$f$
. Parameter values:
$\lambda = 15/4,\,\phi =-5$
and
$\sigma =1$
.
Appendix B. Asymptotic analysis about the wall temperature drop point
In this appendix, we document the behaviour of the model from (2.4) and (2.5) asymptotically close to the point at which the fluid temperature at the wall,
$y=f,$
drops from
$\theta =0$
to
$\theta =-1$
at
$x=\delta .$
To do this, it is instructive to first outline the model behaviour in the warm-wall region for
$x\in [0,\,\delta ].$
We assume inlet conditions at
$x=0$
of the form (3.1) for simplicity.
The velocities in the region
$x\in [0,\,\delta ]$
admit the analytical solutions
$v=0$
and
$u_0=\lambda (1-y^2).$
Denoting
$\theta :=\theta _i$
when
$x\in [0,\,\delta ]$
, the thermal equation from (2.4) reduces to the separable equation
\begin{equation} \sigma \lambda (1-y^2)\frac {\partial \theta _i}{\partial x} = \frac {\partial ^2\theta _i}{\partial y^2}. \end{equation}
Using the boundary and inlet conditions from (2.5), the solution of (B1) is
\begin{equation} \theta _i = \sum _{n=1}^\infty A_n\psi _n(y;\,\gamma _n)\exp \bigg ({-\frac {\gamma _n x}{\lambda \sigma }}\bigg ), \end{equation}
such that
\begin{align} \psi _n &= \frac {1}{\sqrt {y}}\bigg [ {W}_{\frac {\sqrt {\gamma _n}}{4},\,\frac {1}{4}}(\sqrt {\gamma _n}y^2)+\alpha _n{M}_{\frac {\sqrt {\gamma _n}}{4},\,\frac {1}{4}}(\sqrt {\gamma _n}y^2) \bigg ], \end{align}
\begin{align} \alpha _n &= - \frac {\sqrt {\pi }\,\,\left [4\,\varGamma \left (\frac {1}{4}-\frac {\sqrt {\gamma _n}}{4}\right )+\varGamma \left (-\frac {\sqrt {\gamma _n}}{4}-\frac {3}{4}\right )+\sqrt {\gamma _n}\,\varGamma \left (-\frac {\sqrt {\gamma _n}}{4}-\frac {3}{4}\right )\right ]}{\varGamma \left (\frac {1}{4}-\frac {\sqrt {\gamma _n}}{4}\right )\varGamma \left (-\frac {\sqrt {\gamma _n}}{4}-\frac {3}{4}\right )}, \end{align}
Here, W and M are the Whittaker-W and -M functions, respectively, and
$\varGamma$
is the Gamma function. The eigenvalues
$\gamma _n$
and coefficient
$A_n$
are determined via
\begin{align} 0 &= {W}_{\frac {\sqrt {\gamma _n}}{4},\,\frac {1}{4}}(\sqrt {\gamma _n})+\alpha _n{M}_{\frac {\sqrt {\gamma _n}}{4},\,\frac {1}{4}}(\sqrt {\gamma _n}), \end{align}
\begin{align} A_n & = \frac {\phi }{\int _0^1 (1-s^2) \psi _n^2\, \mathrm{d}s}\int _0^1 (1-s^2)^2\psi _n \,\mathrm{d}s. \end{align}
Turning to the behaviour close to the channel wall, we introduce the new coordinates
$0\lt \tilde {Y}:=1-y\ll 1$
and
$0\lt \tilde {X}:=x-\delta \ll 1.$
When
$f=1-o(1),$
we have that
$v\rightarrow 0$
and
$u\sim 2\lambda \tilde {Y}$
, so that (2.4c
) provides
\begin{equation} 2\lambda \sigma \tilde {Y} \frac {\partial \theta }{\partial \tilde {X}}= \frac {\partial ^2 \theta }{\partial \tilde {Y}^2}, \end{equation}
to leading order. The solution of
$\theta$
in the cool-wall region is now necessarily
$\mathcal{O}(1)$
to facilitate the application of the boundary condition
$\theta =-1$
at
$\tilde {Y}=0.$
As
$\tilde {X}\rightarrow 0$
,
$\theta$
must match back to the warm-wall region where
$\theta =\theta _i$
is small. This motivates the expansion
\begin{equation} \theta \sim \varPsi (\eta )-\omega \tilde {Y}, \quad \text{where} \, \omega :=\lim _{x\rightarrow \delta ^-} \frac {\partial \theta _i}{\partial y}\bigg \rvert _{y=1}, \end{equation}
with
$\eta = \tilde {Y}/\tilde {X}^{1/3}$
so that (B3) provides
\begin{equation} -\frac {2\lambda \sigma }{3}\eta ^2\varPsi _\eta = \varPsi _{\eta \eta }. \end{equation}
Using the associated boundary condition
$\varTheta (0)=-1$
and matching condition
$\varTheta \rightarrow 0$
as
$\eta \rightarrow \infty ,$
the solution of (B5a
) is
\begin{equation} \varPsi = -\frac {1}{\varGamma (1/3)}\varGamma _i\bigg (\frac {1}{3},\,\frac {2\lambda \sigma }{9}\eta ^3 \bigg ), \end{equation}
where
$\varGamma _i$
is the incomplete gamma function.
Now, turning back to the original coordinates, the Stefan condition from (2.5f) when
$f=1+o(1)$
is given by
\begin{equation} \frac {\partial f}{\partial t}= -\frac {\partial \theta }{\partial y} \quad \text{at} \, y=1. \end{equation}
Using (B5b
) and integrating (B6a
) in
$t$
, we obtain
\begin{equation} f \sim 1+\bigg [\frac {3}{\varGamma (1/3)}\bigg (\frac {2\lambda \sigma }{9(x-\delta )}\bigg )^{1/3}-\omega \bigg ]t \quad \text{if} \, f\lt 1 \end{equation}
and
$f=1$
otherwise. Therefore, if
$\omega =\mathcal{O}(1)$
, then the first term in the brackets in (B6b
) dominates and ice will form. If the two terms are comparable, then the ice leading edge position
$x:=x_{LE}$
is determined as
\begin{equation} x_{LE} \sim \delta + \frac {2\lambda \sigma }{9}\bigg (\frac {3}{\varGamma (1/3)\omega }\bigg )^3. \end{equation}
In figure 7, we present the ice leading edge obtained numerically from (2.4) and (2.5) versus that in (B7) as a function of
$|\phi |$
, for fixed
$\lambda$
and
$\sigma .$
The bracketed terms in (B6b
) are expected to be of comparable order when
$\phi \gg 1$
(since
$\omega \propto |\phi |)$
, and thus the accuracy of (B7) improves in this limit, as observed.
Appendix C. Numerical results and methodology for the complete closure problem from (4.4)
We now present and discuss the numerical solutions of the scaled closure system from (4.4). We solve this system numerically using finite differences, typically using step lengths
$10^{-4}$
and
$10^{-2}$
in
$X$
and
$\zeta$
, respectively. Fixing the values of
$\sigma ,\,F_0$
and
$\mathcal{C}$
, the functions
$R$
and
$T,$
as well as the value of
$Q$
and
$b,$
were determined via the methods described in § 4.2.3. The requirement that
$b\lt 1$
was ensured; see § 4.1 for further details. Using the expansions from (4.5) as starting conditions at some sufficiently large value of
$-X$
(typically
$X=-10),$
the equations for
$\bar {U}$
and
$\bar {V}$
from (4.4) were solved via a semi-implicit, forward marching finite difference method (Timoshin & Smith Reference Timoshin and Smith1995). This was accompanied by a similar formulation for the thermal equation from (4.4). The numerical solution of
$F$
was obtained via the discretised version of (4.4f
). The finite difference scheme was implemented until a point far downstream where the linear asymptotes in
$F$
became clear.
In figure 8, we present some exemplar numerical solutions and quantities arising from (4.4), for three different parameter sets
$[\sigma ,\,F_0,\,\mathcal{C}]$
which are presented in the caption. For each different parameter set,
$F(X)$
and
$P(X)$
are presented in panels (a,c,e), whilst the longitudinal velocity and temperature at the centreline
$\zeta =0$
,
$\bar {U}(X,\,0)$
and
$\varTheta (X,\,0)$
, as well as the wall shear stress
$- {\partial \bar {U}}/{\partial \zeta }(X,\,1)$
are presented in panels (b,d, f).
Numerical solutions and quantities of the full scaled closure system from (4.4) as functions of
$X$
, for the three parameter sets
$[\sigma ,\,F_0,\,\mathcal{C}]$
with
$a={1}/{2}.$
The top panels for each case show
$F(X)$
and the pressure
$P(X).$
The bottom panels show the longitudinal velocity
$\bar {U}$
and temperature
$\varTheta$
at the centreline
$\zeta =0$
, as well as the wall shear stress,
$\partial \bar {U}/\partial \zeta$
at
$\zeta =1.$
The specific parameters are: (a, b)
$[8,\, {3}/{2},\,{1}/{4}]$
, (c, d)
$[{1}/{5},\, 3,\,4]$
and (e, f)
$[1,\, 6,\,2]$
.
In all results, the observed trends agree closely with the asymptotic findings presented in § 4.1, as well as with the numerical simulations of (2.4) and (2.5) discussed in § 3. As the axial position
$X$
increases from the initial reference point
$(X=-4)$
, the gap width
$F$
initially decreases to a finite minimum before increasing monotonically. Approximately linear growth in
$F$
at large
$-X$
is observed, indicative of a distinct upstream Jeffery–Hamel flow regime (Jeffery Reference Jeffery1915; Hamel Reference Hamel1917). Similarly, the longitudinal water velocity
$U$
exhibits near-linear growth in for large
$-X$
, while the temperature profile asymptotically approaches zero as
$X$
increases. This thermal behaviour is consistent with the numerical results presented in § 3, where the fluid temperature is seen to approach the melt temperature (
$\theta = -1$
) far upstream.
More detailed comments pertaining to figure 8 are now given. The example from figure 8(a,b) uses a relatively large Prandtl number and also has
$CF_0=3/8$
, which is close to the limiting case of lubrication flow discussed in § 4.2.1. We obtain numerically from (4.7) the unique eigenvalue
$b_{\text{s}}=-0.37$
, which agrees with the approximation from (4.10) up to three significant figures, although this is unsurprising in view of the results presented in figure 4. The numerical solution shows that
$F$
grows almost linearly for sufficiently large
$|X|,$
with a smooth transition through its minimum. The temperature perturbation at the centreline decays towards zero as expected, while the centreline velocity and wall shear stress remain nearly constant. The pressure gradient remains favourable throughout, reflecting the expected features of a lubrication-dominated flow.
In contrast to the results discussed previously, those in figure 8(c,d) have a comparatively small value of
$\sigma$
and
$CF_0=12$
, suggesting that solutions will admit plug-flow attributes. We obtain numerically from (4.7)
$b_{\text{s}}=-0.14$
, which agrees excellently with the plug-flow approximation of
$b_{\text{s}}$
from (4.15). The incoming centreline velocity and wall shear stress have reasonably high values; however, the positive slope of
$F$
is comparatively low and the wall shear decreases significantly after the gap minimum but remains positive. The pressure gradient becomes very favourable near the gap minimum, inducing a substantial decrease in the pressure locally as the fluid passes through the gap. Furthermore, the variation in centreline temperature remains rather small after passing through the gap minimum.
In figure 8(e,f), we present results where
$\sigma =1$
and
$CF_0=12.$
Despite the non-extreme value of
$\sigma ,$
the value of
$CF_0$
suggests that plug-flow like features should appear and this is found to be so. In addition, the relatively large negative slope
$F_0$
prescribed upstream is seen to produce large positive slopes upstream, resulting in a significant expansion of the ice gap. The centreline velocity suggests flow acceleration throughout, but however the wall shear stress decreases monotonically with
$X$
which eventually results in a sign change, suggesting slow flow separation there. This separation results in numerical instabilities just beyond
$X=2$
. This slow flow near the ice surface as separation develops explains the significant acceleration of the centreline flow, which is required to conserve the total mass flux across the ice gap minimum. The pressure behaviour is consistent with these observations, as a favourable pressure gradient is encountered ahead of the gap, but followed by an adverse pressure gradient which arises soon after the gap width minimum.
Open access
