1 Introduction
Dynamic wetting or dewetting occurs when one fluid adjacent to a solid substrate is displaced by another one, and is encountered in a variety of natural and industrial problems. It is well known that classical hydrodynamic description of dynamic wetting leads to divergence of viscous stress at the moving contact line (Huh & Scriven Reference Huh and Scriven1971), where the liquid–liquid interface intersects the wall. A range of contact-line models have been proposed to relieve the contact-line singularity, such as wall slip and interface diffusion (Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). The slip model is popular in hydrodynamic theory, which predicts a dynamic contact angle induced by viscous bending of the interface due to the viscous flow close to the contact line (Voinov Reference Voinov1976; Cox Reference Cox1986). On the other hand, an alternative description of moving contact lines is molecular-kinetic theory, which incorporates a local contact-line friction associated with thermally activation of liquid molecules over the solid (Blake & Haynes Reference Blake and Haynes1969; Blake Reference Blake2006; Seveno et al. Reference Seveno, Vaillant, Rioboo, Adão, Conti and De Coninck2009; Blake & De Coninck Reference Blake and De Coninck2011). The idea of molecular-kinetic theory has also been employed to interpret the contact-line pinning/depinning on a substrate with defects (e.g. Rolley & Guthmann Reference Rolley and Guthmann2007).
An interesting feature of a moving contact line for partially wetting fluids is the existence of a maximum speed, beyond which a dynamical wetting transition occurs, leading to entrainment of films, drops or gas bubbles (Snoeijer & Andreotti Reference Snoeijer and Andreotti2013). Film entrainment is the hydrodynamic basis for coating (Weinstein & Ruschak Reference Weinstein and Ruschak2004), which may exhibit complex patterns (Thiele Reference Thiele2014; Wilczek et al. Reference Wilczek, Zhu, Chi, Thiele and Gurevich2017). In contrast, wetting failure associated with gas entrainment is usually detrimental to industrial applications. We focus on wetting transitions associated with liquid entrainment, and recent progress on gas entrainment can be found in Marchand et al. (Reference Marchand, Chan, Snoeijer and Andreotti2012), Vandre, Carvalho & Kumar (Reference Vandre, Carvalho and Kumar2012, Reference Vandre, Carvalho and Kumar2013, Reference Vandre, Carvalho and Kumar2014), Chan et al. (Reference Chan, Srivastava, Marchand, Andreotti, Biferale, Toschi and Snoeijer2013) and Sprittles (Reference Sprittles2015).
Dynamical wetting transition is typically investigated in the situation of dip coating, in which a plate is withdrawn from a liquid bath at a speed
$U$
. The governing parameter is the capillary number
$Ca=\unicode[STIX]{x1D707}U/\unicode[STIX]{x1D70E}$
, where
$\unicode[STIX]{x1D707}$
is the liquid viscosity and
$\unicode[STIX]{x1D70E}$
is the surface tension. For small
$Ca$
, the problem has been studied under the lubrication approximation. Early investigations of film deposition concerned the situation of complete wetting, with particular interest in the speed dependence of the film thickness (Landau & Levich Reference Landau and Levich1942; Wilson Reference Wilson1982). When the liquid is partially wetting, the critical speed for film entrainment is roughly proportional to
$\unicode[STIX]{x1D703}_{e}^{3}$
, where
$\unicode[STIX]{x1D703}_{e}$
is the equilibrium contact angle, as first proposed by de Gennes (Reference de Gennes1986). This scaling was theoretically demonstrated by Eggers (Reference Eggers2004, Reference Eggers2005) for inclined plates and later by Chan, Snoeijer & Eggers (Reference Chan, Snoeijer and Eggers2012) and Qin & Gao (Reference Qin and Gao2018) for vertically aligned plates. At the onset of wetting transition, the meniscus behaves as if the liquid were completely wetting, as the apparent contact angle vanishes. In fact, the wetting transition was predicted to be triggered by a saddle-node bifurcation (Snoeijer et al.
Reference Snoeijer, Andreotti, Delon and Fermigier2007; Chan et al.
Reference Chan, Snoeijer and Eggers2012; Tseluiko, Galvagno & Thiele Reference Tseluiko, Galvagno and Thiele2014; Wilczek et al.
Reference Wilczek, Zhu, Chi, Thiele and Gurevich2017), unless the inclination angle of the plate is sufficiently small (Ziegler, Snoeijer & Eggers Reference Ziegler, Snoeijer and Eggers2009; Galvagno et al.
Reference Galvagno, Tseluiko, Lopez and Thiele2014). The critical condition predicted by lubrication theory was confirmed by recent numerical simulations of Gao et al. (Reference Gao, Li, Feng, Ding and Lu2016), although the experimentally observed value was overestimated (Snoeijer et al.
Reference Snoeijer, Delon, Fermigier and Andreotti2006; Delon et al.
Reference Delon, Fermigier, Snoeijer and Andreotti2008). Beyond the transition, the film evolution is well described by two-dimensional hydrodynamics. The film exhibits complicated morphologies depending on the withdrawal velocity (Snoeijer et al.
Reference Snoeijer, Delon, Fermigier and Andreotti2006; Gao et al.
Reference Gao, Li, Feng, Ding and Lu2016). Specifically, a thick film is obtained for
$Ca$
weakly above the critical value, with the thickness selected by the contact line; a further increase of
$Ca$
gives rise to the traditional Landau–Levich film (Landau & Levich Reference Landau and Levich1942), but connected to a ridge close to the contact line through a capillary shock. Dynamical wetting transition and film deposition by withdrawing a solid cylinder were investigated by Maleki et al. (Reference Maleki, Reyssat, Quéré and Golestanian2007) and Chan, Gueudre & Snoeijer (Reference Chan, Gueudre and Snoeijer2011).
The wetting transition in confined geometries, on the other hand, is rarely considered. The contact-line dynamics in a capillary tube was usually studied in the parametric region below the transition in the case of advancing contact lines (Hoffman Reference Hoffman1975, Reference Hoffman1983; Fermigier & Jenffer Reference Fermigier and Jenffer1991). Quéré (Reference Quéré1991) investigated the displacement of a drop in a capillary tube and found that a liquid film was entrained beyond a critical velocity. In contrast to planar geometries, Callegari, Calvo & Hulin (Reference Callegari, Calvo and Hulin2005) reported a dependence on the film thickness of the dewetting velocity of the entrained liquid film. Recently, experiments on forced dewetting were performed by Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018), who found that the entrainment of liquid films would generate a gas finger and eventually cause the formation of Taylor bubbles (Taylor Reference Taylor1961). This process is due to the axisymmetric confinement and differs from the formation of gas fingers induced by Saffman–Taylor instability in Hele-Shaw cells (Saffman & Taylor Reference Saffman and Taylor1958); in the latter the flow is not confined in the direction parallel to the wall.
In this work, we carry out numerical and theoretical investigations of the forced wetting in a capillary tube. The critical speed for wetting transition is analytically predicted. The simulations not only reproduce the experimental observations of Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018), but also provide detailed flow behaviours, based on which an insightful analysis of the bubble length is performed. The theoretical prediction is in good agreement with the numerical results.
Figure 1. Schematic of the axisymmetric two-phase flow in a capillary, where the liquid is displaced by the gas at a constant flow rate
$Q$
.
2 Governing equations and numerical method
Consider the two-phase flow in a capillary tube of radius
$R$
, as shown in figure 1. The liquid is displaced by the gas with a fixed volume flow rate
$Q$
from right to left, and both phases are incompressible and Newtonian. The liquid is partially wetting and the gas–liquid interface intersects the wall at a microscopic contact angle, measured within the liquid. For high-viscosity liquids, the contact angle is assumed to remain at the equilibrium value
$\unicode[STIX]{x1D703}_{e}$
, independent of the contact-line speed (Eggers Reference Eggers2005). We further assume that the tube radius is much smaller than the capillary length and the velocity is sufficiently slow, such that both gravitational and inertial effects can be neglected. We expect the interface to behave as a stable meniscus when the flow rate is less than a threshold, beyond which a liquid film is entrained on the tube wall, as typically observed in dip-coating geometries (Snoeijer & Andreotti Reference Snoeijer and Andreotti2013).
A diffuse-interface method was employed to deal with the deformation of the gas–liquid interface as well as the moving contact line (Jacqmin Reference Jacqmin2004; Qian, Wang & Sheng Reference Qian, Wang and Sheng2006; Yue, Zhou & Feng Reference Yue, Zhou and Feng2010). A phase-field variable
$\unicode[STIX]{x1D719}$
was introduced to distinguish the two fluids, where
$\unicode[STIX]{x1D719}=1$
in the liquid and
$\unicode[STIX]{x1D719}=-1$
in the gas. The fluids are mixed and
$\unicode[STIX]{x1D719}$
varies continuously across the interface, which is characterized by a small thickness
$\unicode[STIX]{x1D716}$
. The evolution of pressure
$p$
, velocity
$\boldsymbol{v}$
and
$\unicode[STIX]{x1D719}$
are governed by the Stokes–Cahn–Hilliard equations,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle -\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v}+G\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \displaystyle \unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}t+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=\displaystyle \unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FB}^{2}G. & \displaystyle\end{eqnarray}$$
Here, the viscosity of the mixed fluids is given by
$\unicode[STIX]{x1D707}=(1/2)(1+\unicode[STIX]{x1D719})\unicode[STIX]{x1D707}_{l}+(1/2)(1-\unicode[STIX]{x1D719})\unicode[STIX]{x1D707}_{g}$
, which recovers the liquid viscosity
$\unicode[STIX]{x1D707}_{l}$
and the gas viscosity
$\unicode[STIX]{x1D707}_{g}$
away from the interface;
$\unicode[STIX]{x1D6FE}$
is the mobility. The term
$G\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}$
in the momentum equation represents the effect of interfacial forces, where
$G=\unicode[STIX]{x1D6EC}[-\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D719}+(\unicode[STIX]{x1D719}^{2}-1)\unicode[STIX]{x1D719}/\unicode[STIX]{x1D716}^{2}]$
is the bulk chemical potential, and
$\unicode[STIX]{x1D6EC}$
is the mixing energy density, relating to the interfacial tension
$\unicode[STIX]{x1D70E}$
through
$\unicode[STIX]{x1D70E}=2\sqrt{2}\unicode[STIX]{x1D6EC}/(3\unicode[STIX]{x1D716})$
.
We investigate the problem in a reference frame moving to the left with the mean velocity
$U=Q/(\unicode[STIX]{x03C0}R^{2})$
; accordingly, the tube wall moves to the right at the velocity
$U$
. This allows one to track the interface evolution in a fixed computational domain. Since the flow is axisymmetric, we employ a rectangular domain of size
$R\times 40R$
in the
$(r,z)$
plane, which is long enough to avoid any end effect on the interface evolution. Initially, the interface is simply flat and located in the middle of the tube. Parabolic velocity profiles with zero flux across the tube are imposed at the two ends of the domain. The remaining boundary conditions include no-slip condition on the wall and symmetric condition at the axis. In addition, boundary conditions for
$\unicode[STIX]{x1D719}$
are
$\mathbf{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}G=0$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6EC}\mathbf{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}=\left\{\begin{array}{@{}ll@{}}-f_{w}^{\prime }(\unicode[STIX]{x1D719}),\quad & \text{on the wall,}\\ 0,\quad & \text{otherwise},\end{array}\right.\end{eqnarray}$$
where
$f_{w}(\unicode[STIX]{x1D719})=(1/4)\unicode[STIX]{x1D70E}\unicode[STIX]{x1D719}(\unicode[STIX]{x1D719}^{2}-3)\cos \unicode[STIX]{x1D703}_{e}$
is the wall free energy. This form of wall energy can guarantee that the profile of
$\unicode[STIX]{x1D719}$
around the contact line is consistent with that across the interface in the bulk, and the contours of
$\unicode[STIX]{x1D719}$
intersect the wall at the angle
$\unicode[STIX]{x1D703}_{e}$
(Jacqmin Reference Jacqmin2000; Yue et al.
Reference Yue, Zhou and Feng2010). Moreover, the dynamic contact angle in reality may deviate from the equilibrium one; this effect can be realized by further introducing a wall relaxation in (2.4) (Jacqmin Reference Jacqmin2000), which was not adopted in the present work for simplicity.
In the following, we scale all lengths by
$R$
, velocity by
$\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D707}_{l}$
and time by
$\unicode[STIX]{x1D707}_{l}R/\unicode[STIX]{x1D70E}$
. Apart from the contact angle
$\unicode[STIX]{x1D703}_{e}$
, the problem is governed by the following dimensionless numbers: the capillary number
$Ca=\unicode[STIX]{x1D707}_{l}U/\unicode[STIX]{x1D70E}$
, the viscosity ratio
$\unicode[STIX]{x1D707}_{l}/\unicode[STIX]{x1D707}_{g}$
, the Cahn number
$Cn=\unicode[STIX]{x1D716}/R$
and
$S=\sqrt{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D707}_{l}}/R$
. We fixed the viscosity ratio
$\unicode[STIX]{x1D707}_{l}/\unicode[STIX]{x1D707}_{g}=50$
, which is sufficiently large to render the gas phase negligible, as commonly done for receding contact lines. Note that the contact-line movement is realized through interface diffusion, and
$S$
is a dimensionless diffusion length, a microscopic length scale usually encountered in contact-line models. An alternative definition of the diffusion length is
$S^{\ast }=\sqrt{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D707}^{\ast }}/R$
with
$\unicode[STIX]{x1D707}^{\ast }=\sqrt{\unicode[STIX]{x1D707}_{l}\unicode[STIX]{x1D707}_{g}}$
(Yue et al.
Reference Yue, Zhou and Feng2010), relating to the present definition through
$S=(\unicode[STIX]{x1D707}_{l}/\unicode[STIX]{x1D707}_{g})^{1/4}S^{\ast }$
. The value of
$S$
plays an important role on the contact-line velocity. In reality the microscopic length scale of the contact-line singularity is typically a few nanometres, corresponding to
$S\sim O(10^{-5})$
or even smaller for a tube radius of a millimetre (Zhao et al.
Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018). Such a small length is difficult to resolve in diffuse-interface simulations. We have instead employed
$S=0.005$
and 0.01 in our simulations. As shown by (3.9) below, the critical capillary number scales as
$1/\text{ln}\,S^{-1}$
. The large value of
$S$
used here would lead to an overestimation of the critical capillary number; an overestimation of the dewetting velocity of liquid films can be expected as well. Nevertheless, the qualitative picture of interface evolution remains the same, with the calculated process of film deposition occurring at larger values of
$Ca$
in the simulations. Moreover, we chose
$Cn=0.005$
or 0.01 such that the interface is sufficiently thin and the criteria
$Cn<4S^{\ast }$
for the sharp-interface limit is satisfied (Yue et al.
Reference Yue, Zhou and Feng2010). Unless otherwise specified, the results presented below were obtained at
$Cn=0.01$
and
$S=0.01$
.
The governing equations are solved using a finite-element method coupled with a second-order time stepping scheme (Yue et al.
Reference Yue, Zhou, Feng, Ollivier-Gooch and Hu2006; Zhou et al.
Reference Zhou, Yue, Feng, Ollivier-Gooch and Hu2010). The mesh is adaptively refined around the interface to capture the steep variation of
$\unicode[STIX]{x1D719}$
. The method has been tested for a wide range of interfacial flows, including the film deposition on a plate (Gao et al.
Reference Gao, Li, Feng, Ding and Lu2016).
3 Steady meniscus
When the flow rate is lower than a threshold, or
$Ca<Ca_{cr}$
, the fluid interface maintains a steady profile after a transient relaxation, and migrates with the average flow velocity,
$U$
, with respect to the tube wall. Figure 2(a) shows three representative interfacial profiles for
$\unicode[STIX]{x1D703}_{e}=30^{\circ }$
and
$S=0.01$
. As
$Ca$
increases, the interface becomes more curved. The axial distance between the contact line and the meniscus tip,
$\unicode[STIX]{x0394}z_{cl}$
, is shown in the figure 2(b) as a function of
$Ca$
to further characterize the interfacial deformation for
$S=0.005$
and 0.01. It is seen
$\unicode[STIX]{x0394}z_{cl}$
increases monotonically from a value determined by
$\unicode[STIX]{x1D703}_{e}$
, and the steady solution ceases to exist above a critical capillary number,
$Ca_{cr}$
. As expected, a smaller value of
$S$
simply leads to a narrower range of
$Ca$
for the existence of the steady menisci.
Figure 2. (a) Representative interfacial profiles when
$Ca\leqslant Ca_{cr}=0.01$
for
$\unicode[STIX]{x1D703}_{e}=30^{\circ }$
and
$S=0.01$
. The interfaces are shifted such that the tips coincide. (b) Distance between the contact line and the meniscus tip,
$\unicode[STIX]{x0394}z_{cl}$
, as a function of
$Ca$
for
$\unicode[STIX]{x1D703}_{e}=30^{\circ }$
and
$S=0.005$
and 0.01. Inset shows the variation of
$\unicode[STIX]{x1D703}_{ap}$
, relating to
$\unicode[STIX]{x0394}z_{cl}$
through (3.1). The symbols represent the diffuse-interface simulations, and the curves correspond to the lubrication theory with the slip length
$\unicode[STIX]{x1D706}=S$
.
Since
$Ca\ll 1$
, the macroscopic shape of the meniscus is primarily determined by the balance between surface tension and capillary pressure, while the viscous flow only plays an important role in the vicinity of the contact line. Therefore, it is reasonable to approximate the interface away from the contact line by a spherical cap with an apparent contact angle
$\unicode[STIX]{x1D703}_{ap}$
, which is related to
$\unicode[STIX]{x0394}z_{cl}$
through
$$\begin{eqnarray}\unicode[STIX]{x0394}z_{cl}=\frac{1-\sin \unicode[STIX]{x1D703}_{ap}}{\cos \unicode[STIX]{x1D703}_{ap}},\end{eqnarray}$$
and shown in the inset of figure 2(b). The maximum
$\unicode[STIX]{x0394}z_{cl}$
of the steady solution is approximately 1, corresponding to
$\unicode[STIX]{x1D703}_{ap}=0$
, a criterion usually encountered for wetting transition with receding contact lines (Eggers Reference Eggers2004, Reference Eggers2005; Chan et al.
Reference Chan, Snoeijer and Eggers2012).
Following Eggers (Reference Eggers2005), an asymptotic solution of the interfacial deformation and the critical capillary number can be obtained for small contact angles based on lubrication theory. The resulting solution is expected to apply to finite contact angles as well, as demonstrated by Qin & Gao (Reference Qin and Gao2018). In theoretical analyses, it is convenient to employ a Navier slip model rather than Cahn–Hilliard diffusion to regularize the contact-line singularity. With a slip length
$\unicode[STIX]{x1D706}_{s}$
, the ratio of length scales
$\unicode[STIX]{x1D706}\equiv \unicode[STIX]{x1D706}_{s}/R\ll 1$
is a counterpart of
$S$
in our diffuse-interface simulations. Although there are no explicit rules connecting
$\unicode[STIX]{x1D706}$
and
$S$
, we will show that lubrication theory agrees with diffuse-interface simulations when putting
$\unicode[STIX]{x1D706}=S$
, at least for the parameters considered. The thickness of the liquid film is denoted by
$h(x)$
, where
$x\equiv z_{cl}-z$
with
$z_{cl}$
being the position of the contact line. The interfaces close to the contact line and away from it, denoted by
$h_{in}(x)$
and
$h_{out}(x)$
respectively, are treated separately and matched to obtain the full solution. The matching procedure was essentially adopted from Eggers (Reference Eggers2005), with a replacement of the curvature of the outer solution, as shown below.
The inner solution is determined by the balance between the surface tension and the viscous force in a region of size
$O(\unicode[STIX]{x1D706})$
. Since
$\unicode[STIX]{x1D706}\ll 1$
, the curvature of the tube wall can be neglected, and the contact-line solution on a flat plate can be adopted. This solution was first proposed by Duffy & Wilson (Reference Duffy and Wilson1997), and then extended by Eggers (Reference Eggers2005) to allow for slip. For the present purpose, of relevance is the asymptotic form for
$\unicode[STIX]{x1D706}\ll x\ll 1$
,
$$\begin{eqnarray}h_{in}(x)=\unicode[STIX]{x1D6FF}^{1/3}\left[\frac{\unicode[STIX]{x1D703}_{e}^{2}\unicode[STIX]{x1D705}_{y}}{6\unicode[STIX]{x1D706}}x^{2}+\unicode[STIX]{x1D703}_{e}b_{y}x+O(\unicode[STIX]{x1D706})\right],\end{eqnarray}$$
where
$$\begin{eqnarray}\unicode[STIX]{x1D705}_{y}(s_{1};\unicode[STIX]{x1D6FF})=\frac{2^{1/3}\unicode[STIX]{x1D6FD}^{2}}{\unicode[STIX]{x03C0}^{2}\text{Ai}^{2}(s_{1})},\qquad b_{y}(s_{1})=-\frac{2^{2/3}\text{Ai}^{\prime }(s_{1})}{\text{Ai}(s_{1})},\end{eqnarray}$$
are associated with the curvature and slope of the interface close to the contact line, respectively;
$\text{Ai}(s_{1})$
is the Airy function with
$s_{1}$
to be determined and
$\unicode[STIX]{x1D6FF}\equiv 3Ca/\unicode[STIX]{x1D703}_{e}^{3}$
is a reduced capillary number. The value of
$\unicode[STIX]{x1D6FD}^{2}$
associated with
$\unicode[STIX]{x1D705}_{y}$
has the form
$$\begin{eqnarray}\unicode[STIX]{x1D6FD}^{2}=\unicode[STIX]{x03C0}\exp [-1/(3\unicode[STIX]{x1D6FF})]/2^{2/3},\end{eqnarray}$$
as first proposed by Eggers (Reference Eggers2005) using a matching procedure at the region of
$x\gtrsim \unicode[STIX]{x1D706}$
. (Note that
$\unicode[STIX]{x1D709}$
is to be replaced by e
$\unicode[STIX]{x1D709}$
in (31) and (33) of Eggers (Reference Eggers2005).) On the other hand, the outer solution can be well approximated as a spherical cap, which is described as
$(1-h_{out}(x))^{2}+(x+\tan \unicode[STIX]{x1D703}_{ap})^{2}=1/\cos ^{2}\unicode[STIX]{x1D703}_{ap}$
. For
$x\ll 1$
, it can be written as
$$\begin{eqnarray}h_{out}(x)=\unicode[STIX]{x1D703}_{ap}x+\textstyle \frac{1}{2}x^{2}+O(x^{3}),\end{eqnarray}$$
where we have used the condition
$\unicode[STIX]{x1D703}_{ap}\ll 1$
since
$\unicode[STIX]{x1D703}_{ap}<\unicode[STIX]{x1D703}_{e}$
for receding contact lines.
Comparing (3.2) to (3.5), we obtain two matching conditions:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}_{ap}=\unicode[STIX]{x1D6FF}^{1/3}\unicode[STIX]{x1D703}_{e}b_{y}(s_{1}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 1=\frac{\unicode[STIX]{x1D6FF}^{1/3}\unicode[STIX]{x1D703}_{e}^{2}\unicode[STIX]{x1D705}_{y}(s_{1};\unicode[STIX]{x1D6FF})}{3\unicode[STIX]{x1D706}}. & \displaystyle\end{eqnarray}$$
The latter can be rewritten as
$$\begin{eqnarray}\text{Ai}^{2}(s_{1})=\frac{\unicode[STIX]{x1D6FF}^{1/3}\unicode[STIX]{x1D703}_{e}^{2}\exp [-1/(3\unicode[STIX]{x1D6FF})]}{2^{1/3}\times 3\unicode[STIX]{x03C0}\unicode[STIX]{x1D706}},\end{eqnarray}$$
which determines
$s_{1}$
in terms of
$\unicode[STIX]{x1D6FF}$
,
$\unicode[STIX]{x1D703}_{e}$
and
$\unicode[STIX]{x1D706}$
. Then substituting
$s_{1}$
into (3.6) gives
$\unicode[STIX]{x1D703}_{ap}$
and hence
$\unicode[STIX]{x0394}z_{cl}$
. The left-hand side of (3.8) has a global maximum at
$s_{1}=s_{max}=-1.018\,79\ldots .$
When
$Ca$
and hence
$\unicode[STIX]{x1D6FF}$
increase, the right-hand side of (3.6) may exceed the max
$\text{Ai}^{2}(s_{max})$
, and the steady-state solution will be lost. The critical capillary number is accordingly found as
$$\begin{eqnarray}Ca_{cr}=\frac{\unicode[STIX]{x1D703}_{e}^{3}}{9}\left[\ln \left(\frac{Ca_{cr}^{1/3}\unicode[STIX]{x1D703}_{e}}{18^{1/3}\unicode[STIX]{x03C0}\text{Ai}^{2}(s_{max})\unicode[STIX]{x1D706}}\right)\right]^{-1}.\end{eqnarray}$$
Note that the
$\unicode[STIX]{x1D703}_{e}^{3}$
dependence was also proposed by Quéré (Reference Quéré1991), without predicting the logarithmic correction related to
$\unicode[STIX]{x1D703}_{e}$
and
$Ca_{cr}$
. Furthermore, since
$\text{Ai}^{\prime }(s_{max})=0$
, the apparent contact angle
$\unicode[STIX]{x1D703}_{ap}$
vanishes at the threshold, according to (3.6). This is consistent with our numerical simulations, as discussed above.
Figure 3. Formation of a Taylor bubble for
$\unicode[STIX]{x1D703}_{e}=50^{\circ }$
and
$Ca=0.028$
, which is slightly larger than
$Ca_{cr}=0.0273$
. A moving reference frame is employed such that the tube wall moves to the right with a dimensionless velocity
$Ca$
.
A comparison between the lubrication theory and diffuse-interface simulations is presented in figure 2(b) for
$\unicode[STIX]{x1D703}_{e}=30^{\circ }$
. The separation of length scales is arbitrarily chosen as
$\unicode[STIX]{x1D706}=S$
. It is seen that the lubrication theory satisfactorily captures the variation of the interface deformation as well as the critical capillary number, even though the contact angle is not small.
Figure 4. Variation of the contact-line velocity relative to the tube wall,
$Ca_{cl}$
, as a function of
$Ca$
, for representative values of
$\unicode[STIX]{x1D703}_{e}$
.
4 Interfacial dynamics above
$Ca_{cr}$
4.1 Formation of Taylor bubbles
When
$Ca>Ca_{cr}$
, the meniscus becomes unstable and a liquid film is entrained on the tube wall. Figure 3 illustrates the interface evolution for
$\unicode[STIX]{x1D703}_{e}=50^{\circ }$
and
$Ca=0.028$
. Initially, the gas phase penetrates into the liquid, leading to a gas finger. Both the finger tip and the contact line travel at well-defined velocities, but in the opposite directions in the moving frame. While the gas finger gets elongated, the contact line collects the liquid in the film into a dewetting rim, which grows and then collapses into a liquid slug, leading to the formation of a long bubble, known as a Taylor bubble (Taylor Reference Taylor1961). The newly formed meniscus after the rim collapse remains unstable; therefore this process repeats, eventually giving rise to a train of travelling bubbles. The length of the bubble depends on the flow rate and the wettability of the wall. This process of bubble formation has also been reported by recent experiments of Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018). To quantify the bubble length, a systematic description of the contact-line and gas-finger dynamics is required. The front of the gas finger and the region close to the contact line can be studied separately as the liquid film between them remains uniform.
4.2 Moving contact line
Figure 4 shows the contact-line velocity relative to the tube wall,
$Ca_{cl}$
, as a function of
$Ca$
for representative values of
$\unicode[STIX]{x1D703}_{e}$
. The diagonal line
$Ca_{cl}=Ca$
marks a contact-line velocity that equals the mean flow velocity. For the solution of stable menisci,
$Ca_{cl}$
is simply identical to the mean velocity. At the onset of wetting transition, the contact-line velocity drops suddenly and then stabilizes, once the film is entrained, to a value
$Ca^{\ast }$
that depends on
$\unicode[STIX]{x1D703}_{e}$
but apparently not on the driving flow rate. It corresponds to the dewetting velocity of a film on a solid, which has been extensively studied for the case of a plane substrate (Redon, Brochard-Wyart & Rondelez Reference Redon, Brochard-Wyart and Rondelez1991; Snoeijer & Eggers Reference Snoeijer and Eggers2010). Note that Callegari et al. (Reference Callegari, Calvo and Hulin2005) reported a dependence of the dewetting velocity on the film thickness in a capillary tube, where a substantial variation of the film thickness occurs. However, for a fixed contact angle, the range of film thickness can be covered by the present simulations for is much narrower than the experiments of Callegari et al. (Reference Callegari, Calvo and Hulin2005), and the variation of the dewetting velocity is too weak to be resolved.
The present simulations show that wetting transition is supposed to occur at a critical capillary number,
$Ca_{cr}$
, larger than the dewetting velocity of a liquid film,
$Ca^{\ast }$
. In contrast, it was found in the experiments of Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018) that the wetting transition occurred at
$Ca=Ca^{\ast }$
rather than
$Ca_{cr}$
, such that the jump of contact-line velocity was not observed. Interestingly, a similar precritical transition was also encountered in a dip-coating geometry; in the experiments of Snoeijer et al. (Reference Snoeijer, Delon, Fermigier and Andreotti2006) and Delon et al. (Reference Delon, Fermigier, Snoeijer and Andreotti2008), wetting transition was observed to occur at a capillary number well below the value predicted by lubrication theory (Snoeijer et al.
Reference Snoeijer, Andreotti, Delon and Fermigier2007) and numerical simulations (Gao et al.
Reference Gao, Li, Feng, Ding and Lu2016). The underlying mechanisms of these inconsistencies remain unclear. Specifically, the contact-line speed was not well captured in Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018) for
$Ca$
slightly above the critical value, where the present simulations are focused; further experiments in this parameter region would help elucidate whether or not the velocity jump exists in reality.
4.3 Dynamics of gas finger
The dynamics of the gas finger as well as the front of the bubble is characterized by a tip velocity relative to the wall,
$Ca_{b}$
, and a liquid film with thickness
$h_{f}$
. An analysis based on volume conservation reveals that
$$\begin{eqnarray}Ca_{b}=\frac{Ca}{(1-h_{f})^{2}}.\end{eqnarray}$$
The film thickness
$h_{f}$
depends solely on
$Ca_{b}$
. Bretherton (Reference Bretherton1961) first demonstrated that
$h_{f}=1.34Ca_{b}^{2/3}$
, which is only valid for
$Ca_{b}<5\times 10^{-3}$
. For larger
$Ca_{b}$
, Aussillous & Quéré (Reference Aussillous and Quéré2000) proposed an empirical equation through scaling arguments,
$$\begin{eqnarray}h_{f}=\frac{1.34Ca_{b}^{2/3}}{1+2.5\times 1.34Ca_{b}^{2/3}},\end{eqnarray}$$
which fitted the experimental measurements of Taylor (Reference Taylor1961) very well. A similar relation was derived theoretically by Klaseboer, Gupta & Manica (Reference Klaseboer, Gupta and Manica2014), with a modification of the coefficient from 2.5 to
$2.79$
. By varying
$Ca$
and
$\unicode[STIX]{x1D703}_{e}$
, our simulations covered a wide range of
$Ca_{b}$
and
$h_{f}$
, and comparisons with the above expressions are made in figure 5. It is seen in figure 5(a) that the film thickness is best predicted by the expression of Aussillous & Quéré (Reference Aussillous and Quéré2000). Figure 5(b) shows the variation of
$Ca_{b}$
as a function of
$Ca$
, where the continuous curve represents (4.1) with
$h_{f}$
given by (4.2). Once again an excellent agreement between theory and numerical simulations is obtained.
Figure 5. Variations of (a) the film thickness
$h_{f}$
and (b) the bubble velocity
$Ca_{b}$
.
4.4 Bubble length
The length of the bubble,
$L_{b}$
, defined as the distance between the two end points, depends on
$Ca$
and
$\unicode[STIX]{x1D703}_{e}$
. As shown in figure 6(a), a longer bubble can be realized at smaller contact angles and larger driving velocity. A smaller
$\unicode[STIX]{x1D703}_{e}$
corresponds to a slower growth of the rim, enabling a longer time for the gas finger to be elongated; a larger
$Ca$
corresponds to a faster invasion of the gas finger, and again leads to a longer bubble.
An estimation of the bubble length can be provided based on the above scenario. The elongation of the finger occurs at a rate of
$Ca_{b}-Ca^{\ast }$
, and is terminated at the instant
$\unicode[STIX]{x0394}t$
when the rim collapses. Within this period, the liquid volume collected by the dewetting rim is
$V\approx 2\unicode[STIX]{x03C0}h_{f}Ca^{\ast }\unicode[STIX]{x0394}t$
, which is also the rim volume. Here we have assumed that the film thickness is much smaller than the tube radius. On the other hand, the rim collapses when its characteristic length is comparable to the tube radius, corresponding to a volume
$\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}$
(or
$\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}R^{3}$
in dimensional form) with
$\unicode[STIX]{x1D6FC}$
being a fitting coefficient of
$O(1)$
. Therefore, by equating the two forms of the volume, we have
$\unicode[STIX]{x0394}t\approx \unicode[STIX]{x1D6FC}/(2h_{f}Ca^{\ast })$
, and eventually
$$\begin{eqnarray}L_{b}=(Ca_{b}-Ca^{\ast })\unicode[STIX]{x0394}t\approx \frac{\unicode[STIX]{x1D6FC}}{2h_{f}}\left(\frac{Ca_{b}}{Ca^{\ast }}-1\right).\end{eqnarray}$$
Since
$h_{f}$
also depends on
$Ca_{b}$
, equation (4.3) indicates that
$L_{b}$
is a function of
$Ca_{b}$
and
$Ca^{\ast }$
, while
$Ca$
plays an implicit role through (4.1) together with (4.2). The effect of contact angle
$\unicode[STIX]{x1D703}_{e}$
is solely reflected in
$Ca^{\ast }$
.
Figure 6. (a) Variation of the bubble length
$L_{b}$
as a function of
$Ca$
for different contact angles. The lines denote the theoretical prediction (4.3), and the open circles are the simulation results. Inset: comparison of the present theory (line) with the experimental results of Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018) at
$\unicode[STIX]{x1D703}_{e}=68^{\circ }$
(squares). We have fixed
$\unicode[STIX]{x1D6FC}=1$
in (4.3). (b) Variation of the slug length
$L_{s}$
as a function of
$Ca$
.
The theoretical results (4.3) coupled with (4.1) and (4.2) are displayed in figure 6(a) as the solid curves. We adopted the values of
$Ca^{\ast }$
from numerical simulations and simply fixed
$\unicode[STIX]{x1D6FC}=1$
. It is seen that the present theory satisfactorily predicts the bubble length. Due to the large microscopic length used in the present simulations, one may expect the contact-line velocity
$Ca^{\ast }$
to be considerably larger than the experimental value, as also mentioned earlier. According to (4.3), the bubble length
$L_{b}$
would then be significantly underestimated in the present simulations, as compared to the experiments of Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018). We have compared (4.3) with the experimental data at
$\unicode[STIX]{x1D703}_{e}=68^{\circ }$
, by adopting
$Ca^{\ast }\approx 0.014$
extracted from Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018). The present theory reasonably captures the trend of the variation of the bubble length, as shown in the inset of figure 6(a).
4.5 Slug length
The collapse of the rim produces a liquid slug. Due to the continuing accumulation of the film liquid surrounding the leading bubble, the volume of the slug keeps growing until a new film is entrained behind it. The dependence of the eventual slug length
$L_{s}$
on
$Ca$
and
$\unicode[STIX]{x1D703}_{e}$
is displayed in figure 6(b), which exhibits an opposite trend to
$L_{b}$
. For a fixed
$\unicode[STIX]{x1D703}_{e}$
,
$L_{s}$
decreases monotonically with
$Ca$
. As
$Ca$
approaches
$Ca_{cr}$
from above, the evolution of the meniscus towards its equilibrium morphology is characterized by a diverging time scale (Snoeijer et al.
Reference Snoeijer, Andreotti, Delon and Fermigier2007). Therefore,
$L_{s}$
diverges as well. For large
$Ca$
, the bubble length decreases to a constant value that is mildly dependent on
$\unicode[STIX]{x1D703}_{e}$
. In this situation, the slug is roughly a redistribution of the rim fluid, without the accumulation of the film liquid due to the quick deposition of the new film.
5 Conclusion
We have investigated the dynamic wetting transition in a capillary tube using diffuse-interface simulations and theoretical analysis. Our simulations show that the meniscus remains stable for
$Ca<Ca_{cr}$
, and Taylor bubbles are generated beyond the threshold, due to the entrainment of the liquid film and subsequent collapse of the rim near the contact line. Both numerical simulations and lubrication theory predict that the wetting transition occurs at vanishing apparent contact angles. An asymptotic solution of
$Ca_{cr}$
is obtained. At the onset of wetting transition, the contact-line velocity relative to the tube wall drops discontinuously from the average gas velocity displacing the liquid to a smaller dewetting velocity, which is independent of the gas driving velocity. Both the bubble and slug lengths depend, although oppositely, on
$Ca$
and
$\unicode[STIX]{x1D703}_{e}$
. An analytical formula for the bubble length is proposed and fits the present numerical results as well as the experimental measurements of Zhao et al. (Reference Zhao, Pahlavan, Cueto-Felgueroso and Juanes2018).
Acknowledgements
We acknowledge J. Qin and P. Yue for helpful discussions. This work was supported by NSFC (grant nos 11672287, 11422220 and 11621202), the Chinese Academy of Sciences (grant no. KJZD-EW-J01), Youth Innovation Promotion Association CAS and the Strategic Priority Research Program of the Chinese Academy of Sciences (grant no. XDB22040103).
