1. Introduction
The oscillating motion of menisci in thin capillaries is of importance for many applications. One can cite the liquid plugs that obstruct the airways in living organisms for certain pathologies (Baudoin et al. Reference Baudoin, Song, Manneville and Baroud2013), the distribution of fluids in microfluidics (Angeli & Gavriilidis Reference Angeli and Gavriilidis2008) and the oscillations caused by the vapour–liquid mass exchange in heat pipes. This latter application is targeted in the present work as it is relevant to different types of heat pipes. One can cite capillary pumped loops (Zhang, Hou & Sun Reference Zhang, Hou and Sun1998) or loop heat pipes, where the pressure oscillations are observed (Launay et al. Reference Launay, Platel, Dutour and Joly2007) and impact the menisci in the capillary structure. The oscillation of menisci is of special importance for the pulsating heat pipes, also called oscillating heat pipes (Fourgeaud et al. Reference Fourgeaud, Nikolayev, Ercolani, Duplat and Gully2017; Marengo & Nikolayev Reference Marengo and Nikolayev2018; Nikolayev Reference Nikolayev2021), where the liquid films deposited by the oscillating liquid menisci as they recede provide the main channel of the heat and mass transfer. As the film evaporation rate is defined by the local film thickness, one needs to understand the film profile for adequate modelling of the heat pipe. The film evaporation description is the most challenging part because the liquid film can be partially dried out so that triple vapour–liquid–solid contact lines form. Strong heat and mass transfers occur in their vicinity (Janeček & Nikolayev Reference Janeček and Nikolayev2012; Savva, Rednikov & Colinet Reference Savva, Rednikov and Colinet2017) so the contact lines are important to model adequately.
The hydrodynamics of menisci has been extensively studied since the seminal articles of Landau & Levich (Reference Landau and Levich1942), Taylor (Reference Taylor1961) and Bretherton (Reference Bretherton1961). Since their works, a strong effort has been made to understand the dynamics of the Taylor bubbles (i.e. bubbles of length larger than their diameter) and the liquid plugs that separate them. Originally, the hydrodynamics of such a process has been described theoretically within the creeping flow approximation (i.e. for vanishingly small Reynolds numbers) by using the lubrication approach for the liquid film description. The inertial effects have been accounted for by direct numerical simulation (Talimi, Muzychka & Kocabiyik Reference Talimi, Muzychka and Kocabiyik2012). In previous approaches, the liquid film was considered to be continuous, with no dry patches.
The objective of this work is twofold. First, we study the film created by the meniscus oscillation. Second, we want to understand the impact of the film edge, i.e. of the triple contact line, for the simplest case where it is pinned at a surface defect and is thus immobile.
The paper is structured as follows. After an introduction of the model in § 2, the background theory of steady meniscus receding is briefly discussed in § 3. The meniscus oscillation is considered in § 4. The theory is compared to two experimental works involving meniscus oscillation. While the main objective of this paper is to consider the oscillation with no heat and mass transfer, an interesting implication of these results for the film evaporation is discussed in § 5.
2. Model description
Consider a cylindrical capillary tube of an inner radius 
$R$, containing a liquid and its vapour. The vapour–liquid interface is assumed to be axially symmetric (figure 1).
Figure 1. Sketch of the axial cross-section of a cylindrical capillary tube with the liquid film deposited by the receding meniscus. At 
$x=x_{f}$, the film profile solution matches the right branch of the parabola (red solid line), which is a small-slope approximation of the circular meniscus shape (dashed line).
The tube is assumed to be thin enough so the gravity force can be neglected. Following the classical approach (Bretherton Reference Bretherton1961), the vapour–liquid interface can be divided into the film and the meniscus regions. The liquid–vapour interface slope in the film region is assumed to be small so the film can be described with the lubrication theory. The meniscus region is assumed to be controlled by the surface tension only, thus being of constant curvature (shown in figure 1 with a circle of radius 
$R_{m}$).
Because the vapour has a smaller density and viscosity compared with the liquid, the vapour pressure p is assumed to be spatially homogeneous and the vapour-side viscous stress on the interface can be neglected. Under such assumptions, the lubrication theory results in the equation (Nikolayev Reference Nikolayev2010) that describes the interface dynamics in the film region,
\begin{equation} \frac{\partial h}{\partial t}+\frac{\partial}{\partial x}\left( \frac{h^3}{3\mu}\frac{\partial{\rm \Delta} p}{\partial x} \right) ={-}\frac{J}{\rho}, \end{equation}
where 
$h=h(x,t)$ is the local film thickness, and 
$J$ is the mass flux across the interface, defined to be positive at evaporation. Here, 
$\mu$ and 
$\rho$ are the liquid shear viscosity and density, respectively. The lubrication theory is applicable within the assumption 
$|\partial h/\partial x|\ll 1$. The pressure jump 
${\rm \Delta} p=p-p_{l}$ (with 
$p_{l}$, the liquid pressure) across the interface obeys the Laplace equation
\begin{equation} {\rm \Delta} p=\sigma \left( K +\frac{1}{R - h} \right)\simeq \sigma\left( \frac{{\partial ^2}h}{\partial {x^2}}+\frac{1}{R_{m}}\right), \end{equation}
where 
$\sigma$ is the surface tension, 
$K$ is the two-dimensional interface curvature in the axial cross-section shown in figure 1; in the small-slope approximation, 
$K\simeq \partial ^2h/\partial x^2$. Because of this limitation, such a film region theory is not able to describe the meniscus region. The radial contribution to the curvature 
$R_{m}^{-1}$ is assumed to be independent of 
$x$ in the film region (where 
$h$ is much smaller than 
$R$). At a large 
$x$, (2.1) results in an increasingly larger 
$h$, where the viscous forces vanish and the surface tension alone controls the interface, so its curvature 
$\partial ^2h/\partial x^2$ becomes constant. Such a condition corresponds to a parabolic shape in the axial plane. This parabola needs to be joined to the circular meniscus, which results in the condition
\begin{equation} \partial^2 h/\partial x^2|_{x=x_{f}}=R_{m}^{{-}1}, \end{equation}
defined at the ending point 
$(x_{f}, h_{f})$ of the film region.
3. Isothermal problem with infinite film: steady solutions
First, we would like to recall the theory for the case with no phase change (
$J=0$), so (2.1) becomes
\begin{equation} \frac{\partial h}{\partial t} + \frac{\sigma}{3\mu} \frac{\partial}{\partial x}\left( h^3\frac{\partial ^3h}{\partial x^3} \right) =0. \end{equation}
The velocity of the axial meniscus centre (assumed positive for a receding meniscus according to the 
$x$-axis direction choice) is denoted 
$U$. The contact line is not considered so the film is infinite. In this case, it is advantageous to choose the frame of reference linked to the moving meniscus where the axial coordinate becomes 
$x'=x-Ut$. Equation (3.1) can then be rewritten as
\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial}{\partial x'}\left( \frac{\sigma }{3\mu}h^3 \frac{\partial ^3h}{\partial {x'}^3} -Uh \right) =0. \end{equation}
This equation (and this frame of reference) is convenient to use in the present section and in Appendix A because the film is infinite there. However, the prime will be dropped to make notation less heavy. In all other sections, where the contact line is considered, the wall frame of reference will be used as it is more convenient, so the meaning of 
$x$ will be as in (2.1) and (3.1).
The steady version of (3.2) for the case of a constant positive velocity 
$U_{r}$ (that should be used for 
$U$ in this case) is the Landau–Levich equation (Landau & Levich Reference Landau and Levich1942) describing the flat infinite film being deposited by a receding meniscus. Note that the film in a cylindrical capillary (Bretherton Reference Bretherton1961) is described by the same equation because of the approximation (2.2).
The boundary conditions at 
$x\to -\infty$ describe the flat film of the yet unknown thickness 
$h_{r}$,
\begin{equation} h=h_{r},\quad \partial h/\partial x=0. \end{equation} The scaling of this problem (table 1) is based on 
$h_{r}$. The characteristic axial length scale involves the capillary number 
$Ca_{r}=\mu U_{r}/\sigma$ and is chosen in such a way that the dimensionless equation
\begin{equation} \frac{\partial }{\partial \tilde x}\left( { \frac{{\tilde h}^3}{3}\frac{{{\partial ^3}\tilde h}}{{\partial {\tilde x ^3}}}} - \tilde h \right) =0 \end{equation}does not contain any constants; the tilde means hereafter the corresponding dimensionless variable.
Table 1. Dimensional reference values that are used to make the governing equations non-dimensional.
By integrating (3.4) from 
$-\infty$ to 
$\tilde x$ and using the conditions (3.3a,b) (in the dimensionless form, 
$\tilde h(\tilde x\to -\infty )=1$), one obtains
\begin{equation} \frac{\partial^3\tilde h}{\partial\tilde x^3} =3\frac{\tilde h-1}{\tilde h^3}, \end{equation}which is equivalent to the Bretherton (Reference Bretherton1961) equation within a factor 3 that we leave in the equation instead of putting inside the scaling parameters of table 1. This helps us to avoid it in many formulas used below.
Consider now the behaviour at large 
$\tilde x$. Equation (3.5) remains valid until the transition region (between the film and the meniscus), in which 
$1\ll \tilde h < R$. From (3.5), 
${\partial ^3\tilde h}/{\partial \tilde x^3} \simeq 0$. This means that, at large 
$\tilde x$, the second derivative is finite, which is compatible to the condition (2.3). Equation (3.5) can be integrated numerically (see Nikolayev & Sundararaj (Reference Nikolayev and Sundararaj2014) for details). The ending point 
$(\tilde x_{f}, \tilde h_{f})$ of the integration interval is chosen from a condition that 
$\partial ^2\tilde h/\partial \tilde x^2$ reaches a plateau (with a required accuracy). Such a calculation results in a numerical value for this plateau
\begin{equation} \left.\frac{\partial^2\tilde h}{\partial\tilde x^2}\right|_{\tilde x=\tilde x_{f}}=\alpha\simeq 1.3375. \end{equation}
The resulting profile 
$\tilde h(\tilde x)$ can be found in figure 13(a) of Appendix A. Note that 
$\alpha$ is equivalent to the numerical value 0.643 originally found by Bretherton (Reference Bretherton1961); 
$\alpha \simeq 0.643\cdot 3^{2/3}$, where the factor appears because of the different scaling. By comparing (2.3) and (3.6), one finds the expression for the film thickness
\begin{equation} h_{r}=\alpha R_{m}Ca_{r}^{2/3}. \end{equation} The only yet unknown quantity is 
$R_{m}$ that can be found as proposed by Klaseboer, Gupta & Manica (Reference Klaseboer, Gupta and Manica2014). As mentioned before, near the point 
$(\tilde x_{f}, \tilde h_{f})$ the film shape should satisfy the condition (2.3) which means that at 
$h\gg h_{r}$ it asymptotically approaches a parabola
\begin{equation} y=(x-x_{s})^2/(2R_{m})+h_{s}, \end{equation}
where the parameters 
$(x_{s}, h_{s})$ are yet to be determined. To understand their meaning, one recalls that, near its minimum where its curvature is 
$R_{m}^{-1}$, the parabola approximates the circular meniscus profile
\begin{equation} (x-x_{s})^2+(y-h_{s}-R_{m})^2=R_{m}^2. \end{equation}
It is evident now that 
$(x_{s}, h_{s})$ is the circle lowest point. One can obtain 
$(x_{s}, h_{s})$ by fitting the film profile 
$\tilde h(\tilde x)$ near the point 
$(\tilde x_{f}, \tilde h_{f})$ to a parabola. Klaseboer et al. (Reference Klaseboer, Gupta and Manica2014) report that the dimensionless value of 
$\tilde h_{s}$ slightly grows with 
$\tilde h_{f}$. Bretherton (Reference Bretherton1961) gives 
$\tilde h_{s}=2.79$. The asymptotic value 
$\tilde h_{s}=2.90$ is obtained for 
$\tilde h_{f}> 10^6$. However, to obtain a continuous overall interface profile, the matching point 
$(\tilde x_{f}, \tilde h_{f})$ film–parabola should be lower than the point where the parabola–circle transition occurs. This requires 
$h_{f}< R$. For the value 
$\tilde h_{f}\simeq 50$ that satisfies this condition in practical situations, Klaseboer et al. (Reference Klaseboer, Gupta and Manica2014) find 
$\tilde h_{s}=2.5$, which is also the value found from the experimental data fits, as discussed below. In summary, the 
$\tilde h_{s}$ variation is weak and one can consider that the circle is nearly invariant of the specific 
$\tilde h_{f}$ choice. Klaseboer et al. (Reference Klaseboer, Gupta and Manica2014) have proposed the equation
\begin{equation} R_{m}+h_{s}=R,\end{equation}
that centres the circle with respect to the tube and thus links 
$R_{m}$ to 
$R$, cf. figure 1. By using (3.7) in this equation, one finally obtains
\begin{equation} R_{m}= \frac{R}{1 + \alpha \tilde h_{s} Ca_{r}^{2/3}}. \end{equation}By combining (3.7) and (3.11), one can now finalise the film thickness expression
\begin{equation} h_{r}=\frac{\alpha R Ca_{r}^{2/3}}{1 + \alpha \tilde h_{s} Ca_{r}^{2/3}}. \end{equation} The value 
$\tilde h_{s}\simeq 2.5$ has been determined by Aussillous & Quéré (Reference Aussillous and Quéré2000) from the experimental data fits. For 
$Ca_{r}\to 0$, 
$R_{m}\simeq R$ and (3.12) reduces to Bretherton's (Reference Bretherton1961) original expression
\begin{equation} h_{r}=\alpha R Ca_{r}^{2/3}. \end{equation} One can consider the meniscus advancing at a constant velocity 
$U=-U_{a}$ (where 
$U_{a}$ is the modulus of the advancing velocity) over the pre-existing film of thickness 
$h_{r}$. Such a motion has been understood as well. It has been shown (Bretherton Reference Bretherton1961) that the film has a wavy shape (ripples) near the meniscus, cf. the solid curve in figure 13(a) of Appendix A. The wavelength of ripples depends on the ratio 
$U_{a}/U_{r}$ (Maleki et al. Reference Maleki, Reyssat, Restagno, Quéré and Clanet2011; Nikolayev & Sundararaj Reference Nikolayev and Sundararaj2014), where 
$U_{r}$ can be deduced (with (3.13)) from 
$h_{r}$. The meniscus radius for the steady advancing case was determined with (3.10) by Cherukumudi et al. (Reference Cherukumudi, Klaseboer, Khan and Manica2015).
4. Oscillations in the presence of a contact line
A previous study (Nikolayev & Sundararaj Reference Nikolayev and Sundararaj2014) demonstrates the film behaviour for the case of an infinite film. There is no physical criterion imposing its thickness so it is another independent parameter. When the meniscus approaches the leftmost (in the reference of figure 1) position observed during oscillations, the ripples created near the advancing meniscus propagate over the film to infinity, so there is no possible periodical regime. This propagation is amplified by the discrepancy between the imposed film thickness and the film thickness (3.7) defined by the receding meniscus velocity, which is zero at the leftmost point. Thus a discrepancy exists for any imposed film thickness. In practical situations of oscillating motion (Fourgeaud et al. Reference Fourgeaud, Ercolani, Duplat, Gully and Nikolayev2016; Rao et al. Reference Rao, Lefèvre, Czujko, Khandekar and Bonjour2017), the contact line appears because of the film evaporation caused by the tube wall heating. The film completely vaporises beyond the leftmost meniscus position. Before addressing the heating case, in this section we discuss the film shape in the presence of contact line without any heating.
The hydrodynamics of the pinned (static) contact line is simpler than the dynamic case. For this reason one needs to understand it first. This is a purpose of this work. The contact line pinning often occurs in capillaries (Mohammadi & Sharp Reference Mohammadi and Sharp2015). It is caused by the wall heterogeneity (surface defects) that can be either chemical or geometrical (surface roughness). The heterogeneity can be modelled as a spatial variation of surface energy. The result of such a theory (Iliev, Pesheva & Nikolayev Reference Iliev, Pesheva and Nikolayev2014) is that the microscopic contact angle averaged along the contact line can vary between the static advancing 
$\theta _{adv}$ and the static receding 
$\theta _{rec}$ angles while the contact line remains immobile. In our calculation, 
$\theta _{adv} - \theta _{rec}$ (called wetting hysteresis) is assumed to be sufficiently large so the contact line always remains immobile. In experiments, the hysteresis can be as large as 
$50^\circ$ (de Gennes Reference de Gennes1985), which is larger than the angle oscillation magnitude considered below.
At oscillations with the fixed contact line, there are no vortices near it (Ting & Perlin Reference Ting and Perlin1987) and the flow is known to be well described by the lubrication approximation, even down to the nanometric scale (Mortagne et al. Reference Mortagne, Lippera, Tordjeman, Benzaquen and Ondarçuhu2017).
4.1. Relaxing the pressure divergence by the Kelvin effect
The Stokes problem of the straight wedge with a varying opening angle leads to the logarithmic pressure divergence, cf. Appendix B. Such a divergence is integrable and thus does not cause a paradox similar to that of the moving contact line. However, the infinite pressure is non-physical. In addition, the pressure boundary condition at the contact line would be difficult to use in calculation because it requires prior knowledge of the contact angle and its time derivative (cf. (B8)) that need to be determined themselves during the solution procedure. As we consider volatile fluids, the phase change together with the Kelvin effect are introduced. The latter makes the pressure finite everywhere, as shown below. The problem is formulated here for a general case where the tube wall can be superheated or subcooled with respect to the saturation temperature 
$T_{sat}$ corresponding to the imposed vapour pressure 
$p$. The wall superheating is denoted 
${\rm \Delta} T$. The tube wall temperature is thus 
$T_{w}=T_{sat}+{\rm \Delta} T$.
Conventional hypotheses concerning the liquid film mass exchange (Nikolayev Reference Nikolayev2010) are applied. A linear temperature profile in the radial direction is assumed in the thin liquid film, so the energy balance at the interface results in the mass flux
\begin{equation} J = \frac{k(T_{w}-T_{int})}{h\mathcal{L}}, \end{equation}
where 
$T_{int}$ is the temperature of the vapour–liquid interface, 
$k$ is the liquid heat conductivity and 
$\mathcal {L}$ is the latent heat. The film is assumed here to be thin with respect to 
$R$ so the one-dimensional conduction description applies. The evaporation impact on a film is twofold. First, the film thickness decreases with time everywhere along the film, which is described by the balance of the first and the right-hand side term of (2.1). Generally, the film thinning is not strong during an oscillation period (Fourgeaud et al. Reference Fourgeaud, Nikolayev, Ercolani, Duplat and Gully2017).
We focus here on the second effect that appears because of the strength of evaporation in a narrow vicinity of contact line. If the vapour–liquid interface was at a fixed saturation temperature (
$T_{int}=T_{sat}$), the mass flux 
$J$ (4.1) would diverge at the contact line 
$h=0$ as 
$J\sim {\rm \Delta} T/h$, which is non-physical because total evaporated mass (the integral of 
$J$) would be infinite.
The Kelvin effect, i.e. the dependence of 
$T_{int}$ on the interfacial pressure jump 
${\rm \Delta} p$
\begin{equation} T_{int} = T_{sat}\left( {1 + \frac{{\rm \Delta} p}{\mathcal{L} \rho}} \right), \end{equation}
can relax the singularity (Janeček & Nikolayev Reference Janeček and Nikolayev2012), because it allows 
$T_{int}$ to vary along the interface so it can attain the wall temperature 
$T_{w}$ at the contact line so the mass flux
\begin{equation} J(x\to 0)=0. \end{equation}From the temperature continuity, one obtains the condition
\begin{equation} {\rm \Delta} p(x \to 0) ={\rm \Delta} p_{cl}, \end{equation}where the pressure jump at the contact line is
\begin{equation} {\rm \Delta} p_{cl}=\frac{{\mathcal{L}}{\rho}}{T_{sat}} {\rm \Delta} T. \end{equation}One can show that a solution that satisfies this condition can indeed be found (cf. Appendix C).
Equations (4.1) and (4.2) result in
\begin{equation} J = \frac{k}{h\mathcal{L}}\left({\rm \Delta} T- {\rm \Delta} p\frac{ T_{sat}}{\mathcal{L} \rho} \right). \end{equation}With its substitution into (2.1), the governing equation becomes
\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial }{\partial x}\left(\frac{h^3}{3\mu} \frac{\partial{\rm \Delta} p}{\partial x}\right)=\frac{{\rm \Delta} p-{\rm \Delta} p_{cl}}{h}\frac{kT_{sat}}{ ({\mathcal{L}} \rho)^2}. \end{equation}
The problem is now regular (because 
${\rm \Delta} p$ is not divergent anymore), unlike other microscopic approaches (Savva et al. Reference Savva, Rednikov and Colinet2017). As the Kelvin effect alone is capable of relaxing the contact line singularity, the other microscopic scale effects such as hydrodynamic slip, Marangoni effect or interfacial kinetic resistance (Janeček & Nikolayev Reference Janeček and Nikolayev2012) are not crucial anymore. They are not included in our model for the sake of clarity.
The characteristic size of the contact line vicinity where the Kelvin effect is important is 
$\ell _{K}$ (C5), cf. Appendix C for more details. It is nanometric (Janeček et al. Reference Janeček, Andreotti, Pražák, Bárta and Nikolayev2013) and is thus significantly smaller than the characteristic scale of film shape variation that we call macroscopic. For this reason, (4.7) can be understood within a multi-scale paradigm in the spirit of the asymptotic matching techniques (Janeček et al. Reference Janeček, Andreotti, Pražák, Bárta and Nikolayev2013). In the inner region, commonly called the microregion, the first (transient) term is negligible with respect to the Kelvin term (
${\rm \Delta} p$ containing term in the right-hand side). The problem is reduced to that of Appendix C.2. To summarise it, when 
${\rm \Delta} T\neq 0$, a strong interfacial curvature that exists in the microregion can cause a difference between the microscopic contact angle 
$\theta _{micro}$ and the interface slope 
$\theta$ defined at 
$x\to \infty$ within the microregion. In the outer (macroscopic) region, the Kelvin effect is negligible so the 
${\rm \Delta} p$ term on the right-hand side of (4.7) vanishes. For 
${\rm \Delta} T=0$, this equation is that of the isothermal problem (3.1). Both problems can be matched at scale 
$x_{meso}\gg \ell _{K}$, much smaller than the macroscopic scale. The apparent contact angle visible on this latter scale is thus equal to 
$\theta$. It is assumed hereafter that the pinning occurs at a length scale smaller than 
$\ell _{K}$, i.e. there are nanometric defects with sharp borders on which the contact line is pinned.
In § 4, a globally isothermal problem is considered, 
$T_{w}=T_{sat}$, so 
${\rm \Delta} T=0$ and 
${\rm \Delta} p_{cl}=0$, which means 
$\theta =\theta _{micro}$. At such conditions, the mass exchange appearing at the macroscopic scale is very weak so that it can be safely neglected. This does not mean, however, that the mass exchange is absent in the microregion where 
${\rm \Delta} p$ can be large, as mentioned above. The mass flux 
$J$ scales with 
${\rm \Delta} p$ according to (4.6), so phase change occurs. The situation here shares certain similarities with the contact line motion paradox solved by the Kelvin effect (Janeček et al. Reference Janeček, Andreotti, Pražák, Bárta and Nikolayev2013). Consider e.g. an increasing in time 
$\theta$. According to (C7), 
${\rm \Delta} p>0$ in the very contact line vicinity so the condensation occurs there. It is compensated exactly by evaporation farther away from the contact line, so the net mass exchange is zero. It should be noted that the fluid flow associated with the phase change is strongly localised within a nanoscale distance from the contact line comparable to 
$\ell _{K}$.
In conclusion, all results obtained in the present § 4, can be seen as obtained with (3.1) because the microregion details cannot be resolved at the macroscale pictured in the figures below. However, numerical calculations of the regularised equation (4.7) are carried out in reality.
4.2. Oscillation problem statement
The meniscus now oscillates, and the position 
$x_{m}$ of its centre (which is the experimentally measurable quantity) travels periodically with a period 
$P$ and an amplitude 
$A$. One can assume its harmonic oscillation
\begin{equation} x_{m} (t)= x_{i} + A[ 1 - \cos(2{\rm \pi} t/P)], \end{equation}
where 
$x_{i}$ is the initial meniscus centre position. Alternatively, one can take the experimentally measured dependence 
$x_{m}(t)$ while comparing the data with the experiment (cf. § 4.10 below). The contact line is pinned at the position 
$x=0$, and the contact angle 
$\theta$ varies. For the harmonic oscillation case, the meniscus velocity is 
$U(t)= U_0\sin (2{\rm \pi} t/P)$, where the velocity amplitude 
$U_0=2{\rm \pi} A/P$ is convenient to choose as the characteristic velocity to define the capillary number 
$Ca_0=\mu U_0/\sigma$ and to make all the quantities dimensionless (cf. table 1).
Because of the fixed contact line, the frame of reference of the tube wall is chosen. Equation (4.7) (with the substitution of (2.2)) is solved for 
$x\in [0, x_{f}]$. The length 
$x_{f}$ is imposed as explained in §§ 4.3 and 4.5 below.
The boundary conditions are defined as
$$\begin{gather} h(x = 0) = 0, \end{gather}$$
$$\begin{gather}\left.\frac{\partial {\rm \Delta} p}{\partial x} \right|_{ x \to 0} = 0, \end{gather}$$
$$\begin{gather}{\rm \Delta} p( x = x_{f}) = R_{m}^{{-}1}, \end{gather}$$
$$\begin{gather}h( x = x_{f}) = h_{f}, \end{gather}$$
where 
$R_{m}$ and 
$h_{f}$ are discussed in § 4.3. The condition (4.9a) is a geometrical constraint at the contact line. Equation (4.9b) is a weaker form of the condition (4.4) used to provide numerical stability. The boundary conditions (4.9c) and (4.9d) impose the liquid film curvature and thickness at the right end of the integration interval for each time moment.
4.3. Determination of the meniscus curvature
For a small film thickness, one can assume that the meniscus radius 
$R_{m}$ is constant and equal to 
$R$ during oscillation. It is actually a good approximation for a small 
$Ca_0\lesssim 10^{-3}$. However, at a larger 
$Ca_0$, the film thickness impacts 
$R_{m}$ (see § 4.8 below). Since the film thickness depends on the meniscus velocity, so does the meniscus radius 
$R_{m}$, cf. § 3. Therefore, 
$R_{m}$ varies in time. In this section, we generalise to any meniscus dynamics the method for 
$R_{m}$ determination (Klaseboer et al. Reference Klaseboer, Gupta and Manica2014) discussed above for the steady receding case.
Similarly to the steady case of § 3, one needs first to match the film shape 
$h(x)$ to the parabola (3.8), and then the parabola to a circle (3.9). The matching between the film and the parabola means both the continuity and the smoothness (equality of the derivatives)
$$\begin{gather} h_{f}=(x_{f}-x_{s})^2/(2R_{m})+h_{s}, \end{gather}$$
$$\begin{gather}\left.\frac{\partial h}{\partial x} \right|_{ x= x_{f}}= \frac{x_{f}- x_{s}}{R_{m}}, \end{gather}$$
where the left-hand sides come from the film calculation and all the parabola parameters are time-dependent. As in the approach of Klaseboer et al. (Reference Klaseboer, Gupta and Manica2014), (3.10) serves to find 
$R_{m}$. We introduce in addition a relationship of the abscissas of the lowest and rightmost points of a circle that is needed to define 
$x_{s}$ (figure 1)
\begin{equation} R_{m}+ x_{s}=x_{m}.\end{equation} In the present algorithm, 
$x_{f}$ imposed to such a value that the difference 
$x_{m}-x_{f}$ does not vary in time and 
$h_{f}=h(x_{f})$ remains large with respect to the deposited film thickness. As discussed in § 3, the solution is nearly independent of the specific choice of 
$h_{f}$ (and thus of 
$x_{f}$). The set of (2.2), (3.10) and (4.7)–(4.10) is then complete, so the film shape and the unknown parameters (
$R_{m},x_{s}, h_{s}, h_{f}$) can be determined for each 
$t$.
4.4. Initial conditions and solution periodicity
One needs to define now the initial film shape 
$h(x,0)$ at 
$t =0$, which corresponds to the (yet unspecified) leftmost meniscus position 
$x_{i}$ according to (4.8). As an initial film profile, we choose that of equilibrium satisfying the condition 
$\partial h/\partial t=0$ that can be used in (3.1). From the boundary condition (4.9b), one finds straightforwardly 
$\partial ^3 h/\partial x^3=0$, i.e. the parabolic shape
\begin{equation} h(x,0)=\frac{x^2}{2R_{m,i}}+\theta_{i} x, \end{equation}
where 
$\theta _{i}\equiv \theta (t=0)$ is the initial contact angle. It serves as another boundary condition, additional to (4.9a) and (4.9c). By applying (3.10) and (4.10) at 
$t=0$, one gets
$$\begin{gather} R_{m,i}\equiv R_{m}(t=0)=R/(1-\theta_{i}^2/2), \end{gather}$$
$$\begin{gather}x_{i}=R_{m,i}(1-\theta_{i}). \end{gather}$$
These expressions are the small-angle approximations of the expressions 
$R_{m,i}=R/\cos \theta _{i}$ and 
$x_{i}=R_{m,i}(1-\sin \theta _{i})$ because (4.11a) is an approximation of the initially spherical meniscus.
During oscillation, the liquid is driven by the meniscus motion and the free interface remains in the state where 
$\partial h/ \partial t$ is always balanced by the curvature gradient, more precisely, by the second term of (3.1). This occurs because there are no other forces, in particular, no inertia. When the meniscus comes near the leftmost position, 
$U$ decreases and the system approaches the state (4.11) with no curvature gradient, thus, 
$| \partial h/ \partial t|$ decreases to zero or almost zero. It is not a rigorous proof that the state (4.11) belongs to the limit cycle of the system, albeit, it should be quite close to it. This is surely true when the relaxation time 
$t_{rel} \ll P$, which is our case (cf. Appendix A for 
$t_{rel}$ discussion). Indeed, the numerical simulations show that 
$h(x,P)$ is indistinguishable from 
$h(x,0)$, cf. figure 2 below. So do all other parameters (curvature, contact angle, etc.). This finding allows us to simulate a unique period.
Figure 2. Periodic interface shape variation at oscillation for 
$\tilde P=50$, 
$Ca_0=10^{-3}$ and 
$\theta _{i}=20^\circ$. The labels give the times corresponding to each profile, and arrows indicate the meniscus motion direction.
4.5. Numerical implementation
The scales for the main quantities to make them dimensionless are shown in table 1. With such a ‘natural’ scaling three main dimensionless parameters are left: 
$\theta _{i}$, 
$\tilde P$ and 
$Ca_0$. All quantities will be studied in this parametric space. The dimensionless amplitude is linked to the period
\begin{equation} \tilde A= \tilde P/(2{\rm \pi} ), \end{equation}where a dimensionless quantity is denoted with a tilde. There is one more dimensionless parameter
\begin{equation} N=\frac{\mu k T_{sat}}{{\mathcal{L}}\rho\alpha R \,Ca_0} \end{equation}
that describes the magnitude of the Kelvin effect in the microregion. However, it does not impact the interface shape at the film scale provided the characteristic microscopic scale (C5) is chosen to be small enough (cf. § 4.1). The mesh size is exponentially refined near the contact line (as 
$\tilde {x}\to 0$) to capture the contact angle variation without considerably increasing the total number of nodes (Nikolayev Reference Nikolayev2010).
At 
$\tilde t=0$, we choose a value of 
$\tilde x_{f}=10$, for which 
$\tilde h_{f}$ is around 50, cf. the discussion in § 3. At 
$\tilde t>0$, the difference 
$\tilde x_{m}-\tilde x_{f}$ is maintained constant and equal to that at 
$\tilde t=0$. To avoid the discretisation error for the contact angle, the initial interface profile is determined numerically by solving the equilibrium version of (4.7), instead of using the analytical profile (4.11a).
Equation (4.7) is solved numerically with the finite volume method (FVM), which is more stable numerically (Patankar Reference Patankar1980) than a more conventional finite difference method. In one dimension, a finite volume is just a segment. The variables such as 
$h$ and their even-order derivatives are defined at its centre (called a node), while the odd-order derivatives are defined at the segment ends. The FVM has the advantage that the liquid flux is continuous at the segment ends. Nonlinear terms are managed by iteration: they include values from the previous iteration. The numerical algorithm is similar to that used by Nikolayev (Reference Nikolayev2010).
One is interested in amplitudes that are large with respect to the meniscus width 
$x_{m}-x_{s}$, which means large dimensionless periods of oscillation, see (4.12). This signifies that the computational domain size varies considerably during oscillations. The grid thus needs to be adaptive and the calculation time can be of the order of a day on a regular personal computer.
4.6. Interface profile during oscillation
Figure 2 shows the interface profiles at several time moments during oscillation. The meniscus motion follows the harmonic law (4.8). The liquid film is deposited until 
$t=P/2$. For 
$t>P/2$, the meniscus advances over the deposited film. The ripples near the meniscus appear, like during the steady meniscus advance discussed in § 3. The interface profiles 
$\tilde h(\tilde x,0)$ and 
$\tilde h(\tilde x,\tilde P)$ are indistinguishable, which confirms the periodicity of the oscillations.
Figure 3 shows examples of interface profiles at 
$t=P/2$ (when the meniscus is at its rightmost position so the film length attains its maximum) for several values of the initial contact angle 
$\theta _{i}$. All the film profiles are presented in the meniscus reference
\begin{equation} x'=x-2A-x_{i}, \end{equation}
meaning that the meniscus centre is at 
$x'=0$, cf. (4.8). One can see that the interface shape near the meniscus is independent of 
$\theta _{i}$. During the meniscus recession, the film loses information about the contact line, so the film shape near the meniscus is controlled by the meniscus dynamics only. This is not surprising as the flow in the film is expected to occur only in the contact line and meniscus vicinities, but not in the middle of the film. The film profiles exhibit a local minimum 
$h_{min}$ discussed in § 4.7 below. It appears because of the meniscus velocity reduction at the end of a half-period.
Figure 3. Film shapes in the meniscus centre reference at 
$t=P/2$ for 
$Ca_0=10^{-3}$ and different 
$\theta _{i}$ and 
$P$. The quasi-steady profiles 
$\tilde h_{q}$ discussed in § 4.9 are shown for comparison. The scaled contact angle 
$\tilde {\theta }$ defined as 
$\tan \tilde {\theta }=Ca_0^{-1/3}\tan \theta$ is indicated; (a) 
$\tilde P=50$ and (b) 
$\tilde P=150$.
4.7. Contact angle during oscillation
A typical variation of 
$\theta$ during oscillation is plotted in figure 4(a). The initial contact angle 
$\theta _{i}$ is the maximum contact angle achieved during the periodic motion.
Figure 4. Contact angle variation for 
$Ca_0=10^{-3}$, 
$\tilde P=50$. (a) Variation of contact angle during an oscillation for different values of 
$\theta _{i}$. The inset shows enlarged undulating portions of the curves. (b) The interface shape variation near the contact line around 
$t=0.913P$ where 
$\theta _{min}$ is attained for 
$\theta _{i}=20^\circ$ (cf. the inset to figure 4a).
In the beginning of a period, the capillary forces lead to the fast contact angle reduction until the meniscus recedes far enough so the curvature gradient reduces and the contact angle becomes nearly constant for a large part of a period. This nearly constant value is quite insensitive to both 
$\theta _{i}$ and 
$P$. A small ridge of constant curvature forms near the contact line. This phenomenon is similar to the dewetting ridge but of much smaller magnitude because the contact line is pinned. A small ridge can be seen in figure 3(b) in the contact line region for the curve corresponding to 
$\theta _{i}=40^\circ$.
The ridge width slowly grows so 
$\theta$ slowly decreases until the ripples in the near-meniscus region approach the contact line during the backward stroke (figure 4b) at the end of a period. This causes the contact angle oscillations, during which its minimal value 
$\theta _{min}$ is attained
\begin{equation} \theta(t) \geq \theta_{min}. \end{equation}
It depends quite weakly on 
$\theta _{i}$ (figure 4a). The variation of 
$\theta _{min}$ with the system parameters follows the variation of 
$h_{min}$ (which is a minimum of 
$h(x, 0.5P)$ observed near the meniscus). This is illustrated in figure 5, where the variations of 
$h_{min}$ and 
$\theta _{min}$ with the system parameters are compared. Only the variations with 
$\tilde P$ and 
$Ca_0$ are considered (as mentioned above, the dependence on 
$\theta _{i}$ is quite weak). Evidently, the variations of 
$h_{min}$ and 
$\theta _{min}$ with the system parameters are similar, which shows their intrinsic link. At 
$Ca_0\to 0$, the dependence on 
$Ca_0$ is weak, but becomes stronger at large 
$Ca_0$. This minimal value of the contact angle is of importance (cf. § 5 below). Since the motion is periodic, the contact angle 
$\theta _{i}$ is attained at 
$t=P$.
Figure 5. Variations of 
$h_{min}$ and 
$\theta _{min}$ with 
$Ca_0$ for 
$\theta _{i}=20^\circ$ and different 
$\tilde P$; (a) 
$h_{min}$ variation with 
$Ca_0$ and (b) 
$\theta _{min}$ variation with 
$Ca_0$.
4.8. Meniscus curvature during oscillation
The meniscus curvature (figure 6) changes periodically during oscillations. The value of 
$R_{m}$ can be compared to the quasi-steady value 
$R_{m, q}$ given by (3.11) where 
$\tilde h_{s}=2.5$ and 
$Ca_r$ are calculated with the instantaneous meniscus velocity 
$U = U_0\sin (2{\rm \pi} \tilde t/\tilde P)$ during the liquid recession (
$t\leq 0.5P$).
Figure 6. Time evolution of the radius of meniscus curvature during an oscillation. The quasi-steady evolution of 
$R_{m, q}$ at receding (
$t\leq 0.5P$) is shown for comparison: (a) 
$R_{m}$ evolution for different 
$\tilde P$; 
$\theta _{i}=40^\circ$ and 
$Ca_0=10^{-3}$ are fixed, (b) 
$R_{m}$ evolution for different 
$\theta _{i}$; 
$\tilde P =150$ and 
$Ca_0=10^{-3}$ are fixed and (c) 
$R_{m}$ evolution (dashed curves) for different 
$Ca_0$; 
$\tilde P =150$ and 
$\theta _{i}=20^\circ$ are fixed. The solid curves of the respective colour show 
$R_{m, q}(t)$.
For 
$t=0$, 
$R_{m}$ is defined with (4.11b), which differs from the quasi-steady value 
$R_{m, q}=R$ for 
$U=0$. This difference occurs because of the contact line presence. In its absence (pre-wetted tube, the situation equivalent to the limit 
$\theta _{i}\to 0$ in our model), the initial radius would be close to 
$R$ because the wetting film is much thinner than the film considered here.
Within the time scale 
$\sim 5t_{rel}$ (see Appendix A), 
$R_{m}$ relaxes to the quasi-steady value 
$R_{m, q}$. The curvature 
$R_{m}$ remains close to 
$R_{m, q}$ until the deceleration that occurs near the rightmost meniscus position (at 
$t=0.5P$, where 
$U=0$). However, the shape relaxation causes a delay, so the inequality 
$R_{m}< R_{m, q}(U=0)=R$ always holds at the point 
$t=0.5P$. During the backstroke, 
$R_{m}$ varies, finally attaining the initial value (4.11b) that depends only on 
$\theta _{i}$ (figure 6b). Evidently, the amplitude of 
$R_{m}$ oscillation grows with 
$Ca_0$ (figure 6c), as foreseen by (3.11). One also mentions the non-monotonic 
$R_{m}$ variation during the backstroke with a local minimum around 
$t\simeq 0.8-0.9P$ (figure 6c). This minimum appears when the trough of film ripples approaches the contact line close enough and is thus correlated with the contact angle minimum.
4.9. Quasi-steady approach and average film thickness
One can see that the film is thickest in its centre (figure 3), which correlates with the maximum of the meniscus velocity. It is thus interesting to compare the film thickness with its quasi-steady value. Within the quasi-steady approach, the term 
$\partial h/\partial t$ is neglected and the quasi-steady thickness 
$h_{q}(x)$ can be defined as corresponding to the meniscus receding velocity as if it was constant at each time moment (Fourgeaud et al. Reference Fourgeaud, Nikolayev, Ercolani, Duplat and Gully2017; Youn, Han & Shikazono Reference Youn, Han and Shikazono2018). In § 4.8 it has been shown that a simple quasi-steady approach predicts well the meniscus curvature for 
$0< t<0.5P$. It is interesting to see if it is efficient in predicting the film profile. In this section, only the profile at 
$t=0.5P$ is considered, i.e. that with the longest film.
A difficulty appears because the film thickness 
$h(x,0.5P)$ depends on 
$x$. It is clear that 
$h$ depends on the velocity that the meniscus had at film deposition, but at which time moment? The most obvious first option is a moment 
$t$ when the meniscus centre was at the point 
$x$ (i.e. 
$x_{m}(t)=x$). A more sophisticated option is a moment 
$t'>t$ such that
\begin{equation} x_{m}(t')=x+{\rm \Delta} x, \end{equation}
with 
${\rm \Delta} x>0$. In the previous approaches (Fourgeaud et al. Reference Fourgeaud, Nikolayev, Ercolani, Duplat and Gully2017; Youn et al. Reference Youn, Han and Shikazono2018), the first option was used. By using (4.8) one easily finds that this assumption defines 
$U_{q}(x)$ for 
$x\in (x_{i},x_{i}+2A)$. An obvious contradiction occurs at 
$x=x_{i}+2A$ where 
$h_{q}(x_{i}+2A)=0$ because 
$U_{q}(x_{i}+2A)=0$, but the actual film thickness (or rather, the interface height) is 
$R$. Therefore, a more realistic 
${\rm \Delta} x>0$ should be defined. It is reasonable to choose 
${\rm \Delta} x$ to be the meniscus radius, i.e. 
$R_{m}$. However, it varies with time. We propose to use 
${\rm \Delta} x = \langle R_{m} \rangle$, the average value of 
$R_{m}(t)$ (defined using (3.11)) over the first half-period.
For the harmonic meniscus oscillation, one derives from (4.8) and (4.16):
\begin{equation} U_{q}(x)=U_0\left[\frac{ x +{\rm \Delta} x -x_{i}}{ A}\left(2-\frac{ x + {\rm \Delta} x-x_{i}}{A}\right)\right]^{1/2}. \end{equation}
This velocity can now be used to calculate 
$Ca_r$ in (3.12), thus resulting in the quasi-steady thickness 
$h_{q}(x)$. The 
$\tilde h_{q}(\tilde x')$ profiles are shown as solid curves in figure 3. Note that they are independent of 
$\theta _{i}$ (cf. (4.14)), so there is a unique curve in each figure. One can also see the necessity of the 
${\rm \Delta} x$ introduction: if it were not included, the 
$\tilde h_{q}(\tilde x')$ curve would be shifted with respect to 
$\tilde h(\tilde x')$ curves in figure 3. The agreement between the actual and quasi-steady film profiles is good in the central part of the film, which shows that the 
${\rm \Delta} x$ choice is acceptable.
At the beginning of a period, due to the presence of the contact line, 
$R_{m}$ is always larger than 
$R_{m, q}$, which leads to a thinner film: 
$h < h_{q}$, cf. figure 3. This difference becomes relatively less important as 
$P$ increases (compare figures 3a and 3b). To characterise the film at oscillation conditions (e.g. to estimate the film evaporation rate), it is important to know the average film thickness 
$\langle h\rangle$. It is defined as a spatial average over the interval between the contact line and 
$x_{min}$ (abscissa of the point where 
$h_{min}$ is attained),
\begin{equation} \langle h\rangle=\frac{1}{ x_{min}}\int_0^{ x_{min}} h( x, t=0.5 P)\,\mathrm{d}\kern0.06em x. \end{equation} The dependence of 
$\langle h\rangle$ on different parameters can be seen in figure 7. It is an increasing function of 
$\tilde P$ that saturates for 
$\tilde P\to \infty$. It can be compared to the quasi-steady averaged thickness
\begin{equation} \langle h_{q}\rangle=\frac{1}{2 A}\int_{x_{i} - {\rm \Delta} x}^{2 A+x_{i} - {\rm \Delta} x} h_{q}( x)\,\mathrm{d}\kern0.06em x, \end{equation}
which is independent of both 
$A$ and 
$\theta _{i}$ as follows from (4.17). It can be easily calculated without doing any complicated simulations. For small 
$\tilde P$ (i.e. for small 
$\tilde A$), 
$\langle h\rangle <\langle h_{q}\rangle$, mainly because of the contact line vicinity where 
$h(x)< h_{q}(x)$, cf. figure 3. It is not surprising that the saturation value of 
$\langle h\rangle$ increases with 
$\theta _{i}$ (just because of the thicker film near the contact line); however, the increase is weak (figure 7a).
Figure 7. Average thickness of the liquid film at 
$t=0.5P$; 
$\langle h_{q}\rangle$ is shown for comparison. (a) Value of 
$\langle \tilde h\rangle$ as a function of oscillation period for different 
$\theta _{i}$ and for 
$Ca_0=10^{-3}$ and (b) value of 
$\langle h\rangle$ as a function of 
$Ca_0$ for different 
$\tilde P$ and for 
$\theta _{i}=20^\circ$.
In figure 7(b), one can see the dependence of 
$\langle h\rangle$ on 
$Ca_0$. One can see that the quasi-steady average 
$\langle h\rangle _{q}$ (4.19) gives a globally satisfactory approximation of 
$\langle h\rangle$. The increase with 
$Ca_0$ is mainly due to the 
$Ca_0^{2/3}$ factor in (3.7).
4.10. Film shape comparison with experiment
In the experiments of Lips et al. (Reference Lips, Bensalem, Bertin, Ayel, Romestant and Bonjour2010), a capillary tube contains a short liquid plug of pentane in contact with its own vapour at both ends of the tube. One end is connected to a reservoir at constant pressure. The pressure variation at the other end forces the oscillating motion of a liquid plug under isothermal conditions. Such a mode of oscillation leads to the oscillation amplitude increasing in time, which was understood some years later (Signé Mamba et al. Reference Signé Mamba, Magniez, Zoueshtiagh and Baudoin2018). While the amplitude is indeed slightly increasing, the motion is nearly periodic, so the comparison can still be done.
In their experiments, the plug motion is recorded with a high-resolution camera. Both the meniscus velocity and the curvature radius are found from image analysis. The Weber number 
$We_0= 2R\rho U_0^2/\sigma$ is larger than unity (table 2). The Reynolds number 
$Re_0=2R\rho U_0/\mu$ is quite high too and the impact of inertia on the shape of the central meniscus part must be taken into consideration. Indeed, the quasi-steady 
$R_{m}(t)$ evolution and the experimental measurements of Lips et al. (Reference Lips, Bensalem, Bertin, Ayel, Romestant and Bonjour2010) differ (figure 8). One mentions that the measured 
$R_{m}\simeq R$ at 
$t=0, P$, which indicates the complete wetting case. Note the 
$R_{m}$ local minimum around 
$t\simeq 0.8P$. A similar minimum appears in the simulation, see figure 6(c) and the associated discussion.
Figure 8. Liquid plug velocity and left meniscus radius variations during one period of the Lips et al. (Reference Lips, Bensalem, Bertin, Ayel, Romestant and Bonjour2010) experiment (characters) and their polynomial fits (lines) used for the calculation. The experimental data have been made available to us by S. Lips. For comparison, the quasi-steady 
$R_{m,q}$ variation given by (3.11) at receding (
$t\leq 0.5P$) is also shown.
Table 2. Fluid properties at the experimental conditions (1 bar, 
$20\,^\circ$C) of Lips et al. (Reference Lips, Bensalem, Bertin, Ayel, Romestant and Bonjour2010) and key dimensionless numbers.
In spite of high 
$We_0$ and 
$Re_0$ values mentioned above, the thin film can still be considered as controlled by the viscosity only. In the simulation, instead of using the 
$R_{m}$ calculation of § 4.3, the experimental plug velocity and the radius variation shown in figure 8 are used. Under these conditions, the film ripples in the transition region close to the meniscus can be compared to the calculations.
Figure 9 presents several snapshots of plug oscillation. The left column shows the original images of Lips et al. (Reference Lips, Bensalem, Bertin, Ayel, Romestant and Bonjour2010). The liquid film shape in the transition region between the film and the meniscus is enlarged in the middle column. The numerical results are shown in the right column. One can see that the wavy appearance of the interface is truthfully captured by the numerical calculation.
Figure 9. Liquid film shape during meniscus oscillation; two left columns, experimental results by Lips et al. (Reference Lips, Bensalem, Bertin, Ayel, Romestant and Bonjour2010); the right column, numerical results (see the supplementary movie available at https://doi.org/10.1017/jfm.2021.540).
Unfortunately, the quantitative comparison of the film thickness (i.e. the vertical coordinate) is hardly possible since the refraction by the glass capillary is not corrected and the spatial resolution is not sufficient to distinguish the contact line. However, one can compare the axial lengths. The size of one pixel in mm can be obtained from the known outer tube diameter (4 mm) that is visible in the original images. One can compare the axial distance between the local maximum and the local minimum (figure 9) of the film ripple. It can be measured for the images corresponding to 
$t\geq 0.5P$ where the ripple is visible. The distance is almost constant in time. From the experimental images, the axial distance is 
$0.50\pm 0.02$ mm, while from the simulation, it is 
$0.51\pm 0.01$ mm. Evidently, the agreement is excellent.
4.11. Film thickness comparison with experiment
The experiment of Youn et al. (Reference Youn, Han and Shikazono2018) was carried out under adiabatic conditions. They investigated the deposited film thickness of an oscillating meniscus in a cylindrical capillary tube. Two working fluids are selected for the comparison with the present numerical results: water and ethanol. Fluid properties and experimental parameters are summarised in table 3.
Table 3. Fluid properties at the experimental conditions (1 bar, 
$25\,^\circ$C) of Youn et al. (Reference Youn, Han and Shikazono2018). The values of 
$A_0$ and the experimental 
$U$ amplitude are not explicitly given by them and are estimated from their graphs.
In the tests, the capillary tube is partially filled so the syringe piston is in contact with liquid; there is a single meniscus in the tube. The piston is connected to a step motor that imposes the harmonic motion. The other tube end remains open. Initially, the meniscus remains stationary, then, following the piston, starts to oscillate with a constant frequency. Several sensors that measure the film thickness are installed along the tube, within the range of meniscus oscillation. The instantaneous velocity of the meniscus when it passes the sensors is recorded by the high speed camera. Because of the open tube, film evaporation occurs; however, it is much slower than the oscillation velocity so the film thickness variation during the oscillation period is negligibly small. Unfortunately, the film shape near the meniscus and near the contact line was not studied in their experiments and thus cannot be compared with our results.
The numerical calculation has been done as described in § 4 to get the film profiles at 
$t=0.5P$. The experimental film thickness and the numerical results are presented in figure 10 for ethanol and in figure 11 for water. The experimental points are plotted with respect to the meniscus position known from the experimental data.
Figure 10. Liquid film thickness: comparison of the experimental data of Youn et al. (Reference Youn, Han and Shikazono2018), quasi-steady estimation with (3.12) and numerical results for ethanol.
Figure 11. Liquid film thickness: comparison of the experimental data of Youn et al. (Reference Youn, Han and Shikazono2018), quasi-steady estimation with (3.12) and numerical results for water.
The experiment can also be compared with the quasi-steady film thickness calculated with (3.12) (with 
$\tilde h_{s}= 2.5$) where the instantaneous value of the velocity (4.17) is used in 
$Ca_{r}$.
Some experimental parameters are summarised in table 3. The numerical results are in excellent agreement with the quasi-steady data for 2 Hz, where both 
$We_0$ and 
$Re_0$ are moderate. This is not surprising as 
$\tilde P$ is large (see the discussion in § 4.9). Generally, there is a good agreement with the experimental data for the same reason. The discrepancy between the experimental and numerical results is larger for the 6 Hz case where both 
$We_0$ and 
$Re_0$ become large. Both dimensionless numbers are smaller for the ethanol than for the water, and the discrepancy is smaller too. One can conclude that the discrepancy is caused by the inertial effects that become important when 
$Re_0$ attains 500 and 
$We_0$ attains 10.
An additional discrepancy comes from the non-harmonicity of the meniscus velocity profile. As the experimental 
$U(t)$ curves are unavailable, we use the harmonic law (4.8) that corresponds to experimental oscillation period 
$P$ and amplitude 
$A$. Because of the film deposition, the experimental 
$U(t)$ deviates from the harmonic law. This can be observed from table 3. Indeed, the velocity amplitude 
$2{\rm \pi} A/P$ calculated from the oscillation amplitude and the frequency differs substantially from the actual maximum velocity, which points out the non-harmonicity of the oscillations. The comparison would be improved if the experimental 
$U(t)$ were available to us.
5. Combined effect of oscillation and evaporation
In this section we show the implication of the above results for the case of heating conditions (i.e. a positive 
${\rm \Delta} T$). According to (4.6), evaporation occurs for 
${\rm \Delta} T>0$, and the mass flux 
$J(x)\propto h^{-1}(x)$. Instead of solving (4.7) for this case (which is out of the scope of the present article), we apply here the multi-scale reasoning introduced in § 4.1. We consider below a case of a small 
${\rm \Delta} T$ so the evaporation in the film and meniscus regions can be neglected during an oscillation period so the macroscopic results of § 4 still apply. However, because of the above singularity, the effect is not negligible in the microregion (inner region) and leads to a difference between the apparent contact angle 
$\theta$ and the microscopic contact angle 
$\theta _{micro}$, as defined in § 4.1. Their relationship can be expressed as
\begin{equation} \theta = \theta({\rm \Delta} T, \theta_{micro}), \end{equation}cf. Appendix C.2.
When a solution of the full problem exists, the micro- and macroregions can be matched for a given 
${\rm \Delta} T$. They are connected through the formula (5.1), which means that the oscillatory film shape eventually enforces the value of 
$\theta _{micro}$.
In the sense of the multi-scale reasoning, 
$\theta$ must simultaneously satisfy the conditions:
(i) Similarly to the case where
${\rm \Delta} T =0$, the dynamic film shape imposes the $\theta$ value because the contact line is fixed. This implies that the inequality (4.15) should hold for ${\rm \Delta} T\neq 0$ too.(ii) On the other hand,
$\theta$ is defined by (5.1). As shown in Appendix C.2, it is bounded from below: $\theta \geq \theta _{evp\text {-}min} ({\rm \Delta} T)$.
Condition (ii) means that 
$\theta (t)$ remains larger than 
$\theta _{evp\text {-}min}({\rm \Delta} T)$ throughout oscillation. From the inequality (4.15), one thus obtains
\begin{equation} \theta_{min}\geq\theta_{evp\text{-}min}({\rm \Delta} T), \end{equation}
which presents a necessary condition for matching of two regions. With the equality sign 
$\theta _{min}=\theta _{{evp\text {-}min}}({\rm \Delta} T_{max})$, this equation defines a superheating limit 
${\rm \Delta} T_{max}$. Since 
$\theta _{{evp\text {-}min}}({\rm \Delta} T)$ is an increasing function (cf. figure 15), the superheating limit 
${\rm \Delta} T_{max}$ is an upper bound. Thus the inequality (5.2) can hold when 
${\rm \Delta} T < {\rm \Delta} T_{max}$. To determine graphically 
${\rm \Delta} T_{max}$, figure 12 presents an example where 
$\theta (t)$ (bottom and vertical axes, extracted from figure 4a) is plotted together with 
$\theta _{{evp\text {-}min}}({\rm \Delta} T)$ (top and vertical axes, extracted from figure 15). During the oscillation, the minimum value 
$\theta _{min}\simeq 1.8^\circ$ is attained. From the dependence 
$\theta _{evp\text {-}min}({\rm \Delta} T)$, one can deduce that 
${\rm \Delta} T_{max}\simeq 1$ mK. So the solution of the oscillation problem with evaporation is non-existent if 
${\rm \Delta} T>1$ mK.
Figure 12. Apparent contact angle 
$\theta$ calculated for 
$\theta _{micro}=0$ as a function of superheating 
${\rm \Delta} T$ (example of pentane at 1 bar, same as in figure 15) plotted together with the apparent contact angle variation during meniscus oscillation. The latter is a curve for initial value 
$10^\circ$, the same as at the inset of figure 4(a).
The reason for this paradox is the pinned contact line: for the receding contact line, another degree of freedom appears so the contact angle is not constrained any more. The contact line is necessarily depinned when 
$\theta$ attains 
$\theta ({\rm \Delta} T, 0 )$ (actually, a larger value 
$\theta ({\rm \Delta} T, \theta _{rec} )$, but 
$\theta ({\rm \Delta} T, 0 )$ gives a lower bound).
Note that 
$\theta _{min}$ grows slightly with the meniscus velocity amplitude 
$Ca_0$ (figure 5b), so does 
${\rm \Delta} T_{max}$. However, the 
${\rm \Delta} T_{max}$ value remains of the order of mK. As this maximum superheating is considerably smaller than that encountered in practice (where it is rather of several degrees K, see e.g. Fourgeaud et al. Reference Fourgeaud, Nikolayev, Ercolani, Duplat and Gully2017), one can deduce that the contact line receding at evaporation must be accounted for when the meniscus oscillates.
While the calculation has been carried out here only for the pentane case, one can safely state that 
${\rm \Delta} T_{max}$ is much smaller than realistic superheating used in thermal engineering applications for many other fluids.
6. Conclusions
We have analysed the liquid film deposited by an oscillating liquid meniscus in a capillary tube for the case where the contact line is pinned at the farthest meniscus position. The liquid film thickness is not homogeneous because of the varying meniscus velocity. The periodic solution for such a problem has been identified. The average film thickness is one of the most important quantities. It has been analysed depending on the main system parameters, which are the initial contact angle (that at the farthest meniscus position), the period of oscillation and the amplitude of the meniscus velocity represented with the dimensionless capillary number. The average film thickness depends only weakly on the initial contact angle and grows with the oscillation period until saturation. Globally, the average film thickness is well described by the quasi-steady approach.
Both the meniscus curvature and the contact angle vary in time during such a motion. The contact angle remains nearly constant for a large part of the period. This constant value is independent of both the initial contact angle and the oscillation period. The minimal contact angle encountered during oscillation turns out to be an important quantity. It occurs when the largest film ripple approaches the contact line during the meniscus advance. The minimal contact angle weakly depends on both the oscillation period and the initial contact angle and grows with the capillary number. Its value remains, however, small, of the order of several degrees.
Understanding of evaporation that occurs simultaneously with oscillation is important for applications. The strongest impact of evaporation concerns the contact line vicinity where the liquid film is the thinnest. It is shown that the minimal contact angle that occurs during oscillation with the pinned contact line sets an upper bound for the tube superheating, for which a solution for such a problem exists. This upper bound is quite small (e.g. it is 
${\sim }1$ mK for the pentane at 1 bar), much smaller than typical experimental superheating. This shows the necessity of considering the contact line recession during the simultaneous oscillation and evaporation. Such a result is important, e.g. for theoretical modelling of the pulsating heat pipes.
Supplementary movie
Supplementary movie is available at https://doi.org/10.1017/jfm.2021.540.
Acknowledgements
We are indebted to S. Lips for opening the access to his experimental results.
Funding
The present work is supported by the project TOPDESS, financed through the Microgravity Application Program by the European Space Agency. This article is also a part of the PhD thesis of X.Z. co-financed by the CNES and the CEA NUMERICS programme, which has received funding from the European Union Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 800945. An additional financial support of CNES awarded through GdR MFA is acknowledged.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Film relaxation
The meniscus oscillation results in a continuous film shape variation. For this reason, an important quantity is the relaxation time 
$t_{rel}$. It is a characteristic time scale of decrease of a film perturbation caused by the meniscus velocity change. On a time scale 
$\gg t_{rel}$, one expects the meniscus to behave as if the velocity were constant at each time moment (i.e. in a quasi-steady way).
Relaxation of film profile is studied here on an example of a sudden change in the meniscus motion direction, from receding to advancing. The initial profile (red monotonic curve in figure 13a) is that of Landau & Levich (Reference Landau and Levich1942) defined by (3.5) describing the meniscus receding at a constant velocity 
$U_{r}$. In this calculation, 
$R_{m}$ is assumed to remain constant (
$=R$) because 
$Ca_r$ is small. The calculation is performed for 
$Ca_r=10^{-3}$ but the results are independent of 
$Ca_r$ provided it is small enough. The deposited film thickness is given by (3.13).
Figure 13. Film relaxation: profile evolution and relaxation time. (a) Film profile evolution (shown in the frame of reference of the meniscus) from steady recession with the velocity 
$U_{r}$ to steady advance with the velocity 
$U_{a}$ such that 
$|U_{a}|=U_{r}$. The time lag between the two curves is 
${\rm \Delta} \tilde {t} = 0.5$. The inset shows the root mean square deviation from the steady advancing profile and (b) relaxation time as a function of 
$U_{a}/U_{r}$. The dots are the numerical points. The line is the exponential fit.
At 
$t=0$, the meniscus makes a sudden change of its motion direction and for 
$t>0$ advances over the liquid film at a constant velocity 
$U_{a}$. The film is infinite so the meniscus frame of reference and thus (3.2) with 
$U=-U_{a}<0$ are employed. With the scaling of table 1, the dimensionless governing equation is
\begin{equation} \frac{\partial \tilde{h}}{\partial \tilde{t}} + \frac{\partial }{{\partial \tilde x}}\left( {\frac{{\tilde h}^3}{3} \frac{{{\partial ^3}\tilde h}}{{\partial {\tilde x ^3}}}} - \tilde U \tilde h \right) =0, \end{equation}
where 
$\tilde {U}=-U_{a}/U_{r}$. The boundary conditions are (3.3a,b) and (3.6), where both 
$\alpha$ and 
$h_{r}$ are now known.
Equation (A1) is solved numerically, cf. § 4.5 for the details.
The time evolution of the interface profile is shown in figure 13(a). At 
$\tilde t=0$, the profile is given by the Landau & Levich (Reference Landau and Levich1942) profile. As the time increases, the film relaxes to the steady advancing profile (the solid curve in figure 13a), which coincides with the Bretherton (Reference Bretherton1961) rear meniscus profile. By fitting the root mean square deviation (shown in the inset) from the steady advancing profile, one can introduce the relaxation time 
$t_{rel}$; for 
$U_{a}/U_{r}=1$, 
$\tilde t_{rel} \simeq 0.65$.
The calculation shows (figure 13b) that the relaxation time 
$t_{rel}$ decreases with 
$U_{a}$ and remains smaller than 2. As the exponential vanishes after 
$\tilde t\simeq 5\tilde t_{rel}$, one expects that the system behaves independently of a current state after a time lag 
$\tilde t\simeq 10$. This means that, for 
$\tilde {P}>100$ considered above, the transient evolution needs to be considered only in the very beginning and the very end of a period where the contact line affects the overall interface shape. For the remaining part of a period, a quasi-static approach is expected to be valid.
Appendix B. Stokes problem of the straight wedge with a varying angle
We search a solution of the Stokes problem inside a straight two-dimensional liquid wedge (figure 14a) where the opening angle 
$\theta$ varies with the angular velocity 
$\omega$. In polar coordinates (
$r, \varphi$), the Stokes equations for the liquid velocity 
$\boldsymbol v=(v_r, v_\varphi )$ read
$$\begin{gather} \frac{\partial p_l}{\partial r}= \mu\left\{\frac{\partial }{\partial r}\left[\frac{1}{r}\frac{\partial (rv_r)}{\partial r}\right]+\frac{1}{r^2}\frac{\partial^2 v_r}{\partial \varphi^2}-\frac{2}{r^2}\frac{\partial v_\varphi}{\partial \varphi}\right\}, \end{gather}$$
$$\begin{gather}\frac{\partial p_l}{\partial \varphi}= r\mu\left\{\frac{\partial }{\partial r}\left[\frac{1}{r}\frac{\partial (rv_\varphi)}{\partial r}\right]+\frac{1}{r^2}\frac{\partial^2 v_\varphi}{\partial \varphi^2}+\frac{2}{r^2}\frac{\partial v_r}{\partial \varphi}\right\}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol v=0. \end{gather}$$It is well known that this is equivalent to the streamfunction formulation,
\begin{equation} \nabla^4\psi=0, \end{equation}where
$$\begin{gather} v_r =\frac{1}{r}\frac{\partial \psi}{\partial \varphi}, \end{gather}$$
$$\begin{gather}v_\varphi ={-}\frac{\partial \psi}{\partial r} . \end{gather}$$
Equation (B2) admits a solution 
$\psi \propto r^\lambda f(\varphi )$, where 
$\lambda$ is a constant (Moffatt Reference Moffatt1964). On the one hand, 
$\psi$ should scale with 
$\omega$ because it causes the flow. On the other, from (B4), its dimension should be length
$^2$/time. The only choice is thus
\begin{equation} \psi=\omega r^2f(\varphi). \end{equation}
The corresponding function 
$f$ is 
$f(\varphi )=A\cos 2\varphi +B\sin 2\varphi +C\varphi +D$ (Moffatt Reference Moffatt1964), where 
$A,B,C,D$ are constants. They can be determined with the boundary conditions. A zero velocity (
$v_r=v_\varphi =0$) is imposed at the liquid–solid surface 
$\varphi =0$. A radial velocity 
$v_\varphi =\omega r$ is imposed at the liquid–vapour boundary 
$\varphi =\theta$, which is also stress-free (
$\partial v_r / \partial \varphi =0$). The result is
\begin{equation} \psi=\omega r^2\frac{\tan 2\theta(\cos 2\varphi-1)+2\varphi-\sin 2\varphi}{2(\tan 2\theta-2\theta)}. \end{equation}By using it in (B1), one gets
\begin{equation} \left.\begin{array}{c} \displaystyle \dfrac{\partial p_l}{\partial r}=\dfrac{4\mu\omega }{r(\tan 2\theta-2\theta)},\\ \displaystyle \dfrac{\partial p_l}{\partial \varphi}=0, \end{array}\right\} \end{equation}
so 
$p_l\sim \log r$. Note that for a small 
$\theta$, (B7) reduces to
\begin{equation} \frac{\partial p_l}{\partial r}=\frac{3\mu\omega }{2r\theta^3}. \end{equation}
Figure 14. Sketch of the straight liquid wedges with varying angle (a) and curved wedge at evaporation (b).
Appendix C. Lubrication approach to the evaporation problem
C.1. Straight wedge flow
Consider first the case with no phase change. It is described by (2.1), where 
$J=0$. As in the above Stokes problem (Appendix B), one can consider asymptotically (near the contact line) the straight wedge 
$h=\theta x$ with the angular velocity 
$\partial \theta /\partial t=\omega$. The governing equation
\begin{equation} [ \mu\omega x+\frac{\partial}{\partial x}\left( \frac{\theta^3 x^3}{3}\frac{\partial {\rm \Delta} p}{\partial x} \right) =0 \end{equation}results in
\begin{equation} \frac{\partial {\rm \Delta} p}{\partial x}={-}\frac{3\mu\omega}{2x\theta^3}. \end{equation}
With no surprise, this expression agrees with the small 
$\theta$ asymptotics (B8) of the Stokes approach and leads to the logarithmic pressure divergence at the contact line.
Consider now the lubrication theory for the volatile liquids that accounts for the Kelvin effect resulting in (4.7). For the straight wedge with the varying contact angle it becomes
\begin{equation} \mu\omega x + \frac{\partial }{\partial x}\left(\frac{\theta^3 x^3}{3} \frac{\partial{\rm \Delta} p}{\partial x}\right)=\frac{{\rm \Delta} p-{\rm \Delta} p_{cl}}{\theta x}\frac{\mu kT_{sat}}{ ({\mathcal{L}} \rho)^2}. \end{equation}
For the fixed contact angle case (
$\omega =0$), this equation admits an analytical solution (Janeček & Nikolayev Reference Janeček and Nikolayev2012) that satisfies the condition (4.4),
\begin{equation} {\rm \Delta} p={\rm \Delta} p_{cl}\left[1-\frac{\ell_{K}}{x}K_1\left(\frac{\ell_{K}}{x}\right)\right], \end{equation}
where 
$K_1(\cdot )$ is the modified Bessel function of the first order and
\begin{equation} \ell_{K}=\frac{\sqrt{3\mu kT_{sat}}}{{\mathcal{L}} \rho\theta^2} \end{equation}
is a characteristic length of the Kelvin effect. At 
$x\ll \ell _{K}$,
\begin{equation} {\rm \Delta} p={\rm \Delta} p_{cl}\left\{1-\exp\left(-\frac{\ell_{K}}{x}\right) \left[\sqrt{\frac{{\rm \pi} }{2}\frac{\ell_{K}}{x}}+{O}(x^{1/2})\right]\right\}. \end{equation}
For the case of varying contact angle 
$\omega \neq 0$, a solution that satisfies the condition (4.4) can be found as an asymptotic expansion
\begin{equation} {\rm \Delta} p={\rm \Delta} p_{cl}+\frac{3\mu\omega}{\theta^3\ell_{K}^2}x^2+{O}(x^4). \end{equation}C.2. Curved wedge flow caused by evaporation for $\omega =0$
When the substrate is heated, a flow inside the wedge brings the liquid towards the contact line to compensate for the mass loss by evaporation, thus creating the viscous pressure drop described by (C4). It can be seen as a curvature that increases when 
$x\to 0$. The curvature creates a difference between the microscopic contact angle 
$\theta _{micro}$, the actual slope at the contact line, and the interface slope farther away from the contact line, cf. figure 14(b). The characteristic length for this effect is 
$\ell _{K}\sim 10 - 100$ nm (Janeček et al. Reference Janeček, Andreotti, Pražák, Bárta and Nikolayev2013) from the contact line. At a length scale 
$x_{meso}\gg \ell _{K}$ but much smaller than the film-related length scale 
$\sim 10~\mathrm {\mu }$m, one can define the experimentally measurable interface slope 
$\theta$, called the apparent contact angle.
The region 
$x<0< x_{meso}$ is often referred to as the microregion. The numerical calculation of 
$\theta$ is described by Janeček & Nikolayev (Reference Janeček and Nikolayev2012). It is based on the steady version of (4.7) (i.e. with 
$\partial h/\partial t=0$) solved with the boundary conditions (4.9a) and (4.9b). The other two boundary conditions for this fourth-order differential equation are the imposed slope 
$\theta _{micro}$ at the contact line and the condition of zero (on the microregion scale) curvature
\begin{equation} \left.\frac{\partial^2 h}{\partial x^2}\right|_{x= x_{meso}}=0. \end{equation}
Note that the interface slope saturates at 
$x\gg \ell _{K}$ so 
$\theta$ is independent of 
$x_{meso}$.
Figure 15 demonstrates an example of 
$\theta$ as a function of 
$\theta _{micro}$ and 
${\rm \Delta} T$ (cf. (5.1)) for pentane at 1 bar. It turns out that 
$\theta ({\rm \Delta} T, \theta _{micro})$ monotonically grows with both 
$\theta _{micro}$ and intensity of evaporation controlled by 
${\rm \Delta} T$; evidently, 
$\theta ({\rm \Delta} T=0)=\theta _{micro}$.
Figure 15. Value of 
$\theta$ as a function of 
${\rm \Delta} T$ for different 
$\theta _{micro}$ computed for pentane at 1 bar. The curve for 
$\theta _{micro}=0$ corresponds to 
$\theta _{evp\text {-}min}({\rm \Delta} T)$.
One now introduces
\begin{equation} \theta_{evp\text{-}min}({\rm \Delta} T)\equiv \theta({\rm \Delta} T, \theta_{micro}=0), \end{equation}
which is the lower bound of 
$\theta ({\rm \Delta} T, \theta _{micro})$, i.e. its value for the complete wetting case (Janeček et al. Reference Janeček, Andreotti, Pražák, Bárta and Nikolayev2013). Therefore, the 
$\theta _{micro}=0$ curve in figure 15 represents 
$\theta _{evp\text {-}min}({\rm \Delta} T)$. For example, for 
${\rm \Delta} T=$30 mK, one finds 
$\theta _{evp\text {-}min}({\rm \Delta} T=30\ \mathrm {mK})\simeq 8.4^\circ$, see figure 15. This means that, for 
${\rm \Delta} T=$30 mK, it is impossible to have an apparent angle smaller than 
$8.4^\circ$.
