1. Introduction
In the modern world, there is an increasing demand for ultra-thin glass sheets, with numerous applications, including computer and television screens, mobile phones and microchips (Auch et al. Reference Auch, Soo, Ewald and Soo-Jin2002; Plichta, Weber & Habeck Reference Plichta, Weber and Habeck2003; Tanaka Reference Tanaka2013; Choi, Yun & Park Reference Choi, Yun and Park2016; Macrelli, Varshneya & Mauro Reference Macrelli, Varshneya and Mauro2020). The most prevalent manufacturing technique for uniform, flat glass sheets is the float glass process (Nascimento Reference Nascimento2014), in which molten glass floats on a bed of molten metal, typically tin (Pilkington Reference Pilkington1969; Charnock Reference Charnock1970). Although the float glass process can be used to produce ultra-thin glass sheets (Shou et al. Reference Shou, Hongcan, Xin and Yong2019), the inherent asymmetry, whereby one side of the glass is exposed to the air whilst the other rests on the molten metal, can lead to undesirable surface anomalies (Frischat Reference Frischat2002; Tiwari et al. Reference Tiwari, Mondal, Kumar, Agarwal, Dixit and Sharma2011). An alternative method that produces thin glass sheets in a symmetric process is the overflow fusion process, used extensively by Corning Inc. and other glass manufacturers (Lin & Chang Reference Lin and Chang2007b ; Bohun et al. Reference Bohun, Breward, Cummings and Witelski2010; Ellison & Cornejo Reference Ellison and Cornejo2010; Tanaka Reference Tanaka2013).
In overflow fusion, a triangular wedge shaped trough is filled with molten glass via an inlet so that once the trough fills, the glass overflows and flows down the underside of the trough, as shown in figure 1. The two layers meet at the bottom of the trough, where they coalesce to form a single viscous sheet (Dockerty Reference Dockerty1967; Lin & Chang Reference Lin and Chang2007a ; Ellison & Cornejo Reference Ellison and Cornejo2010). The resulting sheet falls under gravity and thins until it reaches the desired thickness, after which it is cooled to form a solid, uniform sheet of thin glass with both surfaces exposed to the same conditions.
Dockerty’s original 1967 patent (Dockerty Reference Dockerty1967), illustrating the sheet forming apparatus.
This process has two qualitatively different flow regimes of interest: the first where the molten glass flows on the underside of the trough, and the second where the glass leaves the trough and falls under gravity. The flow in the first stage falls under the classical regime of lubrication theory (Ockendon & Ockendon Reference Ockendon and Ockendon1995; Myers Reference Myers1998; Lin & Kondic Reference Lin and Kondic2010), while the flow in the second stage is extensional (Dewynne, Ockendon & Wilmott Reference Dewynne, Ockendon and Wilmott1992; Howell Reference Howell1996). These regimes are both characterised by a small aspect ratio (i.e. the flow is ‘long and thin’), with lubrication theory describing thin film flow on a substrate and extensional flow describing uniform ‘plug’ flow between two free surfaces. Both of these stages were explored separately, in relation to overflow fusion, in a 2010 study group report for Corning Inc. (Bohun et al. Reference Bohun, Breward, Cummings and Witelski2010). The main objective of this paper is to show how a flow transitions between these two well-studied regimes. As a result, we not only model overflow fusion, but also explore the solution of a new canonical fluid mechanics problem of general mathematical interest. In particular, the work presented here is relevant to other problems concerning a transition between no-slip and free-slip conditions, such as the transition of an ice sheet from grounded to floating (see e.g. Barcilon & MacAyeal Reference Barcilon and MacAyeal1993; Chugunov & Wilchinsky Reference Chugunov and Wilchinsky1996).
In § 2, we present a mathematical model for the overflow fusion process, focusing on the simple case of zero wedge angle so that in the first stage, the molten glass flows down a vertical wall. In § 3, we determine the leading-order behaviour of the model in the limit of small capillary number. We decompose the domain into four regions of interest, and determine the leading-order solution in each region by matching the solutions in adjoining regions. The downstream behaviour of the leading-order solution reveals appropriate far-field conditions to impose on the full problem. In § 4, we combine the results from each region to obtain a uniformly valid reduced model for the film thickness that holds over the full length of the flow. Finally, in § 5, we discuss the implications of the results, and consider areas for future research.
2. Problem set-up
2.1. Dimensional problem
We consider the simplified problem in which the angle at the apex of the trough is zero, so the outer surface of the trough is simply an infinitesimally thin vertical wall, as shown in figure 2. We use a two-dimensional coordinate system with the
$x$
-axis along the wall and pointing in the direction of gravity
$\boldsymbol{g}$
, the
$z$
-axis pointing in the transverse direction, and the origin at the apex of the trough. The free surface of the glass is denoted by
$z=h(x,t)$
, where
$t$
is time. Due to the symmetry of the problem, we only consider the fluid flow in
$0 \le z \le h(x,t)$
.
Schematic of the dimensional mathematical model, indicating the coordinate system, direction of gravity, and incoming flux.
Although molten glass has been shown to exhibit visco-elastic behaviour (Simmons et al. Reference Simmons, Mohr and Montrose1982, Reference Simmons, Ochoa, Simmons and Mills1988; Bohun et al. Reference Bohun, Breward, Cummings and Witelski2010), we follow previous authors Taroni et al. (Reference Taroni, Breward, Cummings and Griffiths2013), O’Kiely et al. (Reference O’Kiely, Breward, Griffiths, Howell and Lange2018) and Lin & Chang (Reference Lin and Chang2007a
) in modelling glass as a Newtonian viscous fluid. Furthermore, although molten glass has a strongly temperature-dependent viscosity (O’Kiely et al. Reference O’Kiely, Breward, Griffiths, Howell and Lange2019), since our focus is on the fundamental fluid mechanics of the problem, for simplicity we do not include any temperature dependence, and therefore assume that the fluid has constant density
$\rho$
, dynamic viscosity
$\mu$
, and surface tension
$\gamma$
. The fluid velocity
$\boldsymbol{u}=(u,w)$
and pressure
$p$
therefore satisfy the incompressible Navier–Stokes equations
\begin{align} \boldsymbol{\nabla \boldsymbol{\cdot }u} &= 0 , \end{align}
\begin{align} \rho \left ( \frac {\partial {\boldsymbol{u}}}{\partial {t}} + \left ( \boldsymbol{u \boldsymbol{\cdot }\boldsymbol{\nabla }} \right ) \boldsymbol{u} \right ) &= -{\boldsymbol{\nabla }} p + \mu\, {\nabla} ^2 \boldsymbol{u} + \rho \boldsymbol{g}. \\[9pt] \nonumber \end{align}
At the free surface, we impose the kinematic condition and stress balance, namely
\begin{align} w&=h_t + \textit{uh}_x , \end{align}
\begin{align} \gamma \kappa \boldsymbol{n} &= \boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{n}, \\[9pt] \nonumber \end{align}
where partial derivatives are denoted by subscripts,
\begin{equation} \boldsymbol{\sigma } = \begin{pmatrix} -p + 2 \mu u_x && \mu (u_z + w_x) \\ \mu (u_z + w_x) && -p + 2 \mu w_z \end{pmatrix} \end{equation}
is the Newtonian stress tensor,
\begin{align} \boldsymbol{n} = \frac {1}{(1+h_x^2)^{1/2}} \begin{pmatrix} -h_x \\ 1 \end{pmatrix}, \quad \kappa = \frac {h_{\textit{xx}}}{(1+h_x^2)^{3/2}} \end{align}
are the unit normal and curvature of the free surface, respectively, and we have, without loss of generality, taken the atmospheric pressure to be zero.
Along the wall, we impose no flux and no slip, and downstream of the apex along the centreline, we impose symmetry. Thus at
$z=0$
,
\begin{align} u=w=0 \quad\text{for }x\lt 0 , \end{align}
\begin{align} u_z=w = 0\quad\text{for }x\gt 0, \end{align}
and we expect lubrication flow upstream and extensional flow downstream of the apex.
We assume that the flow far upstream is fully developed, with a specified incoming flux
$Q$
between one side of the wall and the free surface, so that the upstream flow satisfies
\begin{align} \int _0^{h(x)} u(x,z)\, \textrm {d}{z} \rightarrow Q \quad \text{as } x \rightarrow -\infty . \end{align}
Solving (2.1), (2.2), (2.5a ) and (2.6), we obtain the upstream far-field behaviour
\begin{align} h(x)\rightarrow h_{{s}}=\left (\frac {3\mu Q}{\rho g}\right )^{1/3}, \quad u(x,z)&\rightarrow \frac {\rho g}{\mu }\left (h_{{s}}z-\frac {z^2}{2}\right ) \quad\text{as }x\rightarrow -\infty , \end{align}
with
$w$
and
$p$
both tending to zero, where
$g=|\boldsymbol{g}|$
.
2.2. Non-dimensionalisation
Before establishing the aspect ratio of the flow in each region, we non-dimensionalise (2.1), (2.2), (2.5) and (2.7) by scaling:
\begin{align} & x = h_{{s}}x^*, \quad\,\, z= h_{{s}}z^*, \quad\,\,\, h = h_{{s}}h^*, \\ & u = U_{{s}}u^*, \quad w = U_{{s}} w^*, \quad p = (\mu U_{{s}}/h_{{s}})p^*, \quad t = (h_{{s}}/U_{{s}})t^*,\nonumber \end{align}
where
$h_{{s}}$
is the upstream film thickness, and
$U_{{s}} = Q/h_{{s}}$
is the upstream average velocity. We assume that the Reynolds number
${Re}=\rho Q/\mu \ll 1$
so that inertia is negligible, which is appropriate for glass manufacturing processes.
Upon dropping stars, we thus obtain the dimensionless Stokes equations
\begin{align} u_x + w_z &= 0, \end{align}
\begin{align} p_x &= u_{\textit{xx}} + u_{zz} + 3, \end{align}
\begin{align} p_z &= w_{\textit{xx}} + w_{zz}, \\[7pt] \nonumber \end{align}
with boundary conditions at the free surface
$z=h(x,t)$
,
\begin{align} w &= h_t + \textit{uh}_x, \end{align}
\begin{align} \frac {1}{\textit{Ca}}\, \frac {h_{\textit{xx}}}{(1+h_x^2)^{3/2}} &= -p -2u_x\,\frac {1+h_x^2}{1-h_x^2}, \end{align}
\begin{align} 4u_xh_x &= (u_z + w_x)\left(1-h_x^2\right), \end{align}
where
$\textit{Ca}=\mu U_{{s}}/\gamma$
, and at
$z=0$
,
\begin{align} u = w = 0 \quad \text{for }x\lt 0, \end{align}
\begin{align} u_z = w = 0 \quad \text{for }x\gt 0, \end{align}
with far-field behaviour
\begin{align} h(x)\rightarrow 1, \quad u(x,z)\rightarrow 3\left (z-\frac {z^2}{2}\right ), \quad w(x,z),\,p(x,z)\rightarrow 0 \quad\text{as }x\rightarrow -\infty . \end{align}
From here on, we assume that the system is in steady state, and from (2.9a ), (2.10a ), (2.11) and (2.12), we derive the net mass conservation equation
\begin{equation} \int _0^{h(x)}u(x,z)\,\mathrm{d}z=1. \end{equation}
By depth integrating (2.9b ) and applying the boundary conditions (2.10c ) and (2.11) (for full details, see Appendix A), we derive the axial force balance
\begin{equation} \frac {\textrm {d}{T}}{\textrm {d}{x}} = u_{z}|_{z=0}- 3h, \end{equation}
where
\begin{equation} T = \left (\int _0^h \boldsymbol{\sigma \boldsymbol{\cdot }i}\, \textrm {d}{z} + \frac {1}{\textit{Ca}}\boldsymbol{t}\right ) \boldsymbol{\cdot } \boldsymbol{i} \end{equation}
is the net tension within the sheet, and
\begin{equation} \boldsymbol{t} = \frac {1}{(1+h_x^2)^{1/2}} \begin{pmatrix} 1 \\ h_x \end{pmatrix} \end{equation}
is the unit tangent to the free surface. The right-hand side of (2.14) represents the net vertical traction acting on the sheet, comprised of the drag at the wall (which is only present for
$x\lt 0$
) and the net gravitational forcing due to the sheet’s weight.
To form a well-posed boundary-value problem, the system (2.9)–(2.12) must be supplemented with downstream boundary conditions. In § 3.4, we show how the downstream behaviour depends on the tension applied to the sheet. We choose to focus on the case where the sheet falls freely under gravity, with no additional forcing being applied. We thus impose
\begin{equation} T \rightarrow \frac {1}{\textit{Ca}} \quad \text{as } x \rightarrow \infty , \end{equation}
which selects a unique solution to the problem. We show in § 3.4 that the corresponding downstream behaviour is
\begin{align} h(x)\sim \frac {8}{3\left (x-x_*\right )^2},\quad u(x,z)\sim \frac {1}{h(x)}, \quad w(x,z)\sim \frac {z\,h'(x)}{h(x)^2} \quad\text{as }x\rightarrow \infty , \end{align}
where
$x_*$
is to be determined as part of the solution in order to satisfy (2.12).
We are left with a single dimensionless parameter in the model (2.9)–(2.18), namely the capillary number
\begin{equation} \textit{Ca}=\frac {\mu U_{{s}}}{\gamma } =\left (\frac {\rho g\mu ^2Q^2}{3\gamma ^3}\right )^{1/3} = \frac {h_{{s}}^2}{3\ell _{{c}}^2}, \end{equation}
where
$\ell _{{c}}=\sqrt {\gamma /\rho g}$
is the capillary length, which is approximately
$4\,\mathrm{mm}$
for the typical parameter values given in table 1. Although typical industrial values of the flux
$Q$
, and hence
$h_{{s}}$
and
$U_{{s}}$
, are not publicly available, we focus on the limit of small
$\textit{Ca}$
. This limit is relevant to the manufacture of ultra-thin glass sheets, when the film flowing down the substrate is much thinner than
$\ell _{{c}}$
. We aim to recover the shape of the free surface
$h(x)$
and the structure of the solution in the limit
$\textit{Ca} \rightarrow 0$
, thus we define the small parameter
\begin{equation} \epsilon = (3\,\textit{Ca})^{1/3}\ll 1, \end{equation}
which we identify in § 3.1 as the aspect ratio of the lubrication flow region, a classic scaling in low capillary number thin film flow (Bretherton Reference Bretherton1961).
Table of typical parameter values for molten glass.
We decompose the flow into four regions, as illustrated in figure 3. For
$x\lt 0$
we have a region of lubrication flow, near
$x=0$
we have a region of fully two-dimensional Stokes flow, and for
$x\gt 0$
we have two regions of extensional flow. The scalings indicated in figure 3 are derived in § 3 using matched asymptotic expansions.
In the lubrication flow region, surface tension, viscosity and gravity all contribute to the flow, which we find varies over a dimensionless length scale of order
$\textit{Ca}^{-1/3}$
. In the Stokes flow region, the boundary condition (2.11) transitions from no-flux/no-slip to a symmetry condition, and the flow adjusts across a region of
$O(1)$
aspect ratio. The flow in this region is governed by the full two-dimensional Stokes equations. Just downstream of
$x=0$
, we have a region of extensional flow whose length we find scales with
$\textit{Ca}^{-1/3}$
, in which surface tension and gravity balance. In this region, the height of the free surface decreases towards zero. As the film gets very thin, viscous effects come into play and, as a result, we enter a second region of extensional flow in which surface tension, gravity and viscous effects all balance. In this region, we find that the film thickness is of order
$\textit{Ca}^{2/5}$
, and the flow varies over a length scale of order
$\textit{Ca}^{-1/5}$
.
We start by considering the upstream lubrication flow. We then determine the solution in each subsequent region by matching to the solution in the region immediately upstream. This approach allows us to determine the solution across the entire flow domain.
3. Flow regions
3.1. Lubrication flow
For
$x \lt 0$
, we have a typical lubrication flow problem. We suppose that the flow is ‘long and thin’ in this region, and scale
\begin{align} x = \tilde {x}/\epsilon,\quad z = \tilde {z},\quad h = \tilde {h},\quad u = \tilde {u},\quad w = \epsilon \tilde {w},\quad p = \tilde {p}/\epsilon, \end{align}
where we use tildes to denote variables in the lubrication regime. The scalings of
$w$
and
$p$
are introduced to maintain dominant balances in the incompressibility condition (2.9a) and momentum equation (2.9b), respectively. Recalling our definition (2.20) of
$\epsilon$
, we see that the flow in this region varies over length scale
$\textit{Ca}^{-1/3}$
, as indicated in figure 3. Substituting into (2.9)–(2.12), we obtain the scaled equations
\begin{align} \tilde {u}_{\tilde {x}} + \tilde {w}_{\tilde {z}} &= 0, \end{align}
\begin{align} \tilde {p}_{\tilde {x}} &=\epsilon ^2\tilde {u}_{\tilde {x}\tilde {x}} +\tilde {u}_{\tilde {z}\tilde {z}} + 3, \end{align}
\begin{align} \tilde {p}_{\tilde {z}}&=\epsilon ^4\tilde {w}_{\tilde {x}\tilde {x}}+\epsilon ^2\tilde {w}_{\tilde {z}\tilde {z}}, \end{align}
with boundary conditions at
$\tilde {z}=\tilde {h}(\tilde {x})$
,
\begin{align} \tilde {w} &= \tilde {u} {\tilde {h}}' , \end{align}
\begin{align} \frac {3{\tilde {h}}''}{(1+\epsilon ^2{{\tilde {h}}'}{}^2)^{3/2}} &= -\tilde {p} -2\epsilon ^2\tilde {u}_{\tilde {x}}\,\frac {1+\epsilon ^2{{\tilde {h}}'}{}^2}{1-\epsilon ^2{{\tilde {h}}'}{}^2} , \end{align}
\begin{align} 4\epsilon ^2\tilde {u}_{\tilde {x}}{\tilde {h}}' &= (\tilde {u}_{\tilde {z}} + \epsilon ^2\tilde {w}_{\tilde {x}})(1-\epsilon ^2{{\tilde {h}}'}{}^2), \end{align}
in which primes denote derivatives, and at
$\tilde {z}=0$
,
\begin{align} \tilde {u} = \tilde {w}=0, \end{align}
with far-field behaviour
\begin{align} \tilde {h}(\tilde {x})\rightarrow 1,\quad \tilde {u}(\tilde {x},\tilde {z})\rightarrow 3\left (\tilde {z}-\frac {\tilde {z}^2}{2}\right ),\quad \tilde {w}(\tilde {x},\tilde {z}),\quad \tilde {p}(\tilde {x},\tilde {z})\rightarrow 0\quad \text{as }x&\rightarrow -\infty . \end{align}
Solving (3.2)–(3.4), we find that the leading-order pressure and longitudinal velocity are given by
\begin{align} \tilde {p}=-3\tilde {h}'', \quad \tilde {u}= 3\left (1+\tilde {h}'''\right )\left (\tilde {h}\tilde {z}-\frac {\tilde {z}^2}{2}\right ), \end{align}
thus the net mass conservation law (2.13) provides the equation for
$\tilde {h}$
, namely
\begin{equation} \tilde {h}^3\left (1 + \tilde {h}'''\right ) = 1. \end{equation}
This is a typical governing equation for a surface-tension-dominated thin film (see e.g. Tuck & Schwartz Reference Tuck and Schwartz1990; Myers Reference Myers1998).
Equation (3.7) must be solved subject to the matching condition
$\tilde {h}(\tilde {x})\rightarrow 1$
as
$\tilde {x}\rightarrow -\infty$
. By linearising about
$\tilde {h}=1$
, we find that there is a two-parameter family of such solutions, with the far-field behaviour
\begin{equation} \tilde {h}(\tilde {x})\sim 1+A\exp \left (\frac {3^{1/3}\tilde {x}}{2}\right )\cos \left (\frac {3^{5/6}(\tilde {x}-\textit{c})}{2}\right )\quad \text{as }\tilde {x}\rightarrow -\infty , \end{equation}
in which the constants
$A$
and
$c$
are a priori unknown. Two more pieces of information are therefore required to close the problem in this region. By matching with the downstream solution, in § 3.4 we derive the boundary conditions
\begin{align} \tilde {h}(0)=\frac {1}{6}\tilde {h}''(0)^3,\quad \tilde {h}'(0)=-\frac {1}{2}\tilde {h}''(0)^2 \end{align}
for
$\tilde {h}$
.
We solve the problem (3.7)–(3.9) numerically using a shooting method, by imposing the conditions (3.9) at
$\tilde {x}=0$
and varying the a priori unknown value of
$\tilde {h}''(0)$
to satisfy
$\tilde {h}(\tilde {x}) \rightarrow 1$
as
$\tilde {x} \rightarrow -\infty$
. By following this procedure, we evaluate
$\tilde {h}''(0)\approx 1.3389$
and obtain the thickness profile shown in figure 4 in the lubrication region. Using our numerical solution, we also determine the values of the undetermined constants in (3.8), namely
$A \approx -1.03$
,
$c \approx 0.282$
.
3.2. Stokes flow
In the Stokes flow region around
$x=0$
, we return to the unscaled dimensionless variables. To match with the solution (3.6) in the lubrication region, we must scale the pressure with
$1/\epsilon$
. We therefore define
$p = \mathfrak{p}/\epsilon$
to obtain the governing equations
\begin{align} u_x + w_z &= 0, \end{align}
\begin{align} \mathfrak{p}_x &= \epsilon u_{\textit{xx}} + \epsilon u_{zz} + 3\epsilon , \end{align}
\begin{align} \mathfrak{p}_z &= \epsilon w_{\textit{xx}} + \epsilon w_{zz}, \\[7pt] \nonumber \end{align}
with boundary conditions at
$z=h(x)$
,
\begin{align} w &= \textit{uh}^{\prime}, \end{align}
\begin{align} \frac {3h^{\prime\prime}}{(1+h^{\prime 2})^{3/2}} &= -\epsilon ^2\mathfrak{p} -2\epsilon ^3 u_x\,\frac {1+h^{\prime 2}}{1-h^{\prime 2}}, \end{align}
\begin{align} 4u_xh' &= (u_z + w_x)(1-h'^2), \\[9pt] \nonumber \end{align}
and at
$z=0$
,
\begin{align} u = w= 0 \quad \text{for }x\lt 0, \end{align}
\begin{align} u_z = w = 0\quad \text{for }x\gt 0. \end{align}
We now expand the dependent variables as asymptotic expansions of the form
$h(x)\sim h_0(x)+\epsilon\, h_1(x)+\cdots$
. From the normal stress balance (3.11b
), we see that
$h_0$
and
$h_1$
must be linear functions of
$x$
, hence by matching with
$\tilde {h}(\tilde {x})$
, we find
\begin{align} h_0(x)=\tilde {h}(0^-), \quad h_1(x)=x\,{\tilde {h}}'(0^-)+a_0, \end{align}
where
$a_0$
is an integration constant. Similarly, from (3.10b
), and (3.10c
), we see that
$\mathfrak{p}_0$
must be constant, and by matching with the solution (3.6), we have
\begin{equation} \mathfrak{p}_0=-3\,\tilde {h}''(0^-). \end{equation}
Then, from the
$O(\epsilon ^2)$
terms in (3.11b
), we find
\begin{equation} h_2(x)=\frac {1}{2}\tilde {h}''(0^-)\,x^2+a_1 x+a_2, \end{equation}
where
$a_1$
and
$a_2$
are further integration constants, which gives
\begin{align} h(x) \sim \tilde {h}(0^-) + \epsilon (x\,\tilde {h}'(0^-)+a_0)+\epsilon ^2\left (\frac {1}{2}\tilde {h}''(0^-)\,x^2+a_1 x+a_2\right ) + \cdots . \end{align}
Now proceeding to
$O(\epsilon )$
in (3.10b
), and (3.10c
), we get the full two-dimensional Stokes equations
\begin{align} {u_0}_x + {w_0}_z &= 0, \end{align}
\begin{align} {\mathfrak{p}_1}_x &= {u_0}_{\textit{xx}} + {u_0}_{zz} + 3, \end{align}
\begin{align} {\mathfrak{p}_1}_z &= {w_0}_{\textit{xx}} + {w_0}_{zz}, \end{align}
subject to, at
$z=h_0$
,
\begin{align} w_0={u_0}_z&=0, \end{align}
and at
$z=0$
,
\begin{align} u_0 = w_0 = 0\quad\text{for }x\lt 0, \end{align}
\begin{align} {u_0}_z = w_0 = 0\quad\text{for }x\gt 0. \end{align}
By matching with the solution (3.6), considering the scaling for
$\tilde {w}$
in (3.1), and using (3.17b
), we determine the upstream far-field behaviour
\begin{align} u_0(x,z)\rightarrow \frac {3}{h_0^3}\left (h_0z-\frac {z^2}{2}\right ), \quad w_0(x,z)\rightarrow 0, \quad {\mathfrak{p}_1}_x(x,z)\rightarrow 3-\frac {3}{h_0^3} \quad \text{as }x\rightarrow -\infty . \end{align}
On the other hand, assuming that the far-field flow is fully developed at the downstream end of the domain, due to the symmetry condition at
$z=0$
, the flow approaches a uniform profile with
\begin{align} u_0(x,z)\rightarrow \frac {1}{h_0}, \quad w_0(x,z)\rightarrow 0, \quad {\mathfrak{p}_1}_x(x,z)\rightarrow 3 \quad \text{as }x\rightarrow \infty . \end{align}
Hence, using (3.14), we have that
\begin{align} \mathfrak{p} \rightarrow -3\,\tilde {h}''(0^-) + \epsilon (3x+a_3) + \cdots \quad \text{as }x\rightarrow \infty , \end{align}
where
$a_3$
is another integration constant.
As we show in the following subsections, it is possible to match with the downstream extensional flow without solving for the pressure and velocity in the Stokes flow region. Nevertheless, for completeness, we show in Appendix B how the solution can be obtained using the Wiener–Hopf method, following the methodology in Luchini (Reference Luchini1991). This analysis shows how the flow adjusts from a no-flux/no-slip condition on the substrate to a symmetry condition along the centreline, and confirms that the solution satisfies the far-field boundary conditions (3.20) and (3.21).
3.3. Capillary–gravity extensional flow
For
$x\gt 0$
, once we exit the two-dimensional Stokes flow region, we return to the ‘long and thin’ scalings employed in § 3.1. The governing equations (3.2) and free-surface boundary conditions (3.3) are unchanged, but the no-flux/no-slip condition (3.4) at
$\tilde {z}=0$
is modified to the symmetry condition
\begin{align} \tilde {u}_{\tilde {z}}=\tilde {w}=0. \end{align}
In this ‘capillary–gravity’ region, the leading-order flow is now extensional, and the film thickness is determined by a balance between surface tension and hydrostatic pressure, with viscous effects subdominant.
We find from (3.2b
), (3.2c
), (3.3c
) and (3.23) that the leading-order velocity
$\tilde {u}$
and pressure
$\tilde {p}$
are independent of
$\tilde {z}$
. Using the net mass conservation law (2.13) to obtain the leading-order velocity, we solve (3.3b
) and match with the far-field pressure (3.22) in the Stokes flow region to obtain
\begin{align} \tilde {u}(\tilde {x})=\frac {1}{\tilde {h}(\tilde {x})}, \quad \tilde {p}(\tilde {x})=3(\tilde {x}-\tilde {h}''(0^-)), \end{align}
thus the normal stress balance (3.3b ) reduces to
\begin{equation} {\tilde {h}}''(\tilde {x})=-\tilde {x}+\tilde {h}''(0^-). \end{equation}
By integrating (3.25) twice and matching with the film thickness (3.16) in the Stokes flow region as
$\tilde {x} \rightarrow 0$
, we find that the film thickness is given by
\begin{equation} \tilde {h}(\tilde {x})=-\frac {\tilde {x}^3}{6}+\frac {1}{2}\tilde {h}''(0^-)\,{\tilde {x}}^2+\tilde {h}'(0^-)\,\tilde {x}+\tilde {h}(0^-) \end{equation}
for
$\tilde {x}\gt 0$
. Thus we find that
$\tilde {h}$
and its first two derivatives are all continuous as we move between the lubrication region and the extensional flow region, although there remains a discontinuity in the third derivative, which is smoothed out over the Stokes flow region.
3.4. Capillary–viscous–gravity extensional flow
The film thickness given by (3.26) will pass through zero and become negative at some finite value of
$\tilde {x}$
, which we denote by
$\tilde {x}_0$
. To prevent this unphysical behaviour, viscous effects must come into play, and we therefore scale our variables so that the extensional viscous term
$\tilde {u}_{\tilde {x}\tilde {x}}$
balances the gravitational force in (3.2b
). We also impose that the surface tension term in (3.3b
) balances the pressure at leading order. As a result, (3.3b
) gives that
$\tilde {p}(\tilde {x}_0)$
must be zero, and hence (3.24) gives
\begin{align} \tilde {h}''(0) = \tilde {x}_0. \end{align}
The resulting scalings are shown in Appendix C to be
\begin{align} & \!\!\!\!\!\!\!\!\!\!\!\tilde {x}-\tilde {x}_0=\epsilon ^{2/5}X, \quad\,\,\,\,\,\tilde {z} = \epsilon ^{6/5} Z,\quad\,\,\,\,\tilde {h} = \epsilon ^{6/5}H, \nonumber\\ & \tilde {u} =\epsilon ^{-6/5} U, \quad \tilde {w}= \epsilon ^{-2/5}W, \quad \tilde {p}=\epsilon ^{2/5}P. \end{align}
Recalling that
$x=\tilde {x}/\epsilon$
, we see that these scalings correspond to
\begin{equation} x = \epsilon ^{-1}\tilde {x}_0+ \epsilon ^{-3/5}X, \end{equation}
thus the flow in the capillary–viscous–gravity region varies over length scale
$\epsilon ^{-3/5}\propto \textit{Ca}^{-1/5}$
, with
$h$
being of size
$\epsilon ^{6/5}\propto \textit{Ca}^{2/5}$
, as indicated in figure 3.
For the choice of scalings (3.28) to be consistent, the film thickness
$\tilde {h}$
in the capillary–gravity region must scale like
$(\tilde {x}-\tilde {x}_0)^3$
as
$\tilde {x} \rightarrow \tilde {x}_0$
. This condition determines that the leading-order solution (3.26) may also be written as
\begin{equation} \tilde {h}(\tilde {x}) = \frac {1}{6}\left (\tilde {x}_0-\tilde {x}\right )^3. \end{equation}
Recalling that
$\tilde {h}$
and its first two derivatives are continuous across
$\tilde {x}=0$
, we have from (3.26) and (3.30) that
\begin{align} \tilde {h}(0) = \frac {1}{6}\tilde {x}_0^3, \quad \tilde {h}'(0) = -\frac {1}{2}\tilde {x}_0^2, \quad \tilde {h}''(0) = \tilde {x}_0, \end{align}
the last of which is consistent with (3.27). These relations provide us with the boundary conditions (3.9) used to close the problem in the lubrication region. As discussed in § 3.1, the numerical solution of the lubrication flow problem determines
$\tilde {h}''(0)=\tilde {x}_0\approx 1.3389$
.
Under the scalings (3.28), the governing equations (3.2) and boundary conditions (3.3) and (3.23) are transformed to
\begin{align} U_X + W_Z &= 0, \end{align}
\begin{align} \epsilon ^{18/5}P_X &= \epsilon ^{18/5}U_{\textit{XX}} + U_{\textit{ZZ}} + 3\epsilon ^{18/5}, \end{align}
\begin{align} P_Z &= \epsilon ^{18/5}W_{\textit{XX}} + W_{\textit{ZZ}}, \end{align}
with boundary conditions at
$Z=H(X)$
,
\begin{align} W &= UH', \end{align}
\begin{align} \frac {3H''}{(1+\epsilon ^{18/5}H^{\prime 2})^{3/2}}&= -P-2U_X\frac {1+\epsilon ^{18/5}H^{\prime 2}}{1-\epsilon ^{18/5}H^{\prime 2}}, \end{align}
\begin{align} 4\epsilon ^{18/5}U_XH' &= (U_Z + \epsilon ^{18/5}W_X)(1-\epsilon ^{18/5}H^{\prime 2}), \end{align}
and at
$Z=0$
,
\begin{align} U_Z = W &= 0. \end{align}
We follow the standard methodology of Howell (Reference Howell1996) to derive the leading-order extensional flow model from the governing equations (3.32)–(3.34) that hold in this region. We find from (3.32b
) and (3.34) that the leading-order velocity
$U$
is independent of
$Z$
, hence the net mass conservation law (2.13) gives
\begin{align} U(X)&=\frac {1}{H(X)}. \end{align}
As a result, (3.32a
) implies that
$W_Z$
is independent of
$Z$
. Hence (3.32c
) gives that the leading-order pressure
$P$
is independent of
$Z$
, thus by (3.33b
) and (3.35) we have
\begin{align} P(X)&=-3H''(X)+\frac {2H'(X)}{H(X)^2}. \end{align}
Substituting the scaled variables (3.28) and symmetry condition (3.34) into the axial force balance (2.14) gives us
\begin{equation} \frac {\textrm {d}}{\textrm {d}{X}}\int _0^{H(X)}(-P+2U_X)\,\textrm {d}{Z} + \frac {H'H''}{(1+\epsilon ^{18/5}H'^2)^{3/2}}=-3H, \end{equation}
and upon substituting in (3.35) and (3.36), we obtain the leading-order governing equation
\begin{align} -4\left (\frac {{H}^{\prime}}{H}\right )' + 3\textit{HH}^{\prime\prime\prime} = -3H \end{align}
for the film thickness
$H(X)$
in this region. This is a typical equation for extensional flow within viscous sheets (Howell Reference Howell1994), with the three terms successively representing viscous effects, surface tension and gravity. We also note that the expression (2.15) for the tension transforms to
\begin{align} T &= \frac {1}{\epsilon ^3}\left (\epsilon ^{18/5}\int _0^{H(X)}(-P+2U_X)\,\textrm {d}{Z} + \frac {3}{( 1+\epsilon ^{18/5}H'^2)^{1/2}} \right )\nonumber \end{align}
\begin{align} &\sim \frac {1}{\epsilon ^3}\left (3 + \epsilon ^{18/5}T_1 + \cdots \right ), \end{align}
where
\begin{align} T_1=-\frac {4H'}{H}+3\textit{HH}''-\frac {3}{2}\,{H'}^2 \end{align}
is the first-order correction to the tension, comprised of the integrated viscous stress, the Laplace pressure due to the curvature of the free surface (again integrated across the film), and a contribution from the longitudinal component of the surface tension acting at the interface. The governing equation (3.38) may then be expressed as the axial force balance
\begin{equation} T_1'=-3H. \end{equation}
We now analyse the possible asymptotic behaviour of
$H(X)$
as
$X\rightarrow \pm \infty$
. We note that the governing equation (3.38) is translation-invariant. Once the arbitrary translation has been removed, there is a two-parameter family of solutions with the required behaviour to match with the solution (3.30) in the capillary–gravity region, namely
\begin{equation} H(X)\sim -\frac {X^3}{6}+c_1 X+c_0+O\big(X^{-2}\big) \quad \text{as }X\rightarrow -\infty , \end{equation}
where the constants
$c_0$
and
$c_1$
are arbitrary. The arbitrary translation has been exploited to remove the quadratic term in (3.42), ensuring that the matching conditions induce corrections of order
$\epsilon ^{4/5}$
(rather than
$\epsilon ^{2/5}$
) in
$\tilde {h}$
in the capillary–gravity region of § 3.3. The behaviour (3.42) corresponds to a dominant balance of the capillary and gravity terms in (3.38).
By searching for self-consistent dominant balances in (3.38), we determine three possible behaviours at the other end of the domain. The first generic behaviour results from a balance between the viscous and capillary terms, such that the solution terminates at a finite value of
$X$
(say
$a$
), with
\begin{equation} H(X)\sim \frac {4}{3}\sqrt {a-X}\quad \text{as }X\nearrow a. \end{equation}
The other generic behaviour, wherein viscous effects dominate, is for the solution to decay exponentially, such that
\begin{equation} H(X)\sim \mathrm{e}^{-b (X-\textit{a})}\quad \text{as }X\rightarrow \infty , \end{equation}
where
$a$
represents an arbitrary translation, and
$b$
is a positive constant. Finally, there is a non-generic far-field behaviour, due to a balance betweeen the viscous and gravity terms in (3.38), in which
$H(X)$
decays algebraically, with
\begin{equation} H(X)\sim \frac {8}{3(X-a)^2}-\frac {256}{21(X-a)^7}+O\big((X-a)^{-12}\big) \quad \text{as }X\rightarrow \infty , \end{equation}
where
$a$
once again represents an arbitrary translation.
Now let us examine the physical implications of each of the possible behaviours presented in (3.43). The first possibility (3.43a
) corresponds to the sheet pinching off at a finite
$X$
. By net mass conservation, the sheet velocity would need to simultaneously tend to infinity such that a unit net flux of liquid exits through the singularity at
$X=a$
. We reject this outcome as physically implausible. For the second possibility (3.43b
), we see that
$H(X)$
tends to zero as
$X\rightarrow \infty$
, but the correction to the tension
$T_1$
given by (3.40) tends to a finite value
$4b$
. In this scenario, a non-zero tension is being applied in the far field, for example by pulling or rolling the sheet, in addition to the downward gravitational force. Only for the final possibility (3.43c
) does
$T_1$
tend to zero, as well as the sheet thickness, as
$X\rightarrow \infty$
.
Henceforth, we focus on the case (3.43c
), which corresponds to the sheet falling freely under gravity without the application of any additional tension. The corresponding unique solution of (3.38) is found numerically by shooting from
$\infty$
, imposing the far-field behaviour (3.43c
) with
$a=0$
. The required translation
$a\approx -2.45$
to produce the behaviour (3.42) at the other end of the domain is then determined a posteriori. From this computed solution, we can read off the values of the constants in (3.42), namely
$c_1\approx -0.211$
and
$c_0\approx 0.931$
. We present the numerical solution to (3.38) subject to (3.43c
) in figure 5. The asymptotic behaviour (3.42) is shown as a green dashed curve, and the first term in (3.43c
) is shown as a red dashed curve, both showing excellent agreement with the computed numerical solution.
We conclude that by focusing on the case of a sheet stretching purely under gravity, we select a unique solution to the problem, and the corresponding downstream behaviour (2.18) anticipated in § 2 is found by translating (3.43c
) back into unscaled variables. We note that the alternative far-field behaviour (3.43b
) could be used instead to model an externally applied tension, thereby introducing an additional control parameter and an additional degree of freedom. This would only affect the leading-order solution in this final region, leaving the upstream flow (for
$\tilde {x}\lt \tilde {x}_0$
) unaffected.
4. Uniformly valid model for film thickness
We combine the governing equations (3.7), (3.25) and (3.38) from the lubrication flow, capillary–gravity and capillary–viscous–gravity regions to find an overall equation for the film thickness, namely
\begin{equation} \tilde {h}+\tilde {h} \tilde {h}''' = \begin{cases} \dfrac {1}{\tilde {h}^2} ,\,& \tilde {x}\lt 0, \\[10pt]\displaystyle \frac {4}{3}\epsilon ^2\left (\frac {\tilde {h}'}{\tilde {h}}\right )^{\prime} ,\,& \tilde {x}\gt 0, \end{cases} \end{equation}
which must be solved subject to the far-field boundary condition (see (2.18))
\begin{align} \tilde {h}(\tilde {x}) \sim \frac {8}{3 (\tilde {x}/\epsilon -x_*)^2} \quad \text{as } \tilde {x} \rightarrow \infty , \end{align}
where
$x_*$
is determined as part of the solution to ensure that
$\tilde {h} \rightarrow 1$
as
$\tilde {x} \rightarrow -\infty$
. Equation (4.1) is valid at leading order across the lubrication and extensional flow regions, and also captures the leading-order contribution from the Stokes flow region in the sense that
$\tilde {h}$
and its first two derivatives are continuous across
$\tilde {x}=0$
.
We determine the solution to (4.1) and (4.2) numerically by imposing that
$\tilde {h}$
and its first two derivatives satisfy (4.2) at
$\tilde {x}=L\gg 1$
; in practice, we find that
$L=10$
suffices. For each value of
$\epsilon$
, we obtain
$x_*$
via a shooting method, varying its value until we get the desired upstream behaviour
$\tilde {h} \rightarrow 1$
as
$\tilde {x}\rightarrow -\infty$
. From the analysis carried out above, we infer the asymptotic behaviour
\begin{equation} x_*\sim \epsilon ^{-1}\big(\tilde {x}_0 +\epsilon ^{2/5}a +\cdots \big) \quad \text{as }\epsilon \rightarrow 0, \end{equation}
where
$\tilde {x}_0\approx 1.3389$
and
$a\approx -2.45$
were obtained in § 3.1 and § 3.4, respectively. The dependence of
$x_*$
on
$\epsilon$
is shown in figure 6, and we see that the two-term asymptotic expansion (4.3) indeed gives the correct behaviour of
$x_*$
as
$\epsilon \rightarrow 0$
.
The value of
$\epsilon x_*$
such that the solution to (4.1), subject to (4.2) imposed at
$\tilde {x} = 10$
, satisfies
$\tilde {h} \rightarrow 1$
as
$\tilde {x} \rightarrow -\infty$
, plotted versus
$\epsilon ^{2/5}$
. The two-term asymptotic expansion of
$x_*$
given by (4.3) is shown as a dashed line.
We present numerical solutions for
$\tilde {h}$
, with various values of
$\epsilon$
, in figure 7. We see that as
$\epsilon$
increases, the qualitative behaviour of the solution remains unchanged. However, the size of the oscillations for
$\tilde {x}\lt 0$
decreases as
$\epsilon$
increases, hence the largest amplitude oscillations in the film thickness are observed as
$\epsilon \rightarrow 0$
. By considering the leading-order solution in the lubrication region, shown in figure 4, we deduce that the maximum film thickness is
$\tilde {h}\approx 1.278$
, which occurs at
$\tilde {x}\approx -1.658$
. We also see that the rate of decay of
$\tilde {h}$
as
$\tilde {x} \rightarrow \infty$
decreases as
$\epsilon$
increases.
5. Conclusions
In this paper, we address the question of how a thin viscous layer of glass transitions from lubrication flow to extensional flow at the apex of an overflow fusion trough in the limit of zero apex angle. In this limit, the problem consists of a thin film of fluid, with specified input flux
$Q$
and resulting thickness
\begin{equation} h_{{s}}=\left (\frac {3\mu Q}{\rho g}\right )^{1/3}, \end{equation}
flowing down a vertical wall under gravity before encountering a change in boundary condition from no-slip to symmetry at the apex of the trough. The resulting mathematical model has a single dimensionless parameter
\begin{equation} \textit{Ca}=\left (\frac {\rho g\mu ^2Q^2}{3\gamma ^3}\right )^{1/3}, \end{equation}
which describes the importance of viscous and gravitational effects relative to surface tension. We derive the leading-order solution to the dimensionless model in the small-
$\textit{Ca}$
limit by decomposing the domain into four regions. Beginning with the upstream lubrication flow region, as the flow propagates downstream, it moves through a region of fully two-dimensional Stokes flow, and then two regions of extensional flow in which different forces balance and control the flow. We derive the leading-order problem in each region, and determine the solution by matching with the flow in the region immediately upstream. We combine these models to arrive at a simple governing equation that captures the leading-order behaviour in each region. We solve this uniformly valid model numerically to determine the shape of the free surface across the whole flow domain.
Using the typical parameter values from table 1, we can determine approximate relationships between
$Q$
,
$h_{{s}}$
and
$\textit{Ca}$
, namely
\begin{equation} Q \approx 6 \times 10^{-7}\, \textit{Ca}^{3/2} \ \text{m}^2 \ \text{s}^{-1} , \quad h_{{s}} \approx 6 \times 10^{-3}\, \textit{Ca}^{1/2} \ \text{m}, \end{equation}
which give a sense of the range of values of these physical parameters for which the small-
$\textit{Ca}$
results are valid. We note that these relationships have been obtained using the assumption of constant viscosity, when in reality the viscosity is highly temperature dependent and can vary over multiple orders of magnitude. Since we focus on the model problem where
$Q$
is the only control parameter, it is of interest how the sheet thickness and length scale in each region of the flow scale with
$Q$
.
As shown in figure 7, the effects of the change in boundary condition at the sheet’s centreline are felt upstream of the apex, resulting in free-surface oscillations that decay over a length scale of order
$\textit{Ca}^{-1/3}\,h_{{s}} \propto Q^{1/9}$
. We find that the largest of these oscillations reaches dimensional height
$\approx 1.28 h_{{s}}$
at distance
$\approx 1.66(3\textit{Ca})^{-1/3}h_{{s}}$
upstream of the apex. Immediately downstream of the apex, the flow enters a regime of extensional flow wherein capillary forces balance gravity. The length scale of the flow in this region is once again
$\textit{Ca}^{-1/3}\,h_{{s}}$
, and we find that the leading-order free surface is given by a cubic that appears to pinch off at a distance
$\approx 1.34(3\textit{Ca})^{-1/3}h_{{s}}$
from the apex of the trough. However, as the height of the free surface approaches zero, the flow enters a new regime wherein surface tension, gravity and viscous forces all balance. Analysis of the capillary–viscous–gravity region reveals that the length scale of this regime is of order
$\textit{Ca}^{-1/5}\,h_{{s}}\propto Q^{1/5}$
, and that the downstream thickness of the fluid sheet is smaller than the upstream film thickness
$h_{{s}}$
by a factor of order
$\textit{Ca}^{2/5}$
. We hence conclude that the film thickness far downstream scales with
$Q^{3/5}$
, as opposed to the upstream scaling
$Q^{1/3}$
, in the small-
$\textit{Ca}$
limit.
We focus on the particular case where the downstream fluid sheet falls freely under gravity, with no additional applied tension. However, we also find that any imposed downstream tension only influences the leading-order flow in this final region, and everything upstream of the apparent pinch point at the bottom of the capillary–gravity region is unaffected. By considering the asymptotic behaviour of the free surface in this region, we determine the corresponding downstream behaviour, which would be essential to determine a numerical solution to the full dimensionless problem (2.9)–(2.12) and (2.17). Although it is outside the scope of this paper, we briefly outline below how such a numerical solution could be obtained.
We first note that the geometry of the problem can be simplified by scaling out the domain height
$h(x)$
(see e.g. O’Kiely et al. Reference O’Kiely, Breward, Griffiths, Howell and Lange2015), at the expense of complicating the governing equations. The infinite domain of this problem should be truncated to a large, finite domain (say
$x \in [-L,L],\ L \gg 1$
), with the upstream velocity and film thickness (2.12) imposed at
$x = -L$
. To model the sheet falling under gravity, with no additional tension applied, we suggest imposing
$\boldsymbol{\sigma \boldsymbol{\cdot }i} = \boldsymbol{0}$
at
$x=L$
, while any additional applied tension can be included by instead using
$\boldsymbol{\sigma \boldsymbol{\cdot }i} = T_1\boldsymbol{i}/h$
at
$x=L$
. It should be noted, however, that any such numerical solution is likely to be tricky when
$\textit{Ca}$
is small, due to the large disparity in length scales exhibited by the solution.
At the apex of the trough, the flow changes from lubrication flow to extensional flow across a region of
$O(1)$
aspect ratio. In this region, the flow is governed by the full two-dimensional Stokes equations, with a free surface that is flat to lowest order. As shown in Appendix B, the resulting problem can be solved using the Wiener–Hopf method. The flow in this region bridges the transition from a no-slip to symmetry condition at the centreline, and smooths out a jump in the third derivative of the film height exhibited by the outer solution. However, we find that in the limit of small
$\textit{Ca}$
considered here, these details do not exert any global influence on the leading-order solution outside the Stokes flow region.
There are many areas of further research that would help us to better understand the overflow fusion process. Since the angle at the bottom of the trough is likely to be
$O(1)$
in practice, the first avenue is to consider a non-zero wedge angle. We anticipate that this change would affect how the flow transitions between the lubrication and extensional regimes, but that the flow far downstream should remain unaffected. It would also be of interest to reintroduce time dependence into the model, and analyse the stability of the steady state that we have found, which is a key concern for Corning Inc. (Bohun et al. Reference Bohun, Breward, Cummings and Witelski2010). The reintroduction of time dependence would also allow us to explore the start-up process. It would also be of interest to include spanwise perturbations to the flow (i.e. into the page in figure 2) and thus analyse the possibility of ribbing or fingering instabilities (Troian et al. Reference Troian, Herbolzheimer, Safran and Joanny1989), which again are of industrial relevance.
Since glass viscosity varies considerably with temperature, future refinements of our model should incorporate heat transfer. Introducing temperature dependence may also result in non-uniform surface tension and hence Marangoni stresses at the free surface, which can affect the formation of instabilities at the advancing front (Kalliadasis, Kiyashko & Demekhin Reference Kalliadasis, Kiyashko and Demekhin2003) or lead to shock formation (Münch Reference Münch2000). It would be enlightening to explore the interplay between a temperature-dependent viscosity and Marangoni stresses, since previous studies of thin film flows with Maragoni effects often neglect this aspect (Bohun et al. Reference Bohun, Breward, Cummings and Witelski2010). Finally, it would be useful to compare the results from our model with a numerical solution of the full dimensionless problem (2.9)–(2.12), as well as with experimental data.
The model derived in this paper provides an important stepping stone towards understanding the transition from lubrication flow to extensional flow during overflow fusion, and thus for optimising this glass manufacturing process. The model provides a framework upon which additional modelling assumptions can be added to further develop mathematical and industrial understanding of the overflow fusion process.
Acknowledgements
We are grateful to J. Abbott (Corning Inc.) and I.M. Griffiths (University of Oxford) for very useful discussions.
Funding
O.W. Hart acknowledges Pembroke College and the Mathematical Institute, University of Oxford, for providing financial support.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Axial force balance derivation
We derive the axial force balance (2.14) using the governing equations (2.9)–(2.11). We begin by using (2.9a ) to write (2.9b ) in the form
\begin{equation} \left (-p+2 u_x\right )_x+\left (u_z+w_x\right )_z+3=0. \end{equation}
We then integrate with respect to
$z$
, and apply the boundary conditions (2.10c
) and (2.11) to obtain
\begin{align} \int _0^{h(x)}\left (-p+2u_{x}\right )_x\,\textrm {d}{z} = -\left .\frac {4u_xh_x}{1-h_x}\right |_{z=h}+u_z|_{z=0}-3h. \end{align}
The Leibniz rule may be used to transform (A2) to
\begin{align} \frac {\textrm {d}}{\textrm {d}{x}}\int _0^{h(x)}(-p+2u_x)\,\textrm {d}{z} = h_x(-p+2u_x)|_{z=h}-\left .\frac {4u_xh_x}{1-h_x}\right |_{z=h}+u_z|_{z=0}-3h, \end{align}
and upon substituting the boundary condition (2.10b ), we recover
\begin{align} \frac {\textrm {d}}{\textrm {d}{x}}\left (\int _0^{h(x)}(-p+2u_x)\,\textrm {d}{z} \right )+ \frac {1}{\textit{Ca}} \frac {h_xh_{\textit{xx}}}{(1+h_x^2)^{3/2}} =u_z|_{z=0}-3h, \end{align}
which, using the definition (2.15) of the net tension in the sheet, can be rewritten to give (2.14).
Appendix B. Solution to the Stokes flow problem
We consider the Stokes flow problem (3.17)–(3.19). We begin by defining new independent variables
$(\xi ,\zeta )=(x,z)/h_0$
, to scale out the constant domain height
$h_0$
, and introducing a streamfunction
$\psi$
such that
$h_0(u_0,w_0)=(\psi _\zeta ,-\psi _\xi )$
, to satisfy the incompressibility condition (3.17a
). By eliminating
$\mathfrak{p}_1$
from (3.17b
) and (3.17c
), we find that
$\psi$
satisfies the two-dimensional biharmonic equation, and we obtain a parameter-free boundary-value problem for
$\psi$
, namely
\begin{equation} {\nabla} ^4 \psi = 0, \end{equation}
for
$(\xi ,\zeta )\in \mathbb{R}\times [0,1]$
, with boundary conditions at
$\zeta =1$
,
\begin{align} \psi = 1, \quad \psi _{\zeta \zeta } = 0, \end{align}
and at
$\zeta =0$
,
\begin{align} \psi = \psi _\zeta =0 \quad \text{for } \xi \lt 0, \end{align}
\begin{align} \psi = \psi _{\zeta \zeta } =0 \quad \text{for } \xi \gt 0. \end{align}
We solve the problem (B1)–(B2) for
$\psi$
using the Wiener–Hopf method (Carrier, Krook & Pearson Reference Carrier, Krook and Pearson2005), following the arguments laid out in Luchini (Reference Luchini1991). We note that a detailed analytic solution to this problem is presented in Barcilon & MacAyeal (Reference Barcilon and MacAyeal1993), and related problems are explored in Richardson (Reference Richardson1970) and Hutter & Olunloyo (Reference Hutter and Olunloyo1980).
Defining the Fourier transform via
\begin{align} \hat {f}(k)=\int _{-\infty }^\infty f(\xi )\,\mathrm{e}^{-\mathrm{i}k\xi }\, \textrm {d}{\xi }, \end{align}
and defining
$\omega = {\nabla} ^2 \psi$
so that
${\nabla} ^2 \omega = 0$
, taking Fourier transforms gives
\begin{align} \hat {\psi }_{\zeta \zeta } -k^2\hat {\psi }&=\hat {\omega }, \end{align}
\begin{align} \hat {\omega }_{\zeta \zeta }-k^2\hat {\omega } &=0, \end{align}
with boundary conditions at
$\zeta =1$
,
\begin{align} \hat {\psi } = 2\unicode{x03C0}\, \delta (k), \quad \hat {\omega } = 0, \end{align}
and at
$\zeta = 0$
,
\begin{align} \hat {\psi }=0, \quad \hat {\omega }=\hat {\omega }_*(k), \end{align}
where
$\delta (k)$
is the Dirac delta function, and
$\omega _*(\xi ) = \omega (\xi ,0)$
is a priori unknown, but zero for
$\xi \gt 0$
, thus
$\hat {\omega }_*(k)$
is holomorphic for
$\textrm {Im}\,{k} \gt 0$
.
Solving (B5)–(B7), we find that
\begin{align} \hat {\psi }(k,\zeta )=\frac {(2-\zeta )\sinh (k\zeta )-\zeta \sinh ((2-\zeta )k)}{4k\sinh ^2k}\,\hat {\omega }_*(k)+2\unicode{x03C0} \zeta\, \delta (k), \end{align}
hence
\begin{align} \hat {u}_*(k) = -K(k)\,\hat {\omega }_*(k) + F(k), \end{align}
where
$u_*(\xi ) = \psi _\zeta (\xi ,0)$
is zero for
$\xi \lt 0$
, thus
$\hat {u}_*(k)$
is holomorphic for
$\textrm {Im}\, k \lt 0$
, and
\begin{align} K(k) = \frac {\sinh {k}\cosh {k}-k}{2k\sinh ^2k}, \quad F(k) = 2\unicode{x03C0}\, \delta (k). \end{align}
We include a minus sign in (B9) when defining
$K(k)$
for consistency with Luchini (Reference Luchini1991) and Barcilon & MacAyeal (Reference Barcilon and MacAyeal1993), and so that
$K(k)$
is positive for real
$k$
.
We seek to determine
$\hat {\omega }_*$
by considering (B9) and exploiting the holomorphicity properties of
$\hat {\omega }_*$
and
$\hat {u}_*$
. Given product and sum decompositions
\begin{align} K(k)&=K_+(k)\, K_-(k), \end{align}
\begin{align} F(k)/K_-(k) &= F_+(k)+F_-(k), \end{align}
where subscript plus or minus denotes holomorphicity for
$\textrm {Im}\,{k} \gt 0$
or
$\textrm {Im}\,{k} \lt 0$
, respectively, (B9) becomes
\begin{align} \hat {u}_*(k)/K_-(k)-F_-(k) = -K_+(k)\,\hat {\omega }_*(k) +F_+(k), \end{align}
where the left-hand side is holomorphic for
$\textrm {Im}\, k \lt 0$
, and the right-hand side is holomorphic for
$\textrm {Im}\, k \gt 0$
. By the standard Wiener–Hopf argument, both sides must be equal to zero. This argument determines
$\hat {\omega }_*$
and thus, via an inverse Fourier transform of (B8),
$\psi$
.
For a given function
$G(k)$
, one possible sum decomposition is given by (Luchini Reference Luchini1991; Carrier et al. Reference Carrier, Krook and Pearson2005)
\begin{align} G_\pm (k) = \mp \frac {1}{2\unicode{x03C0} \mathrm{i}} \int _{-\infty }^\infty \frac {G(\xi )}{k-\xi \pm \mathrm{i}0}\, \textrm {d}{\xi } = \mp \frac {1}{2\unicode{x03C0} \mathrm{i}}\, G(\xi ) \ast \frac {1}{\xi \pm \mathrm{i}0}, \end{align}
where
$\ast$
denotes convolution. This formulation allows us to determine
$G_\pm (k)$
using the convolution theorem and the fact that the Fourier transform of the distribution
$(\xi \pm \mathrm{i}0)^{-1}$
is
$\mp 2 \unicode{x03C0} \mathrm{i}\,H(\pm k)$
, where
$H(x)$
is the Heaviside step function. Due to the properties of the Dirac delta function, (B13) immediately gives
\begin{align} F_\pm (k)=\mp \frac {1}{\mathrm{i}K_-(0)} \frac {1}{k\pm \mathrm{i}0}, \end{align}
and we determine the product decomposition for
$K(k)$
by performing a sum decomposition on
$\log {K(k)}$
and then taking exponentials. The resulting Fourier transforms cannot be determined analytically, but may be found numerically using, for instance, a fast Fourier transform. We omit the computational details here, as they can be found in Luchini (Reference Luchini1991).
Using the decomposition (B13), we have by symmetry that
\begin{align} K_+(0)=K_-(0)=\frac {1}{\sqrt {3}}, \end{align}
and once we obtain
$K_+(k)$
, we determine
$\psi$
using the convolution theorem and Fourier transform of
$(\xi + {\rm i}0)^{-1}$
, namely
\begin{align} \psi (\xi ,\zeta )= \frac {\sqrt {3}}{2\unicode{x03C0} }\int _{-\infty }^{-\xi }\hat {B}(k,\zeta )\,\textrm {d}{k} + \zeta , \end{align}
where
\begin{align} B(\xi ,\zeta )=\frac {1}{K_+(\xi )}\,\frac {(2-\zeta )\sinh (\xi \zeta )-\zeta \sinh ((2-\zeta )\xi )}{4\xi \sinh ^2\xi }. \end{align}
Numerical solution to (4.1) with far-field condition (4.2) imposed at
$\tilde {x}=10$
, choosing
$x_*$
such that
$\tilde {h} \rightarrow 1$
as
$\tilde {x} \rightarrow -\infty$
, for
$\epsilon =0.01$
(green),
$\epsilon =0.1$
(red) and
$\epsilon =0.5$
(blue). The composite leading-order solution (black, dashed) is found by combining (3.30) with the numerical solution to (3.7) subject to boundary conditions (3.31), where
$\tilde {x}_0\approx 1.3389$
.
A density plot of the value of
$\psi (\xi ,\zeta )$
. Streamlines of the flow, given by constant
$\psi$
, are indicated on the plot, with flux
$0.1$
passing between each pair of neighbouring streamlines.
The scaled velocity component
$\psi _\zeta = h_0u_0$
across the thickness of the film, for different fixed values of
$\xi$
(labelled on the plot), namely,
$-2$
(blue),
$-0.1$
(red),
$0.1$
(green),
$0.5$
(orange), and
$2$
(black).
Upon determining
$\psi$
, we are able to visualise the streamlines of the flow, shown in figure 8, and recover the scaled velocity field
$(\psi _\zeta ,-\psi _\xi )=h_0(u_0,w_0)$
. The scaled velocity component
$\psi _\zeta = h_0u_0$
, for a range of values of
$\xi$
, is shown in figure 9. Of particular interest is the change in the velocity between
$\xi =-2$
and
$\xi = -0.1$
, showing that the change in boundary conditions at the origin has an effect on the flow upstream of
$\xi =0$
. These figures demonstrate how the velocity profile transitions from parabolic to uniform, and in particular we find that the scaled velocity field exhibits upstream behaviour
\begin{align} h_0u_0\rightarrow 3\left (\zeta -\frac {\zeta ^2}{2}\right ),\quad h_0w_0\rightarrow 0\quad \text{as }\xi \rightarrow -\infty , \end{align}
and downstream behaviour
\begin{align} h_0u_0\rightarrow 1, \quad h_0w_0\rightarrow 0 \quad \text{as }\xi \rightarrow \infty . \end{align}
Returning to the original variables used in § 3.2 confirms that the flow in the Stokes flow region satisfies the upstream matching condition (3.20), and verifies our assertion (3.21) about the downstream flow behaviour.
Appendix C. Capillary–viscous–gravity scalings
We derive the scalings (3.28) in the capillary–viscous–gravity extensional flow region by considering the governing equations (3.2) and boundary conditions (3.3) that hold in the capillary–gravity extensional flow region. We begin by scaling
\begin{align} \tilde {x}-\tilde {x}_0=\epsilon ^mX, \quad \tilde {z} = \epsilon ^n Z, \quad \tilde {h} = \epsilon ^nH, \end{align}
for some positive exponents
$m$
and
$n$
, which are to be determined. By considering the net conservation of mass (2.13), seeking a dominant balance in (3.2a
), and balancing the pressure gradient with the gravity term in (3.2b
), we establish the scalings
\begin{align} \tilde {u}=\epsilon ^{-n} U, \quad \tilde {w} = \epsilon ^{-m} W, \quad \tilde {p}-\tilde {p}(\tilde {x}_0) = \epsilon ^mP. \end{align}
As a result, (3.2b ) transforms to
\begin{align} P_X &= \epsilon ^{2-n-2m}U_{\textit{XX}} + \epsilon ^{-3n}U_{\textit{ZZ}} + 3, \end{align}
and (3.3b
) transforms to, at
$Z=H(X)$
,
\begin{align} \frac {3\epsilon ^{n-3m}H^{\prime\prime}}{(1+\epsilon ^{2+2n-2m}H^{\prime 2})^{3/2}}&= -\epsilon ^{-m}\,\tilde {p}(\tilde {x_0})-P-2\epsilon ^{2-n-2m}U_X\frac {1+\epsilon ^{2+2n-2m}H^{\prime 2}}{1-\epsilon ^{2+2n-2m}H^{\prime 2}}. \end{align}
By balancing the extensional viscous term
$U_{\textit{XX}}$
with the gravitational force in (C3), and balancing the surface tension term in (C4) with the pressure at leading order, we find that
$\tilde {p}(\tilde {x}_0)$
must be zero, and that
\begin{align} 2-n-2m=0, \quad n-3m=0, \end{align}
from which we recover
\begin{align} m=\frac {2}{5}, \quad n=\frac {6}{5}, \end{align}
which gives the scalings (3.28).
Open access
