1. Introduction
Imagine a liquid sheet, or curtain ejected downwards or obliquely from a long horizontal slot, and let it fall down due to gravity. The first study of curtains was conducted more than sixty years ago: G. I. Taylor (whose work was described in an appendix of the paper by Brown Reference Brown1961) used the slim-jet approximation to calculate the characteristics of downwards-ejected curtains.
Numerous papers have been published since then, of which the following are relevant to the present work: Finnicum, Weinstein & Ruschak (Reference Finnicum, Weinstein and Ruschak1993) and Lhuissier, Brunet & Dorbolo (Reference Lhuissier, Brunet and Dorbolo2016) examined steady oblique curtains affected by an external pressure applied on one side; Benilov (Reference Benilov2019) examined evolving oblique curtains without external pressure; Roche et al. (Reference Roche, Grand, Brunet, Lebon and Limat2006), Dyson et al. (Reference Dyson, Brander, Breward and Howell2009) and Benilov, Barros & O’Brien (Reference Benilov, Barros and O’Brien2016) examined the stability of downwards-ejected curtains; and Benilov (Reference Benilov2023) examined that of oblique ones.
It is now well understood that the dynamics of curtains depends strongly on the Weber number,
${\textit{We}}$
, characterising the strength of inertia relative to surface tension. In particular, it was shown (Keller & Weitz Reference Keller and Weitz1957; Benilov Reference Benilov2019) that, if
${\textit{We}}\lt 1$
, the asymptotic equations for slim curtains admit paradoxical solutions bending upwards – i.e. against gravity. These solutions were discussed further by Weinstein et al. (Reference Weinstein, Ross, Ruschak and Barlow2019), Benilov (Reference Benilov2023) and Della Pia et al. (Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023), and a consensus has emerged that such solutions cannot be observed experimentally due to their instability.
The present paper examines how sufficiently strong viscosity alters the streamwise dynamics of liquid curtains (cross-stream relaxation of viscous stress and similar effects are already well understood).
The paper has the following structure. Section 2 presents the asymptotic equations describing slim curtains (their derivation is presented in Appendix A). Sections 3 and 4 describe the steady and evolving solutions of these equations, respectively; and § 5 outlines how the results of this paper can be tested experimentally and how viscosity modifies the dynamics in the range
${\textit{We}}\lt 1$
(where the unstable upward-bending curtains were thought to exist).
2. Formulation
2.1. Non-dimensional parameters
Let an incompressible liquid of density
$\rho$
, kinematic viscosity
$\nu$
and surface tension
$\sigma$
be ejected from a horizontal, infinitely long outlet of width
$2H$
– see figure 1. Introducing a velocity scale
$U$
and the acceleration due to gravity
$g$
, one can define a streamwise spatial scale
\begin{align} L=\frac {U^{2}}{g}, \end{align}
and three non-dimensional parameters
\begin{equation} \varepsilon =\frac {H}{L},\quad \gamma =\frac {\sigma }{\rho {\textit{HU}}^{2}},\quad \mu =\frac {\nu }{\textit{UL}}, \end{equation}
which are related to the traditional dimensionless groups by
\begin{align} \varepsilon =\frac {1}{{\textit{Fr}}^{2}},\quad \gamma =\frac {1}{\textit{We}},\quad \mu ={\textit{Fr}}^{2}{\textit{Re}}, \end{align}
where
${\textit{Fr}}=U/(gH)^{1/2}$
is the Froude number,
${\textit{We}}=\rho {\textit{HU}}^{2}/\sigma$
the Weber number and
${\textit{Re}}=UH/\nu$
the Reynolds number.
Figure 1. The setting: a two-dimensional liquid curtain ejected from an outlet of width
$2H$
, tilted at an angle
$\alpha _{0}$
;
$(x,z)$
are Cartesian coordinates, and
$(l,\tau )$
are the curvilinear coordinates associated with the curtain’s centreline.
In this paper, it is assumed that
\begin{equation} \varepsilon \ll 1,\quad \gamma \sim 1,\quad \mu \sim 1. \end{equation}
The first of these assumptions has been used in most, if not all, theoretical studies of oblique curtains: it implies that
$L\gg H$
– hence, the curtain can be regarded as slim. The second assumption has also been widely used because it includes the most interesting particular case of critical curtains
$\gamma =1$
(more details given later).
The third assumption, however, differs from the assumption
$\mu \ll 1$
used by previous researchers. As the present results show, stronger viscosity dramatically alters the dynamics of liquid curtains.
2.2. Governing equations and initial conditions
Assumptions (2.4) allow one to reduce the Navier–Stokes equations to a much simpler asymptotic set. The underlying asymptotic analysis can be found in Appendix A, while this section presents its final result.
Liquid curtains can be conveniently modelled using curvilinear coordinates
$ ( l,\tau )$
, where the streamwise variable (arc length)
$l$
and the cross-stream variable
$\tau$
form an orthogonal set (see Benilov Reference Benilov2019, and also Appendix A). Then, the curtain’s centreline corresponds to
$\tau =0$
and its boundaries, to
$\tau =(1/2)W$
where
$W$
is the curtain’s local width. Note that
$W$
depends on the arc length
$l$
and the time
$t$
, as do the Cartesian coordinates of the centreline
$ ( x,z )$
and the local angle
$\alpha$
between the centreline and the horizontal.
As shown in Appendix A, strong viscosity eliminates cross-stream shear, so that the leading-order streamwise velocity
$u$
depends only on
$t$
and
$l$
. The slim-curtain approximation, in turn, makes the cross-stream velocity small.
The following non-dimensional variables will be used:
\begin{align}&\qquad\qquad\qquad\qquad l_{\textit{nd}}=\frac {l}{L},\quad t_{\textit{nd}}=\frac {U}{L}t, \end{align}
\begin{align}& W_{\textit{nd}}=\frac {W}{H},\quad u_{\textit{nd}}=\frac {u}{U},\quad \alpha _{\textit{nd}}=\alpha ,\quad x_{\textit{nd}}=\frac {x}{L},\quad z_{\textit{nd}}=\frac {z}{L}. \end{align}
The non-dimensional unknowns
$ ( W_{\textit{nd}},u_{\textit{nd}},\alpha _{\textit{nd}},x_{\textit{nd}} ,z_{\textit{nd}} )$
are governed by the following asymptotic equations (subscripts
$_{\textit{nd}}$
omitted):
\begin{align}&\qquad\qquad\qquad\qquad\qquad \frac {\partial W}{\partial t}+\dfrac {\partial \left ( uW\right ) }{\partial l}=0, \end{align}
\begin{align}&\quad \frac {\partial u}{\partial t}+u\frac {\partial u}{\partial l}=-\sin \alpha +\frac {4\mu }{W}\dfrac {\partial }{\partial l}\left ( W\frac {\partial u}{\partial l}\right ) -\frac {\partial ^{2}x}{\partial t^{2}}\cos \alpha -\frac {\partial ^{2}z}{\partial t^{2}}\sin \alpha , \end{align}
\begin{align}& 2u\frac {\partial \alpha }{\partial t}+\left ( u^{2}-\frac {2\gamma }{W}-4\mu \frac {\partial u}{\partial l}\right ) \frac {\partial \alpha }{\partial l}=-\cos \alpha +\frac {\partial ^{2}x}{\partial t^{2}}\sin \alpha -\frac {\partial ^{2}z}{\partial t^{2}}\cos \alpha , \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\quad \frac {\partial x}{\partial l}=\cos \alpha ,\quad \frac {\partial z}{\partial l}=\sin \alpha , \end{align}
where the coefficients
$\mu$
and
$\gamma$
are given by (2.2).
The following comments should clarify the physical meaning of (2.7)–(2.10).
-
(i) Equation (2.7) reflects conservation of mass.
-
(ii) Equation (2.8) governs the balance of streamwise forces. The first term on its right-hand side describes the effect of gravity, the second term describes viscosity and the last two terms describe the force of inertia arising due to the use of coordinates associated with the centreline (which, generally, moves with acceleration).
-
(iii) Equation (2.9) governs cross-stream forces. The terms on its right-hand side describe gravity and inertia; the first term on its left-hand side describes the Coriolis force (
$\partial \alpha /\partial t$
is the local angular velocity of the coordinate system’s rotation); and the three terms in parentheses describe the centrifugal force, cross-stream pressure gradient due to surface tension and viscosity. -
(iv) Equations (2.10) reflect geometric relationships between the parameters of the centreline.
It is instructive to consider the horizontal momentum balance, i.e. the following combinations of the governing equations:
\begin{align} \left (2.5\right ) \times u\cos \alpha +\left (2.6\right ) \times W\cos \alpha -\left (2.7\right ) \times W\sin \alpha , \end{align}
which yields, after straightforward algebra
\begin{equation} \frac {\partial \left ( \textit{Wu}\cos \alpha \right ) }{\partial t}+\frac {\partial }{\partial l}\left [ \left ( {\textit{Wu}}^{2}-2\gamma -4\mu W\frac {\partial u}{\partial l}\right ) \cos \alpha \right ] =-W\frac {\partial ^{2}x}{\partial t^{2}}+\textit{Wu}\frac {\partial \alpha }{\partial t}\sin \alpha . \end{equation}
The terms on the right-hand side of (2.12) represent the horizontal components of the force of inertia and Coriolis force. Note that, if (2.12) were written in terms of stationary (inertial) coordinates, its would not have a right-hand side.
2.3. Initial conditions
To set initial conditions for (2.7)–(2.10), one needs to determine their order with respect to time derivatives. To do so, observe that the terms
$\partial ^{2}x/\partial t^{2}$
and
$\partial ^{2}z/\partial t^{2}$
in (2.8)–(2.9) can be both expressed through
$\partial \alpha /\partial t$
and
$\partial ^{2}\alpha /\partial t^{2}$
. Indeed, using (2.10), one can show that
\begin{equation} \frac {\partial }{\partial l}\frac {\partial ^{2}x}{\partial t^{2}}=-\frac {\partial ^{2}\alpha }{\partial t^{2}}\sin \alpha -\left ( \frac {\partial \alpha }{\partial t}\right ) ^{2}\cos \alpha ,\quad \frac {\partial }{\partial l}\frac {\partial ^{2}z}{\partial t^{2}}=\frac {\partial ^{2}\alpha }{\partial t^{2}}\cos \alpha -\left ( \frac {\partial \alpha }{\partial t}\right ) ^{2}\sin \alpha . \end{equation}
Then, differentiating (2.8)–(2.9) with respect to
$l$
, rearranging the result using (2.13) and separating the terms involving
$\sin \alpha$
and
$\cos \alpha$
, one obtains
\begin{align} \frac {\partial }{\partial l}\left [ \!\frac {\partial u}{\partial t} +u\frac {\partial u}{\partial l}-\frac {4\mu }{W}\dfrac {\partial }{\partial l}\left (\! W\frac {\partial u}{\partial l}\!\right )\! \right ] -\left [\! 2u\frac {\partial \alpha }{\partial t}+\left (\! u^{2}-\frac {2\gamma } {W}-4\mu \frac {\partial u}{\partial l}\!\right )\! \frac {\partial \alpha }{\partial l}\!\right ] \frac {\partial \alpha }{\partial l}-\left ( \!\frac {\partial \alpha }{\partial t}\!\right ) ^{2} &=0, \end{align}
\begin{align} \frac {\partial ^{2}\alpha }{\partial t^{2}}+\left [ \!\frac {\partial u}{\partial t}+u\frac {\partial u}{\partial l}-\frac {4\mu }{W}\dfrac {\partial }{\partial l}\left (\! W\frac {\partial u}{\partial l}\!\right )\! \right ] \frac {\partial \alpha }{\partial l} +\frac {\partial }{\partial l}\left [\! 2u\frac {\partial \alpha }{\partial t}+\left (\! u^{2}-\frac {2\gamma }{W}-4\mu \frac {\partial u}{\partial l}\!\right ) \frac {\partial \alpha }{\partial l}\!\right ] &=0. \end{align}
Equations (2.7) and (2.14)–(2.15) form a closed set for
$W$
,
$u$
and
$\alpha$
– whereas
$x$
and
$z$
have been eliminated. Given that the set involves
$\partial W/\partial t$
,
$\partial u/\partial t$
and
$\partial ^{2}\alpha /\partial t^{2}$
, the following initial conditions are needed:
\begin{equation} W=W_{i}(l),\quad u=u_{i}(l),\quad \alpha =\alpha _{i}(l),\quad \frac {\partial \alpha }{\partial t}=\beta _{i}(l)\quad \text{at}\quad t=0, \end{equation}
where
$W_{i}(l)$
…
$\beta _{i}(l)$
are known functions. Physically, conditions (2.16) determine the curtain’s initial thickness, streamwise velocity, slope and local angular velocity.
Note that (2.7) and (2.14)–(2.15) are more convenient computationally than the original set (2.7)–(2.10) (with the latter,
$x$
and
$z$
are to be computed using (2.10), and it is unclear how to accurately compute
$\partial ^{2}x/\partial t^{2}$
and
$\partial ^{2}z/\partial t^{2}$
which appear on the right-hand sides of (2.8)–(2.9)). Theoretically, however, both sets will be used as convenient.
2.4. Boundary conditions
Governing equations (2.7)–(2.10) involve the second spatial derivative of
$u$
(two boundary conditions needed) and the first derivatives of the remaining variables (one condition per variable). It would thus appear that the following boundary conditions suffice:
\begin{align}&\quad\,\, W=2\quad \text{at}\quad l=0, \end{align}
\begin{align}&\quad\,\, \alpha =\alpha _{0}\quad \text{at}\quad l=0, \end{align}
\begin{align}& x=0,\quad z=0\quad \text{at}\quad l=0, \end{align}
\begin{align}&\quad\,\,\, u=1\quad \text{at}\quad l=0, \end{align}
\begin{align}&\, \left \vert \frac {\partial u}{\partial l}\right \vert \lt \infty \quad \text{as}\quad l\rightarrow \infty . \end{align}
Conditions (2.17) and (2.20) reflect the way
$W$
and
$u$
have been non-dimensionalised, (2.18) reflects the orientation of the nozzle and (2.19) implies that the origin of the Cartesian coordinate system is located at the outlet. Finally, condition (2.21) bounds the viscous stress at infinity.
Note also that the governing set is parabolic with respect to
$u$
, but hyperbolic with respect to the other variables. Thus, one should keep in mind that the causality principle allows one to set boundary conditions for all the hyperbolic variables at the same point (as in (2.17)–(2.20)) only if the velocities of all the characteristics are directed away from this point. If one or more are directed towards it, such boundary conditions amount, physically, to imposing restrictions on events that occurred in the past and, mathematically, to constraining the allowable initial conditions.
Thus, to justify conditions (2.17) and (2.20), one needs to calculate the velocities of the characteristics at
$l=0$
.
The velocity
$c_{1}$
of the characteristic associated with
$W$
are readily deduced from (2.7)
\begin{align} c_{1}=u. \end{align}
To determine the other characteristics, keep in (2.15) the second derivatives of
$\alpha$
and omit all the other terms
\begin{align} \frac {\partial ^{2}\alpha }{\partial t^{2}}+2u\frac {\partial ^{2}\alpha }{\partial l\partial t}+\left ( u^{2}-\frac {2\gamma }{W}-4\mu \frac {\partial u}{\partial l}\right ) \frac {\partial ^{2}\alpha }{\partial l^{2}}=0. \end{align}
Now, one can readily calculate the velocities of the characteristics associated with
$\alpha$
\begin{equation} c_{2}=u+\sqrt {\frac {2\gamma }{W}+4\mu \frac {\partial u}{\partial l}},\quad c_{3}=u-\sqrt {\frac {2\gamma }{W}+4\mu \frac {\partial u}{\partial l}}. \end{equation}
Physically,
$c_{1}$
represents propagation of disturbances by advection, whereas
$c_{2}$
and
$c_{3}$
describe viscosity-affected capillary disturbances propagating downstream and upstream, respectively.
Recalling boundary conditions (2.20)–(2.19), one can verify that
$c_{1}$
and
$c_{2}$
are positive – hence, directed away from the boundary, whereas
$c_{3}$
is positive if only if
\begin{equation} 1-4\mu \left ( \frac {\partial u}{\partial l}\right ) _{l=0}-\gamma \gt 0. \end{equation}
In what follows and until § 4.2, a stronger inequality will be assumed to hold
\begin{equation} u-4\mu \frac {\partial u}{\partial l}-\gamma \gt 0\quad \text{for}\quad l\geqslant 0. \end{equation}
This inequality guarantees that boundary conditions (2.20)–(2.19) apply, plus it removes the necessity for an extra boundary condition at infinity (which would be needed if
$c_{3}$
became negative as
$l\rightarrow$
$\infty$
).
The boundary conditions for the alternative set (2.7), (2.14 )–(2.15) are described in Appendix B.
3. Steady curtains
When studying time-independent solutions of (2.7 )–(2.10), it is convenient to replace (2.8) for
$u$
with the horizontal momentum (2.12). Omitting, thus, the time derivatives in (2.7) and (2.12), integrating them with respect to
$l$
and taking into account boundary conditions (2.20), one obtains
\begin{align}&\qquad\qquad W=\frac {2}{u},\\[-12pt]\nonumber \end{align}
\begin{align}& \left ( u-\gamma -\frac {4\mu }{u}\frac {\mathrm{d}u}{\mathrm{d}l}\right ) \cos \alpha =F, \end{align}
where the constant
$F$
is, physically, the horizontal momentum flux.
Omitting the time derivatives in (2.9) and using (3.1)–(3.2) to eliminate
$W$
and
$\mathrm{d}u/\mathrm{d}l$
, respectively, one obtains
\begin{align} \frac {\mathrm{d}\alpha }{\mathrm{d}l}=-\frac {\cos ^{2}\alpha }{Fu}.\\[-12pt]\nonumber \end{align}
In the limit
$\mu \rightarrow 0$
, (3.1)–(3.3) coincide with the steady-curtain equations derived by Finnicum et al. (Reference Finnicum, Weinstein and Ruschak1993) and Benilov (Reference Benilov2019). They also coincide with (9)–(10) of Kistler & Scriven (Reference Kistler and Scriven1994) if the latter are considered in the asymptotic limit used in the present paper – i.e.
${\textit{Re}}\ll 1$
,
${\textit{We}}\sim 1$
.
Equations (3.2) and (3.3) can be reduced to a single equation for
$u$
as a function of
$\alpha$
\begin{equation} \frac {\mathrm{d}u}{\mathrm{d}\alpha }=\frac {Fu^{2}}{4\mu \cos ^{2}\alpha }\left ( \gamma +\frac {F}{\cos \alpha }-u\right ) . \end{equation}
In turn, (2.18) and (2.20) can be reduced to a single boundary condition
\begin{equation} u=1\quad \text{at}\quad \alpha =\alpha _{0}. \end{equation}
To describe a free-falling curtain, one should let
\begin{equation} u\rightarrow \infty \quad \text{as}\quad \alpha \rightarrow -\frac {\pi }{2}, \end{equation}
in which case (3.2) yields a more specific asymptotic
\begin{equation} u=\frac {F}{\frac{\pi }{2}+\alpha }+\mathcal{O}(1)\quad \text{as}\quad \alpha \rightarrow -\frac {\pi }{2}. \end{equation}
Using (3.3), one can show that, in terms of the original variable
$l$
, the above result amounts to
\begin{align} u=\sqrt {2l}+\mathcal{O}(1)\quad \text{as}\quad l\rightarrow \infty , \end{align}
hence, condition (2.21) of finite viscous stress at infinity is satisfied.
Equation (3.4) and conditions (3.5)–(3.6) form a boundary-value problem for
$u(\alpha )$
. Since this is a first-order problem with two boundary conditions, it admits a solution only for certain values of the constant
$F$
, which should thus be treated as an eigenvalue.
Problem (3.4)–(3.6) was solved numerically by shooting: one solution was shot from
$\alpha =\alpha _{0}$
using boundary condition (3.5), and the other was shot from
$\alpha =-{1}/{2}\pi +\Delta \alpha$
(where
$\Delta \alpha$
is a small positive number) using asymptotic (3.7). The value of
$F$
was determined as the one for which the two solutions match at an intermediate point, which effectively means that they coincide.
As shown in the next two subsections, the existence of solution of problem (3.4)–(3.6) depends on whether the ejection angle
$\alpha _{0}$
is positive or negative. In the latter case, it is implied that
$\alpha _{0} \neq \pi /2$
, as vertically ejected steady curtains exist for all parameter values (and have also been observed experimentally by Finnicum et al. Reference Finnicum, Weinstein and Ruschak1993 and Roche et al. Reference Roche, Grand, Brunet, Lebon and Limat2006).
3.1. Oblique curtains ejected downwards (
$-\pi /2\lt \alpha _{0}\leqslant 0$
)
The region on the
$ ( \gamma ,\mu )$
plane where a solution of problem (3.4)–(3.6) exists for
$\alpha _{0}\lt 0$
is shown in figure 2(a). Note that, when the solution approaches the boundary of the existence region, the horizontal momentum flux vanishes,
$F\rightarrow 0$
– which provides a convenient tool for locating this boundary. One can also show that all non-positive values of
$\alpha _{0}$
correspond to the same curve
$F=0$
– hence, have the same existence region on the
$ ( \gamma ,\mu )$
plane (see Appendix C).
Figure 2. Existence of steady oblique curtains ejected downwards (
$-\pi /2\lt \alpha _{0}\leqslant 0$
) and examples thereof. (a) The existence region in the
$( \gamma ,\mu )$
plane (
$\gamma$
characterises surface tension and
$\mu$
viscosity, both are given by (2.2)). (b) Examples of trajectories (note that the
$x$
- and
$z$
-axes have different scales). (c) Streamwise velocity
$u$
vs
$x$
. Points (1)–(3) in panel (a) correspond to curves (1)–(3), respectively, in panels (b–c).
Once the solution
$u(\alpha )$
of boundary-value problem (3.4)–(3.6) is computed, one can readily calculate the curtain’s other characteristics – e.g. its shape. Some examples of such shapes are shown in figure 2(b); one can see that an increase of either
$\gamma$
or
$\mu$
gives rise to an increase of the solution’s curvature near the outlet, while the rest of the curtain becomes near-vertical (this trend can also be deduced through the analysis of the limit
$F\rightarrow 0$
in Appendix C).
The fact that a stronger surface tension (larger
$\gamma$
) bends the curtain stronger is not surprising, whereas the mechanism underlying viscous bending (larger
$\mu$
) is less intuitive. It is likely due to the strong imbalance between the streamwise and cross-stream viscous stresses in a thin curtain.
Figure 2(c) shows the graphs of the streamwise velocity for various values of
$\gamma$
,
$\mu$
and
$\alpha _{0}$
. One can see that, paradoxically, higher-viscosity curtains accelerate faster than the lower-viscosity ones. To resolve the paradox, note that viscosity does slow the curtain, but it also bends it downwards (see figure 2(b)) – which increases the input of gravitational energy, thereby overcompensating the viscous dissipation.
Finally, no solutions have been found if
$\gamma$
exceeds unity, which agrees with the conclusions of Weinstein et al. (Reference Weinstein, Ross, Ruschak and Barlow2019) and Della Pia et al. (Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023) that steady oblique curtains exist only if
${\textit{We}}\gt 1$
. This restriction, however, is weaker than the restriction found in the present work, as the curve in figure 2(a) is located to the left from the vertical straight line
$\gamma =1$
.
As seen later, this conclusion applies to curtains ejected upwards as well.
3.2. Oblique curtains ejected upwards (
$0\lt \alpha _{0}\lt \pi /2$
)
It has turned out that an increase of the ejection angle beyond
$\alpha _{0}=0$
expands the existence region of steady curtains on the
$ ( \gamma ,\mu )$
plane, and also gives rise to a subregion where two solutions exist (as illustrated in figure 3). The two solutions correspond to different values of
$F$
, and the one with the smaller
$F$
is more curved (as illustrated in the right-hand panel of figure 3). Note that the limit
$F\rightarrow 0$
still plays a role: this is how one of the two solutions disappears – hence, the curve corresponding to
$F=0$
is the lower boundary of the two-solution subregion in the
$ ( \gamma ,\mu )$
plane. At its upper boundary, in turn, the two solutions ‘merge’ and disappear, thus forming a fold bifurcation.
Figure 3. The same as in figure 2, but for curtains ejected upwards, at
$\alpha _{0}=75^{\circ }$
. Point (1) in panel (a) corresponds to curve (1) in panels (b–c), but point 2 corresponds to both curves (2u) and (2s) (the former/latter solutions are unstable/stable). The dotted curve in panel (a) is the same as the solid curve in panel (a) of figure 2.
As typically occurs with folds, one of the two merging solutions is unstable, and it is intuitively clear that the more curved curtain should be the unstable one. These two conclusions have been confirmed by deriving the linearised stability problem and exploring its solutions numerically in the
$ ( \gamma ,\mu )$
plane (no further details are given, due to the standard and straightforward nature of the procedures involved).
As before, no solutions have been found for
$\gamma \gt 1$
.
4. Evolving curtains
The fact that steady curtains exist only in a certain region of the parameter space makes one wonder how they evolve if their parameters happen to be outside this region. In §§ 4.1–4.2, this question is clarified for
$\gamma \lt 1$
and
$\gamma \gt 1$
, respectively.
To clarify the evolution, the governing equations were integrated numerically using the method of lines – i.e. the spatial coordinate
$l$
was discretised and the derivatives of the unknowns were approximated using upwind finite differences. The resulting set of ordinary differential equations in
$t$
was solved using an algorithm for integration of stiff problems. This approach was realised in MATLAB using the function ode23tb (in application to the problem at hand, it has turned out to be the fastest of MATLAB’s solvers for stiff ordinary differential equations).
4.1. Curtains with weak surface tension (
$\gamma \lt 1$
)
Assume for simplicity that, far from the outlet, the curtain is horizontal
\begin{align} \alpha \rightarrow 0\quad \text{as}\quad l\rightarrow \infty , \end{align}
spatially homogeneous
\begin{align} \frac {\partial u}{\partial l}\rightarrow 0,\quad \frac {\partial W}{\partial l}\rightarrow 0\quad \text{as}\quad l\rightarrow \infty , \end{align}
and free falling
\begin{align} \frac {\partial ^{2}z}{\partial t^{2}}\rightarrow 1\quad \text{as}\quad l\rightarrow \infty . \end{align}
Such a behaviour is obviously consistent with the governing equation and the zero-viscous-stress condition (2.21).
Figure 4. Evolution of curtains for initial conditions (4.4)–(4.5) with
$\alpha _{0}=-45^{\circ }$
. The values of
$(\gamma ,\mu )$
are indicated above the corresponding panels. (a) A steady state exists (shown by the doted curve); curves (0)–(6) correspond to
$t=0,1,2,3,4,6,8$
, respectively (note that these values are not equispaced). (b) No steady states exist; curves (0)–(6) correspond to
$t=0,1,2,3,5,7,10$
.
Equations (2.7) and (2.14)–(2.15) were integrated numerically with the full set of boundary conditions (described in Appendix B). Various initial conditions have been simulated, but the evolution turned out to depend mostly on whether a steady solution exists for the given values of
$ ( \gamma ,\mu )$
.
Figures 4(a) and 5(a) illustrate situations where a steady solution exists, and figures 4(b) and 5(b) illustrate situations where it does not. In all cases, the initial shape is
\begin{equation} \alpha =\alpha _{0}\operatorname {e}^{-l^{2}/2}\quad \text{at}\quad t=0, \end{equation}
which corresponds to a curtain that is curved near the outlet but otherwise flat, and is ejected at angle
$\alpha _{0}$
(figures 4 and 5 illustrate
$\alpha _{0}\leqslant 0$
and
$\alpha _{0}\gt 0$
, respectively). For the remaining variables, the simplest initial conditions were used,
\begin{equation} W=2,\quad u=1,\quad \frac {\partial \alpha }{\partial t}=0\quad \text{at}\quad t=0 \end{equation}
i.e. the curtain’s initial width, velocity and angular velocity are spatially uniform and equal to their outlet values.
Figure 5. Evolution of curtains for initial conditions (4.4)–(4.5). The values of
$(\gamma ,\mu )$
are indicated above the corresponding panels. (a) Two steady states exist, the stable one is shown by the doted curve; curves (0)–(5) correspond to
$t=0,2,4,6,12,30$
, respectively (note that these values are not equispaced). (b) No steady states exist; curves (0)–(5) correspond to
$t=0,5,10,20,30,40$
.
Figures 4(a) and 5(a) show that, if a steady state exists, the solution tends to it – thus, illustrating its stability. Figures 4(b) and 5(b), in turn, show what happens if no steady state exists: the curtain rotates around the outlet, with its uppermost part wrapping around the outlet tighter and tighter. Note that figures 4(b) and 5(a) show the solutions for the same values of
$\gamma$
and
$\mu$
, but different values of
$\alpha _{0}$
. In the former case, the point
$ ( \gamma ,\mu )$
lies outside the existence region of steady solutions, whereas in the latter the same point lies within.
It can be assumed that the rotation and wrapping would eventually give rise to the teapot effect – i.e. the curtain would touch the lower outside of the outlet and flow down it. The specifics of such evolution depend on the shape of the outlet’s outside.
4.2. Curtains with strong surface tension (
$\gamma \gt 1$
)
If
$\gamma$
exceeds unity, it is impossible to make sure that condition (2.26) holds for all
$l$
and all
$t$
. This implies that the hyperbolic part of the governing set has a characteristic whose velocity may change sign. As a result, setting up the upwind approximation of the derivatives is no longer straightforward: one should decompose the hyperbolic part of the governing system into three equations for the corresponding Riemann invariants, and adjust the finite-difference approximation of one of these equations depending on whether condition (2.26) holds at the current point and time.
There is a simpler option, however: instead of the general case, one can examine the asymptotic limit of near-critical curtains, where
$\gamma$
differs from unity by a small margin. It turns out that this limit is described by an equation of the first order – hence, no Riemann decomposition is needed in this case. Quantitatively, such a model describes only a small fraction of the allowable parameter range, but it should still elucidate the qualitative behaviour of all curtains with
$\gamma \gt 1$
.
4.2.1. Near-critical curtains
Introduce
$\Delta \gamma =\gamma -1$
and assume that
\begin{equation} \left \vert \Delta \gamma \right \vert \ll 1. \end{equation}
Physically, this condition means that the contributions of capillarity and advection to the speed of upstream capillary waves nearly cancel each other.
It can be verified by trial and error that the most inclusive asymptotic regime corresponds to the following scaling:
\begin{align}&\qquad\qquad\qquad\qquad\quad\,\, \mu =\mathcal{O}(\Delta \gamma ), \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad t=\mathcal{O}(1), \end{align}
\begin{align}& l=\mathcal{O}(\Delta \gamma ),\quad x=\mathcal{O}(\Delta \gamma ),\quad z=\mathcal{O}(\Delta \gamma ),\quad \alpha =\mathcal{O}(1), \end{align}
\begin{align}&\quad\, \Delta u=u-1=\mathcal{O}(\Delta \gamma ),\quad \Delta W=W-2=\mathcal{O} (\Delta \gamma ). \end{align}
Physically, condition (4.7) implies that viscosity is on par with the combined effect of capillarity/advection; (4.8) indicates that the time scale remains the same as before, whereas (4.9) shows that the spatial scale becomes small (observe that
$\alpha$
changes by an order of unity over a small region near the outlet). Conditions (4.10), in turn, indicate that the curtain’s velocity and thickness do not change much in this region, with the implication that they do so over a larger distance (to be discussed later).
Keeping in set (2.7)–(2.10) the leading-order terms only, one can verify that the equations for
$x$
and
$z$
decouple from the rest of the set, the equation for
$W$
and boundary conditions (2.17) yield
$\Delta W=-2\Delta u$
, the two remaining equations become
\begin{align}& 2\frac {\partial \alpha }{\partial t}+\left ( \Delta u-\Delta \gamma -4\mu \frac {\partial u}{\partial l}\right ) \frac {\partial \alpha }{\partial l}=-\cos \alpha . \end{align}
\begin{align}&\qquad\qquad \frac {\partial \Delta u}{\partial l}=4\mu \frac {\partial ^{2}\Delta u}{\partial l^{2}}-\sin \alpha , \end{align}
and boundary conditions (2.20)–(2.21) yield
\begin{align}&\quad\,\, \Delta u=0\quad \text{at}\quad l=0, \end{align}
\begin{align}& \left \vert \frac {\partial \Delta u}{\partial l}\right \vert \lt \infty \quad \text{as}\quad l\rightarrow \infty . \end{align}
To formulate boundary conditions for
$\alpha$
, observe that (4.11) is a first-order hyperbolic equation and the velocity of its characteristic is
\begin{align} c_{3}=\Delta u-\Delta \gamma -4\mu \frac {\partial \Delta u}{\partial l}, \end{align}
which is the near-critical limit of characteristic (2.24) of the general set. If
$c_{3}$
is negative at the outlet, no boundary condition is needed for
$\alpha (t,l)$
at
$l=0$
, otherwise
$\alpha (t,0)$
is determined by the orientation of the nozzle. These two options can be written as a single equality: using (4.13) to evaluate
$c_{3}$
at
$l=0$
, one obtains
\begin{equation} \left ( \alpha -\alpha _{0}\right ) \operatorname {H}\left ( -\Delta \gamma -4\mu \frac {\partial \Delta u}{\partial l}\right ) =0\quad \text{at}\quad l=0, \end{equation}
where
$\operatorname {H}$
is Heaviside’s step function.
One can also assume for simplicity that
$c_{3}$
is positive at infinity,
\begin{equation} \Delta u-\Delta \gamma -4\mu \frac {\partial \Delta u}{\partial l}\gt 0\quad \text{as}\quad l\rightarrow \infty , \end{equation}
with the implications that no boundary condition is needed for
$\alpha$
as
$l\rightarrow \infty$
. Requirement (4.17) is not very restrictive; in particular, it holds for the experiment of Della Pia et al. (Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023).
Subject to requirement (4.17) and given a suitable initial condition for
$\alpha$
, (4.11)–(4.12) and boundary conditions (4.13)–(4.16) fully determine
$\alpha$
and
$u$
.
4.2.2. Solutions of the near-critical problem
4.2.2.1. Steady curtains
It can be readily verified that (4.11)–(4.12) and boundary conditions (4.13), (4.16) are satisfied for
\begin{align}&\qquad\qquad\qquad\qquad\qquad \alpha =\arctan \left ( \tan \alpha _{0}-\frac {l}{F}\right )\! , \end{align}
\begin{align}& \Delta u=-\frac {1}{4\mu }\operatorname {e}^{l/(4\mu )}\int _{0}^{l}\left [ \sqrt {F^{2}+\left ( F^{2}\tan \alpha _{0}-l\right ) ^{2}}\operatorname {sgn} F+\Delta \gamma \right ] \operatorname {e}^{-l^{\prime }/(4\mu )}\mathrm{d} l^{\prime }, \end{align}
where
$F$
is an undetermined constant. It remains to satisfy boundary condition (4.14); substituting expression (4.19) into it, one obtains
\begin{equation} \left ( \operatorname {sgn}F\right ) \int _{0}^{\infty }\sqrt {F^{2}+\left ( F\tan \alpha _{0}-l^{\prime }\right ) ^{2}}\operatorname {e}^{-l^{\prime }/(4\mu )}\mathrm{d}l^{\prime }=-4\mu \Delta \gamma . \end{equation}
Given specific values of
$\alpha _{0}$
,
$\mu$
, and
$\Delta \gamma$
, this equation determines
$F$
.
The following comments are in order.
-
(i) If
$\Delta \gamma \lt 0$
, solutions (4.18)–(4.20) are the near-critical limit of the steady solutions of the general set (2.7)–(2.10) examined in § 3. -
(ii) If
$\Delta \gamma \geqslant 0$
, the general equations do not have steady solutions, but the near-critical equations do. Note also that, if
$\Delta \gamma$
is positive, (4.20) yields a negative
$F$
– in which case solution (4.18) describes upward-bending curtains. These are the high-viscosity near-critical analogues of the solutions found by Benilov (Reference Benilov2019) and shown to be unstable by Benilov (Reference Benilov2023). Their physical meaning – or rather the lack thereof – will be discussed in § 5.1.
4.2.2.2. Flat rotating curtains
Equation (4.11) and boundary condition (4.16) admit a time-dependent solution describing a flat curtain rotating around the outlet
\begin{equation} \alpha =\arctan \frac {\operatorname {e}^{C-t}-1}{2\operatorname {e}^{\left ( C-t\right ) /2}}, \end{equation}
where
$C$
is an arbitrary constant. Solving then the boundary-value problem (4.12)–(4.14), one obtains the following expression for the streamwise velocity:
\begin{equation} \Delta u=-l\sin \alpha (t). \end{equation}
It can be verify that, if the curtain’s initial slope is not too large,
\begin{align} \sin \alpha \lt \frac {\Delta \gamma }{4\mu }\quad \text{at}\quad t=0, \end{align}
solution (4.21) satisfies boundary condition (4.16).
It follows from (4.21) that
\begin{align}&\quad \frac {\partial \alpha }{\partial t}\gt 0\quad \text{for all }t, \end{align}
\begin{align}& \alpha \rightarrow -\frac {\pi }{2}\quad \text{as}\quad t\rightarrow \infty , \end{align}
i.e. the curtain is rotating clockwise and eventually becomes vertical.
It has turned out that, if
$\Delta \gamma \geqslant 0$
, all solutions of problem (4.11)–(4.16) with initial conditions satisfying assumption (4.17) exhibit the same behaviour.
4.2.2.3. Numerical solutions
Typical evolution of curtains governed by problem (4.11)–(4.16) subject to assumption (4.17) is illustrated in figure 6 for
\begin{equation} \Delta \gamma =1,\quad \mu =0.5, \end{equation}
and for the following initial conditions:
-
(i)
(4.27)which is a perturbed version of the flat-curtain solution (4.21);
\begin{equation} \alpha =-\frac {\pi }{12}\big ( 1+5l\operatorname {e}^{-l^{2}/2}\big ) \quad \text{at}\quad t=0, \end{equation}
-
(ii)
(4.28)which would be a steady solution if
\begin{equation} \alpha =\arctan \left ( -2+\sqrt {3}-\frac {l}{0.6692}\right ) \quad \text{at}\quad t=0, \end{equation}
$\Delta \gamma$
in (4.26) were changed from
$\Delta \gamma =1$
to
$\Delta \gamma =-1$
. This initial condition is intended to loosely mimic the experiment of Della Pia et al. (Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023), in which the evolution started from a steady curtain with
$\gamma \lt 1$
, then the ejection velocity was reduced resulting in a flow with
$\gamma \gt 1$
.
Evidently, in both cases, the curtain rotates clockwise towards the vertical, but does not pass it. The latter feature distinguishes them from curtains with
$\gamma \lt 1$
– compare figure 6 with figures 4(b) and 5(b).
4.2.3. The curtain’s structure near the outlet and at infinity
-
(i) Note that, in both solutions illustrated in figure 6, the curtain’s slope at
$l=0$
varies with time and, thus, does not necessarily match the orientation of the nozzle. More generally, the fact that (4.11) may not require a boundary condition at
$l=0$
implies that a boundary layer may exists near the outlet, where the flow changes its direction from that prescribed by the nozzle to that arising naturally in the outside region. This boundary layer, in fact, is visible in the experiments and computations of Della Pia et al. (Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023) (see their figures 10 and 11(d), respectively).One can further conjecture that, if
$\gamma$
exceeds unity by an order-one margin (which is the case in the simulations of Della Pia et al. Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023), the length of this boundary layer is comparable to its width – hence, the slim-curtain approximation used in this paper is inapplicable. Near-critical curtains, however, should be tractable similarly to the treatment of non-viscous near-critical curtains by Benilov (Reference Benilov2021). -
(ii) Substituting solution (4.22) into assumption (4.17), one can show that the latter holds only if the curtain is initially ejected at a negative angle. One can furthermore argue that any curtain whose large-distance part is oriented downwards
(4.29)satisfies assumption (4.17). In such a curtain, the velocity grows with the distance from the outlet
\begin{equation} \alpha (0,l)\lt 0\quad \text{as}\quad l\rightarrow \infty , \end{equation}
(4.30)hence, assumption (4.17) holds.
\begin{align} u\rightarrow \infty \quad \text{as}\quad l\rightarrow \infty , \end{align}
As for upward-ejected curtains – for which condition (4.29) does not hold – the near-critical (4.11) requires a boundary condition at infinity. To derive such, one should recall that (4.11) applies only at small distances,
$l\sim \Delta \gamma$
(as follows from scaling (4.9)), then consider the large-distance problem, and eventually match the two solutions.
Qualitatively, one can conjecture that upward -jected curtains either rotate clockwise and eventually flow downwards (like their downwards-ejected counterparts) – or rotate counter-clockwise towards the upward direction and eventually become unstable.
5. Discussion
5.1. Upward-bending curtains
As mentioned in the Introduction, some of the previous models predict curtains which bend upwards and keep flowing upwards until the initial supply of kinetic energy runs out and a stagnation point is formed (Keller & Weitz Reference Keller and Weitz1957; Benilov Reference Benilov2019). It was later shown, however, that such solutions are unstable with respect to perturbation growing downstream and eventually blowing up near the stagnation point (Benilov Reference Benilov2023).
There are also models derived exclusively for the near-outlet region (e.g. Benilov Reference Benilov2021, and the near-critical approximation of the present paper), which do not include the stagnation point. The upward-bending solutions appear stable in this case, but one should still be able to show their instability once the near-outlet (inner) solution is matched to the large-distance (outer) one. It is also telling that stagnation points destabilise the flow in many other settings (e.g. Friedlander & Vishik Reference Friedlander and Vishik1991; Benilov & Lapin Reference Benilov and Lapin2014; Dubovskaya & Benilov Reference Dubovskaya and Benilov2021).
Interestingly, the general model (2.7)–(2.10) – which is applicable to the whole curtain – does not admit upward-bending steady curtains at all, so the question of their stability is irrelevant. It is also worth noting that the non-existence proof (given in Appendix D) is based on the divergence of a certain integral at the stagnation point.
One should conclude that it is the stagnation point that eliminates upward-bending curtains, even though it manifests itself differently in different models: either as instability of the steady solution or its non-existence. A alternative argument in favour of non-existence of upward-bending curtains is given by Weinstein et al. (Reference Weinstein, Ross, Ruschak and Barlow2019), based on the fact that the inviscid curtain equations are singular at points where the local Weber number equals unity.
One way or another, upward-bending curtains cannot be observed in an experiment.
5.2. How can the results of this paper be tested experimentally?
As mentioned before, the constraint imposed on the existence of steady stable curtains by high viscosity and surface tension is stronger than the previously known constraint due to high surface tension alone,
$\gamma \lt 1$
(Weinstein et al. Reference Weinstein, Ross, Ruschak and Barlow2019; Benilov Reference Benilov2023; Della Pia et al. Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023). It is still instructive to illustrate both constraints graphically for a specific liquid and in dimensional form – say, in terms of the ejection velocity
$U$
and half-width
$H$
of the outlet.
Table 1. The density
$\rho$
, kinematic viscosity
$\nu$
and surface tension
$\sigma$
of water (Lindstrom & Mallard Reference Lindstrom and Mallard1997; Wagner & Pruß Reference Wagner and Pruß2002), ethylene glycol (Bohne, Fischer & Obermeier Reference Bohne, Fischer and Obermeier1984; MEGlobal 2024) and glycerol (Adamenko et al. Reference Adamenko, Bulavin, Ilyin, Zelinsky and Moroz2006; Cheng Reference Cheng2008), all at
$20\,^{\circ }\mathrm{C}$
.
Figure 7. Existence of steady oblique curtains ejected downwards (
$-\pi /2\lt \alpha _{0}\leqslant 0)$
, on the
$(H,U)$
plane (
$H$
is the half-width of the outlet,
$U$
the ejection velocity) for: (a) water, (b) ethylene glycol, (c) glycerol. This figure is the dimensional equivalent of figure 2(a). The dashed line corresponds to the curve
$\gamma =1$
, where
$\gamma$
is given by (2.2). The black dots correspond to the rows of table 2.
This has been done for water, ethylene glycol and glycerol (in all cases for oblique downwards-ejected curtains,
$-\pi /2\lt \alpha _{0}\leqslant 0$
) – see figure 7. The physical parameters of the three liquids are listed in table 1.
Interestingly, the two curves are almost indistinguishable for water, but are separated by a noticeable gap for the other two liquids. The difference is due to the higher viscosity of ethylene glycol and glycerol, which exceed that of water by factors of approximately 20 and 1100, respectively. To further illustrate the implications of the difference in viscosity, table 2 lists, and figure 7 shows, examples of steady oblique curtains that exist in less viscous liquid(s) but not in more viscous one(s).
Table 2. Examples of existence (‘yes’) and non-existence (‘no’) of steady oblique curtains ejected downwards (
$-\pi /2\lt \alpha _{0}\leqslant 0$
), for various liquids. A ‘yes’ implies that
$\gamma$
does not exceed a certain critical value
$\varGamma (\mu )$
, which is a stronger requirement than the previously known condition
$\gamma \lt 1$
(
$\gamma$
characterises surface tension and
$\mu$
viscosity, both are given by (2.2)). Here,
$H$
is the half-width of the outlet,
$U$
the ejection velocity. The rows of this table correspond to the black dots in figure 7, with both ordered from top to bottom.
Thus, to observe the effect of viscosity on the dynamics of curtains, one should experiment with a liquid whose viscosity is at least comparable to that of ethylene glycol. One should also use a sufficiently wide nozzle (observe that, as
$H\rightarrow 0$
, the solid and dashed curves in all three panels of figure 7 become indistinguishable).
6. Summary and concluding remarks
The following results have been reported in this paper.
(i) An oblique curtain can be steady only if it satisfies a constraint of the form
\begin{equation} \gamma \lt \varGamma (\mu ), \end{equation}
where
$\gamma$
and
$\mu$
are given by expressions (2.2) and characterise the effects of surface tension and viscosity, respectively. For all curtains ejected downwards, the function
$\varGamma (\mu )$
(plotted in figure 2(a)) does not depend on the ejection angle
$\alpha _{0}$
– but once
$\alpha _{0}$
becomes positive,
$\varGamma (\mu )$
starts to increase (see an example in figure 3(a)). In all cases, however, constraint (6.1) is more restrictive than the previously known condition for the existence of steady stable curtains,
$\gamma \lt 1$
(or, equivalently,
${\textit{We}}\gt 1$
).
(ii) Oblique curtains not satisfying constraint (6.1) evolve: they typically rotate around the outlet, but the finer details of the evolution depend on whether
$\gamma$
exceeds unity (in addition to violating condition (6.1)): curtains with
$\gamma \gt 1$
rotate towards the downward direction (as in figure 6), whereas those with
$\varGamma (\mu )\lt \gamma \lt 1$
rotate beyond it (as in figures 4(b) and 5(b)). Since in both cases the curtain does not follow the natural arch-like trajectory, both can be viewed as instances of the ‘teapot effect’ (e.g. Kistler & Scriven Reference Kistler and Scriven1994) – but only the latter can result in the classical teapot pattern, where the jet or curtain runs down the lower outside of the container holding the liquid.
Both of the above scenarios involve strong bending of the curtain near the outlet. One might speculate that the high-curvature region can destabilise the whole curtain, rendering its evolution volatile (the slim-curtain approximation used in this paper is inapplicable in this case, hence its prediction reflects only the global tendency). The most likely scenario involves dislodgement of one of the curtain’s contact lines from the edge of the outlet and formation of a ‘bulge’ where the flow abruptly changes direction – as suggested by the experiments and simulations of Kistler & Scriven (Reference Kistler and Scriven1994) and Della Pia et al. (Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023). After the bulge, the curtain is either vertical (see figures 10 and 11(d) of Della Pia et al. Reference Della Pia, Antoniades, Ioannidis, Wejko, Barlow, Chiatto, Weinstein and de Luca2023) or directed at an angle beyond the vertical (see figures 18, 20 and 22, 23 of Kistler & Scriven Reference Kistler and Scriven1994). The existence of a steady oblique curtain in the latter case can be explained by bulge-induced changes in the curtain parameters:
$H$
increases,
$U$
decreases – and although the flux
$UH$
remains the same,
$\gamma$
becomes larger and
$\mu$
smaller. As a result, a steady oblique curtain may come into existence.
Declaration of interests
The author reports no conflicts of interests.
Appendix A. Derivation of asymptotic (2.7)–(2.10)
A.1. The Navier–Stokes equations and boundary conditions
Consider the setting described in the beginning of § 2 and introduce the fluid velocity
$\boldsymbol{u}$
, its streamwise and cross-stream components
$u$
and
$w$
, respectively, and the pressure
$p$
. Introduce also orthogonal curvilinear coordinates
$ ( l,\tau )$
related to the Cartesian coordinates
$ ( x,z )$
by
\begin{align} x=x(l,\tau ,t),\quad z=z(l,\tau ,t), \end{align}
where
$t$
is the time. Physically,
$ ( l,\tau )$
are the streamwise and cross-stream coordinates, respectively.
In addition to the non-dimensional variables (2.5)–(2.6), introduce
\begin{align} \tau _{\textit{nd}}=\frac {\tau }{H},\quad w_{\textit{nd}}=\frac {w}{U},\quad p_{\textit{nd}}=\frac {p}{\rho U^{2}}. \end{align}
In terms of the non-dimensional variables, the Navier–Stokes equations are (subscripts
$_{\textit{nd}}$
omitted)
\begin{align}&\,\, \left ( \frac {D\boldsymbol{u}}{Dt}\right ) _{l}+\dfrac {1}{h_{l}}\dfrac {\partial p}{\partial l}=\frac {\mu }{\varepsilon ^{2}}\left ( \boldsymbol{\nabla} ^{2}\boldsymbol{u}\right ) _{l}-\dfrac {1}{h_{l}}\dfrac {\partial z}{\partial l}, \end{align}
\begin{align}& \left ( \frac {D\boldsymbol{u}}{Dt}\right ) _{\tau }+\dfrac {1}{\varepsilon h_{\tau } }\dfrac {\partial p}{\partial \tau }=\frac {\mu }{\varepsilon }\left ( \boldsymbol{\nabla} ^{2}\boldsymbol{u}\right ) _{\tau }-\frac {1}{\varepsilon h_{\tau }}\dfrac {\partial z}{\partial \tau }, \end{align}
\begin{align}&\qquad\qquad \frac {\partial \left ( uh_{\tau }\right ) }{\partial l}+\frac {1}{\varepsilon }\frac {\partial \left ( wh_{l}\right ) }{\partial \tau }=0, \end{align}
where
$\varepsilon$
and
$\mu$
are given by (2.2),
\begin{equation} h_{l}=\sqrt {\left ( \frac {\partial x}{\partial l}\right ) ^{2}+\left ( \frac {\partial z}{\partial l}\right ) ^{2}},\quad h_{\tau }=\frac {1}{\varepsilon }\sqrt {\left ( \frac {\partial x}{\partial \tau }\right ) ^{2}+\left ( \frac {\partial z}{\partial \tau }\right ) ^{2}}, \end{equation}
are the Lamé coefficients and
\begin{align} \left ( \frac {D\boldsymbol{u}}{Dt}\right ) _{l}&=\frac {\partial u}{\partial t}+\frac {1}{h_{l}}\left [ u-\frac {1}{h_{l}}\left ( \frac {\partial x}{\partial l}\frac {\partial x}{\partial t}+\frac {\partial z}{\partial l}\frac {\partial z}{\partial t}\right ) \right ] \left ( \frac {\partial u}{\partial l}+\frac {w}{\varepsilon h_{\tau }}\dfrac {\partial h_{l}}{\partial \tau }\right ) +\frac {1}{\varepsilon h_{\tau }} \nonumber \\&\quad\times\left [ w-\frac {1}{\varepsilon h_{\tau }}\left ( \frac {\partial x}{\partial \tau }\frac {\partial x}{\partial t}+\frac {\partial z}{\partial \tau }\frac {\partial z}{\partial t}\right ) \right ] \frac {\partial u}{\partial \tau }+\frac {w}{h_{l}h_{\tau }}\left [ \frac {1}{\varepsilon }\left ( \dfrac {\partial x}{\partial l}\frac {\partial ^{2}x}{\partial t\,\partial \tau }+\dfrac {\partial z}{\partial l}\frac {\partial ^{2}z}{\partial t\,\partial \tau }\right ) \right . \nonumber \\&\quad -\left . \frac {1}{\varepsilon ^{3}h_{\tau }^{2}}\left ( \frac {\partial x}{\partial \tau }\frac {\partial x}{\partial t}+\frac {\partial z}{\partial \tau }\frac {\partial z}{\partial t}\right ) \left ( \dfrac {\partial x}{\partial l}\frac {\partial ^{2}x}{\partial \tau ^{2}}+\dfrac {\partial z}{\partial l} \frac {\partial ^{2}z}{\partial \tau ^{2}}\right ) -w\frac {\partial h_{\tau } }{\partial l}\right ]\! ,\end{align}
\begin{align} \left ( \frac {D\boldsymbol{u}}{Dt}\right ) _{\tau }&=\frac {\partial u_{\tau } }{\partial t}+\frac {1}{h_{l}}\left [ u_{l}-\frac {1}{h_{l}}\left ( \frac {\partial x}{\partial l}\frac {\partial x}{\partial t}+\frac {\partial z}{\partial l}\frac {\partial z}{\partial t}\right ) \right ] \frac {\partial u_{\tau }}{\partial l}\nonumber \\&\quad +\frac {1}{h_{\tau }}\left [ u_{\tau }-\frac {1}{\varepsilon h_{\tau }}\left ( \frac {\partial x}{\partial \tau }\frac {\partial x}{\partial t}+\frac {\partial z}{\partial \tau }\frac {\partial z}{\partial t}\right ) \right ] \left ( \frac {1}{\varepsilon }\frac {\partial u_{\tau }}{\partial \tau }+\frac {u_{l}} {h_{l}}\frac {\partial h_{\tau }}{\partial l}\right ) \nonumber \\&\quad +\frac {u_{l}}{h_{l}h_{\tau }}\left [ \frac {1}{\varepsilon }\left ( \dfrac {\partial x}{\partial \tau }\frac {\partial ^{2}x}{\partial t\,\partial l}+\dfrac {\partial z}{\partial \tau }\frac {\partial ^{2}z}{\partial t\,\partial l}\right ) \right . \nonumber \\&\quad -\left . \frac {1}{\varepsilon h_{l}^{2}}\left ( \frac {\partial x}{\partial l}\frac {\partial x}{\partial t}+\frac {\partial z}{\partial l}\frac {\partial z}{\partial t}\right ) \left ( \dfrac {\partial x}{\partial \tau }\frac {\partial ^{2}x}{\partial l^{2}}+\dfrac {\partial z}{\partial \tau }\frac {\partial ^{2}z}{\partial l^{2}}\right ) -\frac {u_{l}}{\varepsilon } \frac {\partial h_{l}}{\partial \tau }\right ]\! ,\end{align}
\begin{align} \big( \boldsymbol{\nabla} ^{2}\boldsymbol{u}\big) _{l}&=\frac {1}{h_{l}}\frac {\partial }{\partial l}\left \{ \frac {1}{h_{l}h_{\tau }}\left [ \varepsilon ^{2} \frac {\partial \left ( uh_{\tau }\right ) }{\partial l}+\varepsilon \frac {\partial \left ( wh_{l}\right ) }{\partial \tau }\right ] \right \} \nonumber \\ &\quad +\frac {1}{h_{\tau }}\frac {\partial }{\partial \tau }\left \{ \frac {1}{h_{l} h_{\tau }}\left [ \frac {\partial \left ( uh_{l}\right ) }{\partial \tau }-\varepsilon \frac {\partial \left ( wh_{\tau }\right ) }{\partial l}\right ] \right \} ,\end{align}
\begin{align} \big(\boldsymbol{\nabla} ^{2}\boldsymbol{u}\big) _{\tau }&=\frac {1}{h_{\tau }}\frac {\partial }{\partial \tau }\left \{ \frac {1}{h_{\tau }h_{l}}\left [ \frac {\partial \left ( uh_{\tau }\right ) }{\partial l}+\frac {1}{\varepsilon }\frac {\partial \left ( wh_{l}\right ) }{\partial \tau }\right ] \right \} \nonumber \\ &\quad +\frac {1}{h_{l}}\frac {\partial }{\partial l}\left \{ \frac {1}{h_{\tau }h_{l} }\left [ \varepsilon \frac {\partial \left ( wh_{\tau }\right ) }{\partial l}-\frac {\partial \left ( uh_{l}\right ) }{\partial \tau }\right ] \right \} .\end{align}
Let the boundaries of the curtain be described by
\begin{align} \tau =\pm \frac {W}{2}, \end{align}
where
$W(t,l)$
is the curtain’s local thickness. Then the non-dimensional kinematic and dynamic conditions are
\begin{align} &\frac {\partial W}{\partial t}+\dfrac {1}{h_{l}}\left [ u-\dfrac {1}{h_{l} }\left ( \frac {\partial x}{\partial l}\frac {\partial x}{\partial t} +\frac {\partial z}{\partial l}\frac {\partial z}{\partial t}\right ) \right ] \dfrac {\partial W}{\partial l}=\pm \dfrac {2}{\varepsilon h_{\tau }}\left [ w-\dfrac {1}{\varepsilon h_{\tau } }\left ( \frac {\partial x}{\partial \tau }\frac {\partial x}{\partial t} +\frac {\partial z}{\partial \tau }\frac {\partial z}{\partial t}\right ) \right ] \nonumber \\ &\quad \text{if}\quad \tau =\pm \frac {W}{2}, \end{align}
\begin{align} &\left [ 2\varepsilon ^{2}\left ( \frac {1}{h_{l}}\frac {\partial u}{\partial l}+\frac {w}{\varepsilon h_{l}h_{\tau }}\frac {\partial h_{l}}{\partial \tau }\right ) -\frac {\varepsilon ^{2}p-\varepsilon ^{3}\gamma c_{\pm }}{\mu }\right ] \dfrac {1}{h_{l}}\dfrac {\partial W}{\partial l}=\pm \left [ \frac {h_{l}}{h_{\tau }}\frac {\partial }{\partial \tau }\left ( \frac {u}{h_{l}}\right )\right .\nonumber \\ &\quad +\left .\frac {\varepsilon h_{\tau }}{h_{l}}\frac {\partial }{\partial l}\left ( \frac {w}{h_{\tau }}\right ) \right ] \dfrac {2}{h_{\tau } }\quad \text{if}\quad \tau =\pm \frac {W}{2}, \end{align}
\begin{align} &\left [ \frac {h_{l}}{h_{\tau }}\frac {\partial }{\partial \tau }\left ( \frac {u}{h_{l}}\right ) +\frac {\varepsilon h_{\tau }}{h_{l}}\frac {\partial }{\partial l}\left ( \frac {w}{h_{\tau }}\right ) \right ] \dfrac {1}{2h_{l}}\dfrac {\partial W}{\partial l}=\pm \left [ 2\left ( \frac {1}{\varepsilon h_{\tau }}\frac {\partial w} {\partial \tau }+\frac {u}{h_{l}h_{\tau }}\frac {\partial h_{\tau }}{\partial l}\right )\right .\nonumber \\ & -\left.\frac {p-\varepsilon \gamma c_{\pm }}{\mu }\right ] \dfrac {2}{h_{\tau }}\quad \text{if}\quad \tau =\pm \frac {W}{2}, \end{align}
where
$\gamma$
is given by (2.2) and the curvatures of the curtain’s boundaries are
\begin{equation} c_{\pm }=-\frac {1}{h_{l}h_{\tau }}\left \{ \varepsilon \frac {\partial }{\partial l}\left [ \frac {\frac{h_{\tau }}{h_{l}}\frac{\partial W}{\partial l}} {2\sqrt {\left (\frac{\varepsilon }{h_{l}}\right ) ^{2}+\left (\frac{1}{h_{\tau }}\right ) ^{2}}}\right ] \pm \frac {1}{\varepsilon }\frac {\partial }{\partial \tau }\left [ \frac {\frac{h_{l}}{h_{\tau }}}{\sqrt {\left (\frac{\varepsilon }{h_{l}}\right ) ^{2}+\left (\frac{1}{h_{\tau }}\right ) ^{2}}}\right ] \right \} . \end{equation}
Let the coordinates
$ ( l,\tau )$
have unit Jacobian. Non-dimensionally, this condition and the orthogonality condition are
\begin{equation} \frac {\partial x}{\partial l}\frac {\partial z}{\partial \tau }-\frac {\partial x}{\partial \tau }\frac {\partial z}{\partial l}=\varepsilon ,\quad \frac {\partial x}{\partial l}\frac {\partial x}{\partial \tau }+\frac {\partial z}{\partial l}\frac {\partial z}{\partial \tau }=0. \end{equation}
Given boundary conditions at the outlet and initial conditions, (A3)–(A16) fully describe the flow characteristics and the Cartesian coordinates as functions of
$l$
,
$\tau$
and
$t$
.
A.2. Asymptotic analysis
As shown by Benilov (Reference Benilov2019), the solution of (A16) can be sought as power series in
$\varepsilon \tau$
, which yields
\begin{align} x&=x^{(0)}-\varepsilon \tau \sin \alpha -\frac {\varepsilon ^{2}\tau ^{2}}{2} \frac {\partial \alpha }{\partial l}\sin \alpha +\mathcal{O}\bigl(\varepsilon ^{3} \tau ^{3}\bigr), \end{align}
\begin{align} z&=z^{(0)}+\varepsilon \tau \cos \alpha +\frac {\varepsilon ^{2}\tau ^{2}}{2} \frac {\partial \alpha }{\partial l}\cos \alpha +\mathcal{O}\bigl(\varepsilon ^{3} \tau ^{3}\bigr), \end{align}
where the centreline’s coordinates
$ ( x^{(0)}(t,l),z^{(0)}(t,l) )$
and slope
$\alpha (t,l)$
are inter-related by
\begin{equation} \frac {\partial x^{(0)}}{\partial l}=\cos \alpha ,\quad \frac {\partial z^{(0)} }{\partial l}=\sin \alpha . \end{equation}
It is convenient to replace
$ ( u,w )$
with
$ ( \hat {u} ,\hat {w} )$
such that
\begin{equation} u=\frac {1}{h_{l}}\left ( \frac {\partial x}{\partial l}\frac {\partial x}{\partial t}+\frac {\partial z}{\partial l}\frac {\partial z}{\partial t}\right ) +\hat {u},\quad w=\frac {1}{\varepsilon h_{\tau }}\left ( \frac {\partial x}{\partial \tau }\frac {\partial x}{\partial t}+\frac {\partial z}{\partial \tau }\frac {\partial z}{\partial t}\right ) +\varepsilon \hat {w}. \end{equation}
Physically,
$ ( \hat {u},\hat {w} )$
are the velocities of the flow relative to the motion of the curtain as a whole, and
$\varepsilon$
that appears next to
$\hat {w}$
reflects the weakness of the cross-stream flow by comparison with the streamwise one.
Expand the flow characteristics in powers of
$\varepsilon \tau$
\begin{align} \hat {u}(t,l,\tau )&=\hat {u}^{(0)}(t,l)+\varepsilon \tau \,\hat {u}^{(1)} (t,l)+\varepsilon ^{2}\tau ^{2}\hat {u}^{(2)}(t,l)+\mathcal{O}\bigl(\varepsilon ^{3}\tau ^{3}\bigr),\nonumber \\ \hat {w}(t,l,\tau )&=\tau \hat {w}^{(1)}(t,l)+\varepsilon \tau ^{2}\hat {w} ^{(2)}(t,l)+\mathcal{O}\bigl(\varepsilon ^{2}\tau ^{3}\bigr),\nonumber \\ p(t,l,\tau )&=p^{(0)}(t,l)+\varepsilon \tau p^{(1)}(t,l)+\varepsilon ^{2}\tau ^{2}p^{(2)}(t,l)+\mathcal{O}\bigl(\varepsilon ^{3}\tau ^{3}\bigr). \end{align}
Substituting these expansions and (A17)–(A18) into expressions (A20), the continuity (A5) and boundary conditions (A13)–(A15), one obtains
\begin{equation}\begin{aligned} \hat {u}^{(1)}&=-\frac {\partial \alpha }{\partial l}\hat {u}^{(0)}, &\hat {u}^{(2)}&=2\frac {\partial \hat {u}^{(0)}}{\partial l}\frac {1}{\tau _{\pm }} \dfrac {\partial \tau _{\pm }}{\partial l}+\frac {1}{2}\left [ \frac {\partial ^{2}\hat {u}^{(0)}}{\partial l^{2}}-\left ( \frac {\partial \alpha }{\partial l}\right ) ^{2}\hat {u}^{(0)}\right ]\! ,\nonumber \\ \hat {w}^{(1)}&=-\frac {\partial \hat {u}_{l}^{(0)}}{\partial l}, &\hat {w}^{(2)}&=-\frac {\partial \alpha }{\partial l}\frac {\partial \hat {u}^{(0)} }{\partial l}.\nonumber \\ p^{(0)}&=-2\mu \frac {\partial \hat {u}^{(0)}}{\partial l}, &p^{(1)} &=2\mu \left ( \frac {\partial ^{2}\alpha }{\partial l\partial t}+\frac {\partial ^{2}\alpha }{\partial l^{2}}\hat {u}^{(0)}-\frac {\partial \alpha }{\partial l}\frac {\partial \hat {u}^{(0)}}{\partial l}\right ) -\frac {\partial \alpha }{\partial l}\frac {2\gamma }{W}, \end{aligned}\end{equation}
while the rest of the governing set can be reduced to
\begin{align} &W\left ( \frac {\partial \hat {u}^{(0)}}{\partial t}+\frac {\partial ^{2}x^{(0)} }{\partial t^{2}}\cos \alpha +\frac {\partial ^{2}z^{(0)}}{\partial t^{2}} \sin \alpha +\hat {u}^{(0)}\frac {\partial \hat {u}^{(0)}}{\partial l}\right ) =4\mu \frac {\partial }{\partial l}\left ( W\frac {\partial \hat {u}^{(0)}}{\partial l}\right ) -W\sin \alpha ,\nonumber \\&\,\, -\frac {\partial ^{2}x^{(0)}}{\partial t^{2}}\sin \alpha +\frac {\partial ^{2}z^{(0)}}{\partial t^{2}}\cos \alpha +2\hat {u}^{(0)}\frac {\partial \alpha }{\partial t}+\frac {\partial \alpha }{\partial l}\left ( \hat {u}^{(0)2} -\frac {2\gamma }{W}-4\mu \frac {\partial \hat {u}^{(0)}}{\partial l}\right ) =-\cos \alpha ,\nonumber \\&\qquad\qquad\qquad\qquad\qquad\qquad\quad \frac {\partial W}{\partial t}+\dfrac {\partial \left ( \hat {u}^{(0)}W\right ) }{\partial l}=0. \end{align}
Finally, omitting the superscript
$^{(0)}$
and hats from the above equations and (A19), one obtains (2.7)–(2.10), as required.
Appendix B. Boundary conditions for (2.7), (2.14)–(2.15)
Introduce
$\beta (t,l)$
and
$D(t,l)$
such that
\begin{align}&\qquad\qquad\qquad \frac {\partial \alpha }{\partial t}=\beta , \end{align}
\begin{align}& \frac {\partial u}{\partial t}+u\frac {\partial u}{\partial l}-\frac {4\mu } {W}\dfrac {\partial }{\partial l}\left ( W\frac {\partial u}{\partial l}\right ) =C, \end{align}
and rewrite (2.14)–(2.15) in the form
\begin{align} \frac {\partial \beta }{\partial t}&=-C\frac {\partial \alpha }{\partial l} -\frac {\partial }{\partial l}\left [ 2u\beta +\left ( u^{2}-\frac {2\gamma } {W}-4\mu \frac {\partial u}{\partial l}\right ) \frac {\partial \alpha }{\partial l}\right ]\! , \end{align}
\begin{align} \frac {\partial C}{\partial l}&=\beta ^{2}+\left [ 2u\beta +\left ( u^{2} -\frac {2\gamma }{W}-4\mu \frac {\partial u}{\partial l}\right ) \frac {\partial \alpha }{\partial l}\right ] \frac {\partial \alpha }{\partial l}. \end{align}
Equations (2.7) and (B1)–(B3) determine the evolution of
$W$
,
$\alpha$
,
$u$
and
$\beta$
, while (B4) determines
$C$
. Counting the derivatives with respect to
$l$
in these equations, one understands that
$W$
,
$\alpha$
,
$\beta$
and
$C$
need one boundary condition each, and
$u$
needs two. All these conditions except those for
$\beta$
and
$C$
are delivered by (2.17), (2.18), (2.20) and (2.21). The boundary condition for
$\beta$
can be obtained by differentiating (2.18) with respect to
$t$
and recalling (B1)
\begin{align} \beta =0\quad \text{at}\quad l=0. \end{align}
The condition for
$C$
follows from (B2) and the boundary conditions for the other variables
\begin{align} C=-\sin \alpha _{0}\quad \text{at}\quad l=0. \end{align}
Appendix C. Limit
$F\rightarrow 0$
of boundary-value problem (3.4)–(3.6)
Replace
$\alpha$
with
\begin{equation} \phi =\frac {1}{F}\left ( \alpha +\frac {\pi }{2}\right ), \end{equation}
and rewrite (3.4), boundary condition (3.5) and asymptotics (3.7) (which can be used instead of boundary condition (3.6)) in the form
\begin{align}& \frac {\mathrm{d}u}{\mathrm{d}\phi }=-\frac {F^{2}u^{2}}{4\mu \sin ^{2}F\phi }\left ( u-\gamma -\frac {F}{\sin F\phi }\right )\! , \end{align}
\begin{align}&\qquad u=1\quad \text{at}\quad \phi =\frac {1}{F}\left ( \alpha _{0}+\frac {\pi } {2}\right )\! , \end{align}
\begin{align}&\qquad\,\, u=\frac {1}{\phi }+\mathcal{O}(1)\quad \text{as}\quad \phi \rightarrow 0. \end{align}
Assuming that
$F\ll 1$
and keeping in (C2)–(C3) the leading-order terms only, one obtains
\begin{align}& \frac {\mathrm{d}u}{\mathrm{d}\phi }=-\frac {u^{2}}{4\mu \phi ^{2}}\left ( u-\gamma -\frac {1}{\phi }\right )\! , \end{align}
\begin{align}&\qquad u=1\quad \text{as}\quad \phi \rightarrow \infty , \end{align}
whereas (C4) remains the same.
Since (C5) is of the first order, its solution cannot generally satisfy both boundary conditions (C4) and (C6). Thus, either
$\mu$
or
$\gamma$
should be treated as an eigenvalue. As a result, boundary-value problem (C4)–(C6) determines a curve on the
$ ( \gamma ,\mu )$
plane; it was computed and is shown in figure 2(a).
Observe that the ejection angle
$\alpha _{0}$
does not appear in in problem (C4)–(C6) – hence, the above-mentioned curve where
$F=0$
is the same for all values of
$\alpha _{0}$
.
Appendix D. Non-existence of upward-bending curtains
If upward-bending curtains exist, they are described by (3.4) and boundary condition (3.5) at the outlet – plus an additional condition is required at the terminal (stagnation) point where the fluid runs out of its initial supply of kinetic energy. Denoting the curtain’s slope at this point by
$\alpha _{t}$
, one should require
\begin{equation} u=0\quad \text{at}\quad \alpha =\alpha _{t}. \end{equation}
Divide (3.4) by
$u^{2}$
and integrate it from
$\alpha =\alpha _{0}$
to
$\alpha =\alpha _{t}$
. Evaluating the integral on the left-hand side using boundary conditions (3.5) and (D1), one obtains
\begin{equation} -\infty =\int _{\alpha _{0}}^{\alpha _{t}}\frac {F}{4\mu \cos ^{2}\alpha }\left ( \gamma +\frac {F}{\cos \alpha }-u\right ) \mathrm{d}\alpha . \end{equation}
This equality can hold only if its right-hand is infinite – hence,
$\alpha _{t}=\pi /2$
. Note, however, that
\begin{align} \frac {F}{4\mu \cos ^{2}\alpha }\left ( \gamma +\frac {F}{\cos \alpha }-u\right ) \sim \frac {F^{2}}{4\mu \cos ^{3}\alpha }\rightarrow +\infty \quad \text{as} \quad \alpha \rightarrow \frac {\pi }{2}-0, \end{align}
hence, the integral on the right-hand side of (D2) is
$+\infty$
, and (D2) does not hold anyway.
Open access
