Xiaolong Zhang

Vadim S Nikolayev

Liquid film dynamics with immobile contact line during meniscus oscillation

Keywords: Bubble dynamics, Contact lines, Thin films, Lubrication theory, Evaporation

This paper presents a theoretical analysis of the liquid film dynamics during the oscillation of a meniscus between a liquid and its vapour in a cylindrical capillary. By using the theory of Taylor bubbles, the dynamic profile of the deposited liquid film is calculated within the lubrication approximation accounting for the finiteness of the film length, i.e. for the presence of the contact line. The latter is assumed to be pinned on a surface defect and thus immobile; the contact angle is allowed to vary. The fluid flow effect on the curvature in the central meniscus part is neglected. This curvature varies in time because of the film variation and is determined as a part of the solution. The film dynamics depends on the initial contact angle, which is the maximal contact angle attained during oscillation. The average film thickness is studied as a function of system parameters. The numerical results are compared to existing experimental data and to the results of the quasi-steady approximation. Finally, the problem of an oscillating meniscus is considered accounting for the superheating of the capillary wall with respect to the saturation temperature, which causes evaporation. When the superheating exceeds a quite low threshold, oscillations with a pinned contact line are impossible anymore and the contact line receding caused by evaporation needs to be accounted for.

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. 2013), the distribution of fluids in microfluidics (Angeli & Gavriilidis 2008) 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 et al. 1998) or loop heat pipes, where the pressure oscillations are observed (Launay et al. 2007) and impact the menisci in the capillary structure. The oscillation of menisci is of special importance for the pulsating heat pipes, called also oscillating heat pipes (Fourgeaud et al. 2017;Marengo & Nikolayev 2018;Nikolayev 2021), 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 2012;Savva et al. 2017) so the contact lines are important to model adequately. tube wall

x m h s x s R m parabola (x f ,h f )
Figure 1: Sketch of the axial cross-section of a cylindrical capillary tube with the liquid film deposited by the receding meniscus. At 𝑥 = 𝑥 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 hydrodynamics of menisci has been extensively studied since the seminal articles of Landau & Levich (1942), Taylor (1961) and Bretherton (1961). Since their works, a strong effort has been made to understand the dynamics of the Taylor bubbles (i.e. bubbles of the 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 et al. 2012). 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 sec. 2, the background theory of steady meniscus receding is briefly discussed in sec. 3. The meniscus oscillation is considered in sec. 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 sec. 5.

Model description

Consider a cylindrical capillary tube of an inner radius 𝑅, containing a liquid and its vapour. The vapour-liquid interface is assumed to be axially symmetric (Fig. 1).

The tube is assumed to be thin enough so the gravity force can be neglected. Following the classical approach (Bretherton 1961), 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 Fig. 1 with a circle of radius 𝑅 m ).

Because the vapour has a smaller density and viscosity compared with the liquid, the vapour pressure 𝑝 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 2010) that describes the interface dynamics in the film region, 𝜕ℎ 𝜕𝑡

+ 𝜕 𝜕𝑥 ℎ 3 3𝜇 𝜕Δ𝑝 𝜕𝑥 = - 𝐽 𝜌 , (2.1)
where ℎ = ℎ(𝑥, 𝑡) is the local film thickness, and 𝐽 is the mass flux across the interface, defined to be positive at evaporation. Here, 𝜇 and 𝜌 are the liquid shear viscosity and density, respectively. The lubrication theory is applicable within the assumption |𝜕ℎ/𝜕𝑥| 1. The pressure jump Δ𝑝 = 𝑝 -𝑝 l (with 𝑝 l , the liquid pressure) across the interface obeys the Laplace equation

Δ𝑝 = 𝜎 𝐾 + 1 𝑅 -ℎ 𝜎 𝜕 2 ℎ 𝜕𝑥 2 + 1 𝑅 m , (2.2)
where 𝜎 is the surface tension, 𝐾 is the two-dimensional interface curvature in the axial cross-section shown in Fig. 1; in the small-slope approximation, 𝐾 𝜕 2 ℎ/𝜕𝑥 2 . Because of this limitation, such a film region theory is not able to describe the meniscus region. The radial contribution to the curvature 𝑅 -1 m is assumed to be independent of 𝑥 in the film region (where ℎ is much smaller than 𝑅). At a large 𝑥, Eq. (2.1) results in an increasingly larger ℎ, where the viscous forces vanish and the surface tension alone controls the interface, so its curvature 𝜕 2 ℎ/𝜕𝑥 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

𝜕 2 ℎ/𝜕𝑥 2 | 𝑥=𝑥 f = 𝑅 -1 m .
(2.3) defined at the ending point (𝑥 f , ℎ f ) of the film region.

Isothermal problem with infinite film: steady solutions

First, we would like to recall the theory for the case with no phase change (𝐽 = 0), so Eq. (2.1) becomes 𝜕ℎ 𝜕𝑡

+ 𝜎 3𝜇 𝜕 𝜕𝑥 ℎ 3 𝜕 3 ℎ 𝜕𝑥 3 = 0.
(3.1)

The velocity of the axial meniscus centre (assumed positive for a receding meniscus according to the 𝑥-axis direction choice) is denoted 𝑈. 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 𝑥 = 𝑥 -𝑈𝑡. Eq. (3.1) can then be rewritten as

𝜕ℎ 𝜕𝑡 + 𝜕 𝜕𝑥 𝜎 3𝜇 ℎ 3 𝜕 3 ℎ 𝜕𝑥 3 -𝑈ℎ = 0. (3.2)
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 𝑥 will be as in (2.1) and (3.1). The steady version of Eq. (3.2) for the case of a constant positive velocity 𝑈 r (that should be used for 𝑈 in this case) is the Landau-Levich equation (Landau & Levich 1942) describing the flat infinite film being deposited by a receding meniscus. Note that the film in a cylindrical capillary (Bretherton 1961) is described by the same equation because of the approximation (2.2).

The boundary conditions at 𝑥 → -∞ describe the flat film of the yet unknown thickness ℎ r , ℎ = ℎ r , 𝜕ℎ/𝜕𝑥 = 0.

(3.

3)

The scaling of this problem (Table 1) is based on ℎ r . The characteristic axial length scale involves the capillary number 𝐶𝑎 r = 𝜇𝑈 r /𝜎 and is chosen in such a way that the dimensionless equation

ℎ r 𝐶𝑎 -1/3 r 𝛼𝑅𝐶𝑎 1/3 0 film thickness ℎ ℎ r 𝛼𝑅𝐶𝑎 2/3 0 time 𝑡 ℎ r 𝐶𝑎 -1/3 r 𝑈 -1 r 𝛼𝑅𝐶𝑎 1/3 0 𝑈 -1 0 velocity 𝑈 𝑈 r 𝑈 0 = 2𝜋 𝐴/𝑃
𝜕 𝜕 x h3 3 𝜕 3 h 𝜕 x3 -h = 0, (3.4)
does not contain any constants; the tilde means hereafter the corresponding dimensionless variable. By integrating Eq. (3.4) from -∞ to x and using the conditions (3.3) (in the dimensionless form, h( x → -∞) = 1), one obtains

𝜕 3 h 𝜕 x3 = 3 h -1 h3 , (3.5)
which is equivalent to the Bretherton 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 x. Eq. (3.5) remains valid until the transition region (between the film and the meniscus), in which 1 h < 𝑅. From Eq. (3.5), 𝜕 3 h/𝜕 x3 0. This means that, at large x, the second derivative is finite, which is compatible to the condition (2.3). Eq. (3.5) can be integrated numerically (see Nikolayev & Sundararaj (2014) for details). The ending point ( xf , hf ) of the integration interval is chosen from a condition that 𝜕 2 h/𝜕 x2 reaches a plateau (with a required accuracy). Such a calculation results in a numerical value for this plateau

𝜕 2 h 𝜕 x2 x= xf = 𝛼 1.3375 (3.6)
The resulting profile h( x) can be found in Fig. 13a of Appendix A. Note that 𝛼 is equivalent to the numerical value 0.643 originally found by Bretherton (1961); 𝛼 0.643 • 3 2/3 , where the factor appears because of the different scaling. By comparing Eqs. (2.3, 3.6), one finds the expression for the film thickness

ℎ r = 𝛼𝑅 m 𝐶𝑎 2/3 r . (3.7)
The only yet unknown quantity is 𝑅 m that can be found as proposed by Klaseboer et al. (2014). As mentioned before, near the point ( xf , hf ) the film shape should satisfy the condition (2.3) which means that at ℎ ℎ r it asymptotically approaches a parabola

𝑦 = (𝑥 -𝑥 s ) 2 /(2𝑅 m ) + ℎ s , (3.8)
where the parameters (𝑥 s , ℎ s ) are yet to be determined. To understand their meaning, one recalls that near its minimum where its curvature is 𝑅 -1 m , the parabola approximates the circular meniscus profile

(𝑥 -𝑥 s ) 2 + (𝑦 -ℎ s -𝑅 m ) 2 = 𝑅 2 m .
(3.9)

It is evident now that (𝑥 s , ℎ s ) is the circle lowest point. One can obtain (𝑥 s , ℎ s ) by fitting the film profile h( x) near the point ( xf , hf ) to a parabola. Klaseboer et al. (2014) report that the dimensionless value of hs slightly grows with hf . Bretherton (1961) gives hs = 2.79. The asymptotic value hs = 2.90 is obtained for hf > 10 6 . However, to obtain a continuous overall interface profile, the matching point ( xf , hf ) film-parabola should be lower than the point where the parabola-circle transition occurs. This requires ℎ f < 𝑅. For the value hf 50 that satisfies this condition in practical situations, Klaseboer et al. (2014) find hs = 2.5, which is also the value found from the experimental data fits, as discussed below. In summary, the hs variation is weak and one can consider that the circle is nearly invariant of the specific hf choice. Klaseboer et al. (2014) have proposed the equation

𝑅 m + ℎ s = 𝑅 (3.10)
that centres the circle with respect to the tube and thus links 𝑅 m to 𝑅, cf. (3.12)

The value hs 2.5 has been determined by Aussillous & Quéré (2000) from the experimental data fits. For 𝐶𝑎 r → 0, 𝑅 m 𝑅 and Eq. (3.12) reduces to Bretherton's original expression ℎ r = 𝛼𝑅𝐶𝑎 2/3 r .

(3.13) One can consider the meniscus advancing at a constant velocity 𝑈 = -𝑈 a (where 𝑈 a is the modulus of the advancing velocity) over the pre-existing film of thickness ℎ r . Such a motion has been understood as well. It has been shown (Bretherton 1961) that the film has a wavy shape (ripples) near the meniscus, cf. the solid curve in Fig. 13a of Appendix A. The wavelength of ripples depends on the ratio 𝑈 a /𝑈 r (Maleki et al. 2011;Nikolayev & Sundararaj 2014), where 𝑈 r can be deduced (with Eq. 3.13) from ℎ r . The meniscus radius for the steady advancing case was determined with Eq. (3.10) by Cherukumudi et al. (2015).

Oscillations in the presence of a contact line

A previous study (Nikolayev & Sundararaj 2014) 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 Fig. 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. 2016;Rao et al. 2017), 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 2015). 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 et al. 2014) is that the microscopic contact angle averaged along the contact line can vary between the static advancing 𝜃 adv and the static receding 𝜃 rec angles while the contact line remains immobile. In our calculation, 𝜃 adv -𝜃 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 • (de Gennes 1985), 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 1987) and the flow is known to be well described by the lubrication approximation, even down to the nanometric scale (Mortagne et al. 2017).

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. Eq. B 8) 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 to be 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 𝑇 sat corresponding to the imposed vapour pressure 𝑝. The wall superheating is denoted Δ𝑇. The tube wall temperature is thus 𝑇 w = 𝑇 sat + Δ𝑇.

Conventional hypotheses concerning the liquid film mass exchange (Nikolayev 2010) 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

𝐽 = 𝑘 (𝑇 w -𝑇 int ) ℎL , (4.1)
where 𝑇 int is the temperature of the vapour-liquid interface, 𝑘 is the liquid heat conductivity and L is the latent heat. The film is assumed here to be thin with respect to 𝑅 so the onedimensional 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 Eq. (2.1). Generally, the film thinning is not strong during an oscillation period (Fourgeaud et al. 2017).

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 (𝑇 int = 𝑇 sat ), the mass flux 𝐽 (4.1) would diverge at the contact line ℎ = 0 as 𝐽 ∼ Δ𝑇/ℎ, which is non-physical because total evaporated mass (the integral of 𝐽) would be infinite.

The Kelvin effect, i.e. the dependence of 𝑇 int on the interfacial pressure jump Δ𝑝

𝑇 int = 𝑇 sat 1 + Δ𝑝 L 𝜌 (4.2)
can relax the singularity (Janeček & Nikolayev 2012), because it allows 𝑇 int to vary along the interface so it can attain the wall temperature 𝑇 w at the contact line so the mass flux

𝐽 (𝑥 → 0) = 0. (4.3)
From the temperature continuity, one obtains the condition

Δ𝑝(𝑥 → 0) = Δ𝑝 cl , (4.4)
where a constant pressure jump at the contact line is introduced as

Δ𝑝 cl = L 𝜌 𝑇 sat Δ𝑇 . (4.5)
One can show that a solution that satisfies this condition can indeed be found (cf. Appendix C).

Eqs. (4.1, 4.2) result in

𝐽 = 𝑘 ℎL Δ𝑇 -Δ𝑝 𝑇 sat L 𝜌 . (4.6)
With its substitution into Eq. (2.1), the governing equation becomes 𝜕ℎ 𝜕𝑡

+ 𝜕 𝜕𝑥 ℎ 3 3𝜇 𝜕Δ𝑝 𝜕𝑥 = Δ𝑝 -Δ𝑝 cl ℎ 𝑘𝑇 sat (L 𝜌) 2 . (4.7)
The problem is now regular (because Δ𝑝 is not divergent anymore), unlike other microscopic approaches (Savva et al. 2017). 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 2012) 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 ℓ K (C 4), cf. Appendix C for more details. It is nanometric (Janeček et al. 2013) and is thus significantly smaller than the characteristic scale of film shape variation that we call macroscopic. For this reason, Eq. (4.7) can be understood within a multi-scale paradigm in the spirit of the asymptotic matching techniques (Janeček et al. 2013). In the inner region, commonly called the microregion, the first (transient) term is negligible with respect to the Kelvin term (Δ𝑝 containing term in the r.h.s.). The problem is reduced to that of Appendix C.2. To summarise it, when Δ𝑇 ≠ 0, a strong interfacial curvature that exists in the microregion can cause a difference between the microscopic contact angle 𝜃 micro and the interface slope 𝜃 defined at 𝑥 → ∞ within the microregion. In the outer (macroscopic) region, the Kelvin effect is negligible so the Δ𝑝 term on the r.h.s. of Eq. (4.7) vanishes. For Δ𝑇 = 0, this equation is that of the isothermal problem (3.1). Both problems can be matched at scale 𝑥 meso ℓ K , much smaller than the macroscopic scale. The apparent contact angle visible on this latter scale is thus equal to 𝜃. It is assumed hereafter that the pinning occurs at a length scale smaller than ℓ K , i.e. there are nanometric defects with sharp borders on which the contact line is pinned.

In sec. 4, a globally isothermal problem is considered, 𝑇 w = 𝑇 sat , so Δ𝑇 = 0 and Δ𝑝 cl = 0, which means 𝜃 = 𝜃 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 Δ𝑝 can be large, as mentioned above. The mass flux 𝐽 scales with Δ𝑝 according to Eq. (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. 2013). Consider e.g. an increasing in time 𝜃. According to Eq. (C 6), Δ𝑝 > 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 ℓ K .

In conclusion, all results obtained in the present sec. 4, can be seen as obtained with Eq. (3.1) because the microregion details cannot be resolved at the macroscale pictured in the figures below. However, the numerical calculations of the regularised Eq. (4.7) are carried out in reality.

Oscillation problem statement

The meniscus now oscillates, and the position 𝑥 m of its centre (which is the experimentally measurable quantity) travels periodically with a period 𝑃 and an amplitude 𝐴. One can assume its harmonic oscillation

𝑥 m (𝑡) = 𝑥 i + 𝐴[1 -cos(2𝜋𝑡/𝑃)], (4.8)
where 𝑥 i is the initial meniscus centre position. Alternatively, one can take the experimentally measured dependence 𝑥 m (𝑡) while comparing the data with the experiment (cf. sec. 4.10 below). The contact line is pinned at the position 𝑥 = 0, and the contact angle 𝜃 varies.

For the harmonic oscillation case, the meniscus velocity is 𝑈 (𝑡) = 𝑈 0 sin(2𝜋𝑡/𝑃), where the velocity amplitude 𝑈 0 = 2𝜋 𝐴/𝑃 is convenient to choose as the characteristic velocity to define the capillary number 𝐶𝑎 0 = 𝜇𝑈 0 /𝜎 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. Eq. (4.7) (with the substitution of Eq. 2.2) is solved for 𝑥 ∈ [0, 𝑥 f ]. The length 𝑥 f is imposed as explained in secs. 4.3, 4.5 below.

The boundary conditions are defined as

ℎ (𝑥 = 0) = 0, (4.9a) 𝜕Δ𝑝 𝜕𝑥 𝑥→0 = 0, (4.9b) Δ𝑝 (𝑥 = 𝑥 f ) = 𝑅 -1 m , (4.9c) ℎ (𝑥 = 𝑥 f ) = ℎ f , (4.9d)
where 𝑅 m and ℎ f are discussed in sec. 4.3. The condition (4.9a) is a geometrical constraint at the contact line. Eq. (4.9b) is a weaker form of the condition (4.4) used to provide numerical stability. The boundary conditions (4.9c, 4.9d) impose the liquid film curvature and thickness at the right end of the integration interval for each time moment.

Determination of the meniscus curvature

For a small film thickness, one can assume that the meniscus radius 𝑅 m is constant and equal to 𝑅 during oscillation. It is actually a good approximation for a small 𝐶𝑎 0 10 -3 . However, at a larger 𝐶𝑎 0 , the film thickness impacts 𝑅 m (see sec. 4.8 below). Since the film thickness depends on the meniscus velocity, so does the meniscus radius 𝑅 m , cf. sec. 3. Therefore, 𝑅 m varies in time. In this section, we generalise to any meniscus dynamics the method for 𝑅 m determination (Klaseboer et al. 2014) discussed above for the steady receding case.

Similarly to the steady case of sec. 3, one needs first to match the film shape ℎ(𝑥) 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)

ℎ f = (𝑥 f -𝑥 s ) 2 /(2𝑅 m ) + ℎ s , (4.10a) 𝜕ℎ 𝜕𝑥 𝑥=𝑥 f = 𝑥 f -𝑥 s 𝑅 m , (4.10b)
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. (2014), Eq. (3.10) serves to find 𝑅 m .

We introduce in addition a relationship of the abscissas of the lowest and rightmost points of a circle that is needed to define 𝑥 s (Fig. 1):

𝑅 m + 𝑥 s = 𝑥 m . (4.10c)
In the present algorithm, 𝑥 f imposed to such a value that the difference 𝑥 m -𝑥 f does not vary in time and ℎ f = ℎ(𝑥 f ) remains large with respect to the deposited film thickness. As discussed in sec. 3, the solution is nearly independent of the specific choice of ℎ f (and thus of 𝑥 f ). The set of Eqs. (2.2), (3.10) and (4.7)-(4.10) is then complete, so the film shape and the unknown parameters (𝑅 m , 𝑥 s , ℎ s , ℎ f ) can be determined for each 𝑡.

Initial conditions and solution periodicity

One needs to define now the initial film shape ℎ(𝑥, 0) at 𝑡 = 0, which corresponds to the (yet unspecified) leftmost meniscus position 𝑥 i according to Eq. (4.8). As an initial film profile, we choose that of equilibrium satisfying the condition 𝜕ℎ/𝜕𝑡 = 0 that can be used in Eq. (3.1). From the boundary condition (4.9b), one finds straightforwardly 𝜕 3 ℎ/𝜕𝑥 3 = 0, i.e. the parabolic shape

ℎ(𝑥, 0) = 𝑥 2 2𝑅 m,i + 𝜃 i 𝑥, (4.11a)
where 𝜃 i ≡ 𝜃 (𝑡 = 0) is the initial contact angle. It serves as another boundary condition, additional to (4.9a) and (4.9c). By applying Eqs. (3.10) and (4.10) at 𝑡 = 0, one gets

𝑅 m,i ≡ 𝑅 m (𝑡 = 0) = 𝑅/(1 -𝜃 2 i /2), (4.11b) 𝑥 i = 𝑅 m,i (1 -𝜃 i ).
(4.11c)

These expressions are the small-angle approximations of the expressions 𝑅 m,i = 𝑅/cos 𝜃 i and

𝑥 i = 𝑅 m,i (1 -sin 𝜃 i ) because Eq. (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 𝜕ℎ/𝜕𝑡 is always balanced by the curvature gradient, more precisely, by the second term of Eq. (3.1). This occurs because there are no other forces, in particular, no inertia. When the meniscus comes near the leftmost position, 𝑈 decreases and the system approaches the state (4.11) with no curvature gradient, thus, |𝜕ℎ/𝜕𝑡| 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 𝑡 rel 𝑃, which is our case (cf. Appendix A for 𝑡 rel discussion). Indeed, the numerical simulations show that ℎ(𝑥, 𝑃) is indistinguishable from ℎ(𝑥, 0), cf. Fig. 2 below. So do all other parameters (curvature, contact angle, etc.). This finding allows us to simulate a unique period.

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: 𝜃 i , P and 𝐶𝑎 0 . All quantities will be studied in this parametric space. The dimensionless amplitude is linked to the period à = P/(2𝜋), (4.12) where a dimensionless quantity is denoted with a tilde. There is one more dimensionless parameter

𝑁 = 𝜇𝑘𝑇 sat L 𝜌𝛼𝑅 𝐶𝑎 0 (4.13)
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 (C 4) is chosen to be small enough (cf. sec. 4.1). The mesh size is exponentially refined near the contact line (as x → 0) to capture the contact angle variation without considerably increasing the total number of nodes (Nikolayev 2010). At t = 0, a value of xf = 10, for which hf is around 50, cf. the discussion in sec. 3. At t > 0, the difference xmxf is maintained constant, and equal to that at t = 0. To avoid the discretisation error for the contact angle, the initial interface profile is determined numerically by solving the equilibrium version of Eq. (4.7), instead of using the analytical profile (4.11a).

Eq. (4.7) is solved numerically with the finite volume method (FVM), which is more stable numerically (Patankar 1980) than a more conventional finite difference method. In one dimension, a finite volume is just a segment. The variables such as ℎ and their evenorder 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 (2010).

One is interested in amplitudes that are large with respect to the meniscus width 𝑥 m -𝑥 s , which means large dimensionless periods of oscillation, see Eq. (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 PC.

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 𝑡 = 𝑃/2. For 𝑡 > 𝑃/2, the meniscus advances over the deposited film. The ripples near the meniscus appear, like during the steady meniscus advance discussed in sec. 3. The interface profiles h( x, 0) and h( x, P) are indistinguishable, which confirms the periodicity of the oscillations.

Fig. 3 shows examples of interface profiles at 𝑡 = 𝑃/2 (when the meniscus is at its rightmost position so the film length attains its maximum) for several values of the initial contact angle 𝜃 i . All the film profiles are presented in the meniscus reference

𝑥 = 𝑥 -2𝐴 -𝑥 i , (4.14)
meaning that the meniscus centre is at 𝑥 = 0, cf. Eq. (4.8). One can see that the interface shape near the meniscus is independent of 𝜃 i . During the meniscus receding, 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 ℎ min discussed in sec. 4.7 below. It appears because of the meniscus velocity reduction at the end of a half-period.

Contact angle during oscillation

A typical variation of 𝜃 during oscillation is plotted in Fig. 4a. The initial contact angle 𝜃 i is the maximum contact angle achieved during the periodic motion.

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 𝜃 i and 𝑃. 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 Fig. 3b in the contact line region for the curve corresponding to 𝜃 i = 40 • .

The ridge width slowly grows so 𝜃 slowly decreases until the ripples in the near-meniscus region approach the contact line during the backward stroke (Fig. 4b) at the end of a period. This causes the contact angle oscillations, during which its minimal value 𝜃 min is attained 𝜃 (𝑡) 𝜃 min .

(4.15) It depends quite weakly on 𝜃 i (Fig. 4a). The variation of 𝜃 min with the system parameters follows the variation of ℎ min (which is a minimum of ℎ(𝑥, 0.5𝑃) observed near the meniscus). This is illustrated in Fig. 5, where the variations of ℎ min and 𝜃 min with the system parameters are compared. Only the variations with P and 𝐶𝑎 0 are considered (as mentioned above, the dependence on 𝜃 i is quite weak). Evidently, the variations of ℎ min and 𝜃 min with the system parameters are similar, which shows their intrinsic link. At 𝐶𝑎 0 → 0, the dependence on 𝐶𝑎 0 is weak, but becomes stronger at large 𝐶𝑎 0 . This minimal value of the contact angle is of importance (cf. sec. 5 below). Since the motion is periodic, the contact angle 𝜃 i is attained at 𝑡 = 𝑃.

Meniscus curvature during oscillation

The meniscus curvature (Figs. 6) changes periodically during oscillations. The value of 𝑅 m can be compared to the quasi-steady value 𝑅 m,q given by Eq. (3.11) where hs = 2.5 and 𝐶𝑎 𝑟 are calculated with the instantaneous meniscus velocity 𝑈 = 𝑈 0 sin(2𝜋t/ P) during the liquid receding (𝑡 0.5𝑃).

For 𝑡 = 0, 𝑅 m is defined with Eq. (4.11b), which differs from the quasi-steady value 𝑅 m,q = 𝑅 for 𝑈 = 0. This difference occurs because of the contact line presence. In its absence (pre-wetted tube, the situation equivalent to the limit 𝜃 i → 0 in our model), the initial radius would be close to 𝑅 because the wetting film is much thinner than the film considered here. Within the time scale ∼ 5𝑡 rel (see Appendix A), 𝑅 m relaxes to the quasi-steady value 𝑅 m,q . The curvature 𝑅 m remains close to 𝑅 m,q until the deceleration that occurs near the rightmost meniscus position (at 𝑡 = 0.5𝑃, where 𝑈 = 0). However, the shape relaxation causes a delay, so the inequality 𝑅 m < 𝑅 m,q (𝑈 = 0) = 𝑅 always holds at the point 𝑡 = 0.5𝑃. During the backstroke, 𝑅 m varies, finally attaining the initial value (4.11b) that depends only on 𝜃 i (Fig. 6b). Evidently, the amplitude of 𝑅 m oscillation grows with 𝐶𝑎 0 (Fig. 6c), as foreseen by Eq. (3.11). One also mentions the non-monotonic 𝑅 m variation during the backstroke with a local minimum around 𝑡 0.8 -0.9𝑃 (Fig. 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 (Figs. 3), which correlates with the maximum of the meniscus velocity. It is thus interesting to compare the film thickness with its quasisteady value. Within the quasi-steady approach, the term 𝜕ℎ/𝜕𝑡 is neglected and the quasisteady thickness ℎ q (𝑥) can be defined as corresponding to the meniscus receding velocity as if it was constant at each time moment (Fourgeaud et al. 2017;Youn et al. 2018). In sec. 4.8 it has been shown that a simple quasi-steady approach predicts well the meniscus curvature for 0 < 𝑡 < 0.5𝑃. It is interesting to see if it is efficient in predicting the film profile. In this section, only the profile at 𝑡 = 0.5𝑃 is considered, i.e. that with the longest film.

A difficulty appears because the film thickness ℎ(𝑥, 0.5𝑃) depends on 𝑥. It is clear that ℎ depends on the velocity that the meniscus had at film deposition, but at which time moment?

The most obvious first option is a moment 𝑡 when the meniscus centre was at the point 𝑥 (i.e. 𝑥 m (𝑡) = 𝑥). A more sophisticated option is a moment 𝑡 > 𝑡 such that

𝑥 m (𝑡 ) = 𝑥 + Δ𝑥 (4.16)
with Δ𝑥 > 0. In the previous approaches (Fourgeaud et al. 2017;Youn et al. 2018), the first option was used. By using Eq. (4.8) one easily finds that this assumption defines 𝑈 q (𝑥) for 𝑥 ∈ (𝑥 i , 𝑥 i + 2𝐴). An obvious contradiction occurs at 𝑥 = 𝑥 i + 2𝐴 where ℎ q (𝑥 i + 2𝐴) = 0 because 𝑈 q (𝑥 i + 2𝐴) = 0, but the actual film thickness (or rather, the interface height) is 𝑅. Therefore, a more realistic Δ𝑥 > 0 should be defined. It is reasonable to choose Δ𝑥 to be the meniscus radius, i.e. 𝑅 m . However, it varies with time. We propose to use Δ𝑥 = 𝑅 m , the average value of 𝑅 m (𝑡) (defined using Eq. 3.11) over the first half-period.

For the harmonic meniscus oscillation, one derives from Eqs. (4.8,4.16):

𝑈 q (𝑥) = 𝑈 0 𝑥 + Δ𝑥 -𝑥 i 𝐴 2 - 𝑥 + Δ𝑥 -𝑥 i 𝐴 1/2 . (4.17)
This velocity can now be used to calculate 𝐶𝑎 𝑟 in Eq. (3.12), thus resulting in the quasisteady thickness ℎ q (𝑥). The hq ( x ) profiles are shown as solid curves in Figs. 3. Note that they are independent of 𝜃 i (cf. Eq. 4.14), so there is a unique curve in each figure. One can also see the necessity of the Δ𝑥 introduction: if it were not included, the hq ( x ) curve would be shifted with respect to h( x ) curves in Figs. 3. The agreement between the actual and quasi-steady film profiles is good in the central part of the film, which shows that the Δ𝑥 choice is acceptable.

At the beginning of a period, due to the presence of the contact line, 𝑅 m is always larger than 𝑅 m,q , which leads to a thinner film: ℎ < ℎ q , cf. Figs. 3. This difference becomes relatively less important as 𝑃 increases (compare Figs. 3a and3b).

To characterise the film at oscillation conditions (e.g. to estimate the film evaporation rate), it is important to know the average film thickness ℎ . It is defined as a spatial average over the interval between the contact line and 𝑥 min (abscissa of the point where ℎ min is attained),

ℎ = 1 𝑥 min ∫ 𝑥 min 0 ℎ(𝑥, 𝑡 = 0.5𝑃)d𝑥. (4.18)
The dependence of ℎ on different parameters can be seen in Figs. 7. It is an increasing function of P that saturates for P → ∞. It can be compared to the quasi-steady averaged thickness

ℎ q = 1 2𝐴 ∫ 2𝐴+𝑥 i -Δ𝑥 𝑥 i -Δ𝑥 ℎ q (𝑥)d𝑥 (4.19)
which is independent of both 𝐴 and 𝜃 i as follows from Eq. (4.17). It can be easily calculated without doing any complicated simulations. For small P (i.e. for small Ã), ℎ < ℎ q , mainly because of the contact line vicinity where ℎ(𝑥) < ℎ q (𝑥), cf. Figs. 3. It is not surprising that the saturation value of ℎ increases with 𝜃 i (just because of the thicker film near the contact line); however, the increase is weak (Fig. 7a). In Fig. 7b, one can see the dependence of ℎ on 𝐶𝑎 0 . One can see that the quasi-steady average ℎ q Eq. ( 4.19) gives a globally satisfactory approximation of ℎ . The increase with 𝐶𝑎 0 is mainly due to the 𝐶𝑎 2/3 0 factor in Eq. (3.7).

Film shape comparison with experiment

In the experiments of Lips et al. (2010), 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. 2018). 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 𝑊 𝑒 0 = 2𝑅𝜌𝑈 2 0 /𝜎 is larger than unity (table 2). The Reynolds number Re 0 = 2𝑅𝜌𝑈 0 /𝜇 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 𝑅 m (𝑡) evolution and the experimental measurements of Lips et al. (2010) differ (Fig. 8). One mentions that the measured 𝑅 m 𝑅 at 𝑡 = 0, 𝑃, which indicates the complete wetting case. Note the 𝑅 m local minimum around 𝑡 0.8𝑃. A similar minimum appears in the simulation, see Fig. 6c and the associated discussion.

In spite of high 𝑊 𝑒 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 𝑅 m calculation of sec. 4.3, the experimental plug velocity and the radius variation shown in Fig. 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. (2010). 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.

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 (Fig. 9) of the film ripple. It can be measured for the images corresponding to 𝑡 0.5𝑃 where the ripple is visible. The distance is almost constant in time. From the experimental images, the axial distance is 0.50 ± 0.02 mm, while from the simulation, it is 0.51 ± 0.01 mm. Evidently, the agreement is excellent. (2018), quasi-steady estimation with Eq. (3.12) and numerical results for ethanol.

U/U 0 R m /R t/P U experimental R m,q R m experimental

Film thickness comparison with experiment

The experiment of Youn et al. (2018) 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.

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 sec. 4 to get the film profiles at 𝑡 = 0.5𝑃. The experimental film thickness and the numerical results are presented in Fig. 10 for ethanol and in Fig. 11 for water. The experimental points are plotted with respect to the meniscus position known from the experimental data.

The experiment can also be compared with the quasi-steady film thickness calculated with Eq. (3.12) (with hs 2.5) where the instantaneous value of the velocity (4.17) is used in 𝐶𝑎 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 𝑊 𝑒 0 and Re 0 are moderate. This is not surprising as P is large (see the discussion in sec. 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 𝑊 𝑒 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 𝑊 𝑒 0 attains 10.

An additional discrepancy comes from the non-harmonicity of the meniscus velocity profile. As the experimental 𝑈 (𝑡) curves are unavailable, we use the harmonic law (4.8) that corresponds to experimental oscillation period 𝑃 and amplitude 𝐴. Because of the film deposition, the experimental 𝑈 (𝑡) deviates from the harmonic law. This can be observed from table 3. Indeed, the velocity amplitude 2𝜋 𝐴/𝑃 calculated from the oscillation amplitude and (2018). 𝐴 0 and the experimental 𝑈 amplitude are not explicitly given by them and are estimated from their graphs.

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 𝑈 (𝑡) were available to us.

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 Δ𝑇). According to Eq. (4.6), evaporation occurs for Δ𝑇 > 0, and the mass flux 𝐽 (𝑥) ∝ ℎ -1 (𝑥). Instead of solving Eq. (4.7) for this case (which is out of the scope of the present article), we apply here the multi-scale reasoning introduced in sec. 4.1. We consider below a case of a small Δ𝑇 so the evaporation in the film and meniscus regions can be neglected during an oscillation period so the macroscopic results of sec. 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 𝜃 and the microscopic contact angle 𝜃 micro , as defined in sec. 4.1. Their relationship can be expressed as

𝜃 = 𝜃 (Δ𝑇, 𝜃 micro ), (5.1) cf. Appendix C.2.
When a solution of the full problem exists, the micro-and macroregions can be matched for a given Δ𝑇. They are connected through the formula (5.1), which means that the oscillatory film shape eventually enforces the value of 𝜃 micro .

In the sense of the multi-scale reasoning, 𝜃 must simultaneously satisfy the conditions: (i) Similarly to the case where Δ𝑇 = 0, the dynamic film shape imposes the 𝜃 value because the contact line is fixed. This implies that the inequality (4.15) should hold for Δ𝑇 ≠ 0 too.

(ii) On the other hand, 𝜃 is defined by Eq. (5.1). As shown in Appendix C.2, it is bounded from below: 𝜃 𝜃 evp-min (Δ𝑇).

Condition (ii) means that 𝜃 (𝑡) remains larger than 𝜃 evp-min (Δ𝑇) throughout oscillation. From the inequality (4.15), one thus obtains 𝜃 min 𝜃 evp-min (Δ𝑇),

(5.2) which presents a necessary condition for matching of two regions. With the equality sign 𝜃 min = 𝜃 evp-min (Δ𝑇 max ), this equation defines a superheating limit Δ𝑇 max . Since 𝜃 evp-min (Δ𝑇) is an increasing function (cf. Fig. 15), the superheating limit Δ𝑇 max is an upper bound. Thus the inequality (5.2) can hold when Δ𝑇 < Δ𝑇 max . To determine graphically Δ𝑇 max , Fig. 12 presents an example where 𝜃 (𝑡) (bottom and vertical axes, extracted from Fig. 4a) is plotted together with 𝜃 evp-min (Δ𝑇) (top and vertical axes, extracted from Fig. 15). During the oscillation, the minimum value 𝜃 min 1.8 • is attained. From the dependence 𝜃 evp-min (Δ𝑇), one can deduce that Δ𝑇 max 1 mK. So the solution of the oscillation problem with evaporation is non-existent if Δ𝑇 > 1mK.

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 𝜃 attains 𝜃 (Δ𝑇, 0) (actually, a larger value 𝜃 (Δ𝑇, 𝜃 rec ) but 𝜃 (Δ𝑇, 0) gives a lower bound).

Note that 𝜃 min grows slightly with the meniscus velocity amplitude 𝐶𝑎 0 (Fig. 5b), so does Δ𝑇 max . However, the Δ𝑇 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. (2017)), 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 Δ𝑇 max is much smaller than realistic superheating used in thermal engineering applications for many other fluids.

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 Figure 12: Apparent contact angle 𝜃 calculated for 𝜃 micro = 0 as a function of superheating Δ𝑇 (example of pentane at 1bar, same as in Fig. 15) plotted together with the apparent contact angle variation during meniscus oscillation. The latter is a curve for initial value 10 • , the same as at the inset of Fig. 4a.

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 ∼ 1 mK for the pentane at 1 bar), much smaller than typical experimental superheating. This shows the necessity of considering the contact line receding during the simultaneous oscillation and evaporation. Such a result is important, e.g. for theoretical modelling of the pulsating heat pipes. a film perturbation caused by the meniscus velocity change. On a time scale 𝑡 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 Fig. 13a) is that of Landau & Levich (1942) defined by Eq. (3.5) describing the meniscus receding at a constant velocity 𝑈 r . In this calculation, 𝑅 m is assumed to remain constant (= 𝑅) because 𝐶𝑎 𝑟 is small. The calculation is performed for 𝐶𝑎 𝑟 = 10 -3 but the results are independent of 𝐶𝑎 𝑟 provided it is small enough. The deposited film thickness is given by Eq. (3.13).

At 𝑡 = 0, the meniscus makes a sudden change of its motion direction and for 𝑡 > 0 advances over the liquid film at a constant velocity 𝑈 a . The film is infinite so the meniscus frame of reference and thus Eq. (3.2) with 𝑈 = -𝑈 a < 0 are employed. With the scaling of Table 1, the dimensionless governing equation is

𝜕 h 𝜕 t + 𝜕 𝜕 x h3 3 𝜕 3 h 𝜕 x3 -Ũ h = 0, ( A 1)
where Ũ = -𝑈 a /𝑈 r . The boundary conditions are Eqs. (3.3, 3.6), where both 𝛼 and ℎ r are now known. Eq. (A 1) is solved numerically, cf. sec. 4.5 for the details. The time evolution of the interface profile is shown in Fig. 13a. At t = 0, the profile is given by the Landau & Levich profile. As the time increases, the film relaxes to the steady advancing profile (the solid curve in Fig. 13a), which coincides with the Bretherton 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 𝑡 rel ; for 𝑈 a /𝑈 r = 1, trel 0.65.

The calculation shows (Fig. 13b) that the relaxation time 𝑡 rel decreases with 𝑈 a and remains smaller than 2. As the exponential vanishes after t 5t rel , one expects that the system behaves independently of a current state after a time lag t 10. This means that, for 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 (Fig. 14a) where the opening angle 𝜃 varies with the angular velocity 𝜔. It is well known that this is equivalent to the streamfunction formulation,

∇ 4 𝜓 = 0, ( B 2)
where For the fixed contact angle case (𝜔 = 0), this equation admits an analytical solution (Janeček & Nikolayev 2012) that satisfies the condition (4.4):

𝑣
Δ𝑝 = Δ𝑝 cl 1 - ℓ K 𝑥 𝐾 1 ℓ K 𝑥 , ( C 3)
where 𝐾 1 (•) is the modified Bessel function of the first order and C.2. Curved wedge flow caused by evaporation for 𝜔 = 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 Eq. (C 3). It can be seen as a curvature that increases when 𝑥 → 0. The curvature creates a difference between the microscopic contact angle 𝜃 micro , the actual slope at the contact line, and the interface slope farther away from the contact line, cf. Fig. 14b. The characteristic length for this effect is ℓ K ∼ 10 -100 nm (Janeček et al. 2013) from the contact line. At a length scale 𝑥 meso ℓ K but much smaller than the film-related length scale ∼ 10 𝜇m, one can define the experimentally measurable interface slope 𝜃, called the apparent contact angle.

ℓ K = √ 3𝜇𝑘𝑇 sat L 𝜌𝜃 2 (C 4)
The region 𝑥 < 0 < 𝑥 meso is often referred to as the microregion. The numerical calculation of 𝜃 is described by Janeček & Nikolayev (2012). It is based on the steady version of Eq. (4.7) (i.e. with 𝜕ℎ/𝜕𝑡 = 0) solved with the boundary conditions (4.9a, 4.9b). The other two boundary conditions for this fourth-order differential equation are the imposed slope 𝜃 micro at the contact line and the condition of zero (on the microregion scale) curvature

𝜕 2 ℎ 𝜕𝑥 2 𝑥=𝑥 meso = 0 (C 7)
Note that the interface slope saturates at 𝑥 ℓ K so 𝜃 is independent of 𝑥 meso . Figure 15 demonstrates an example of 𝜃 as a function of 𝜃 micro and Δ𝑇 for pentane at 1 bar, cf, Eq. (5.1). It turns out that 𝜃 (Δ𝑇, 𝜃 micro ) monotonically grows with both 𝜃 micro and intensity of evaporation controlled by Δ𝑇; evidently, 𝜃 (Δ𝑇 = 0) = 𝜃 micro .

One now introduces 𝜃 evp-min (Δ𝑇) ≡ 𝜃 (Δ𝑇, 𝜃 micro = 0), (C 8) which is the lower bound of 𝜃 (Δ𝑇, 𝜃 micro ), i.e. its value for the complete wetting case (Janeček et al. 2013). Therefore, the 𝜃 micro = 0 curve in Fig. 15 represents 𝜃 evp-min (Δ𝑇). For example, for Δ𝑇 =30 mK, one finds 𝜃 evp-min (Δ𝑇 = 30 mK) 8.4 • , see Fig. 15. This means that, for Δ𝑇 =30 mK, it is impossible to have an apparent angle smaller than 8.4 • .

  Fig. 1. By using Eq. (3.7) in this equation, one finally obtains

Figure 2 :Figure 3 :

 23 Figure 2: Periodic interface shape variation at oscillation for P = 50, 𝐶𝑎 0 = 10 -3 and 𝜃 i = 20 • . The labels give the times corresponding to each profile, and arrows indicate the meniscus motion direction.

Figure 4 :Figure 5 :

 45 Figure 4: Contact angle variation for 𝐶𝑎 0 = 10 -3 , P = 50. (a) Variation of contact angle during an oscillation for different values of 𝜃 i . The inset shows enlarged undulating portions of the curves. (b) The interface shape variation near the contact line around 𝑡 = 0.913𝑃 where 𝜃 min is attained for 𝜃 i = 20 • (cf. the inset to Fig. 4a).

Figure 6 :

 6 Figure 6: Time evolution of the radius of meniscus curvature during an oscillation. The quasi-steady evolution of 𝑅 m,q at receding (𝑡 0.5𝑃) is shown for comparison: (a) 𝑅 m evolution for different P; 𝜃 i = 40 • and 𝐶𝑎 0 = 10 -3 are fixed, (b) 𝑅 m evolution for different 𝜃 i ; P = 150 and 𝐶𝑎 0 = 10 -3 are fixed and (c) 𝑅 m evolution (dashed curves) for different 𝐶𝑎 0 ; P = 150 and 𝜃 i = 20 • are fixed. The solid curves of the respective colourshow 𝑅 m,q (𝑡).

Figure 7 :

 7 Figure 7: Average thickness of the liquid film at 𝑡 = 0.5𝑃; ℎ q is shown for comparison. (A) Value of h as a function of oscillation period for different 𝜃 i and for 𝐶𝑎 0 = 10 -3 and (b) value of ℎ as a function of 𝐶𝑎 0 for different P and for 𝜃 i = 20 • .

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Liquid plug velocity and left meniscus radius variations during one period of the Lips et al. (2010) 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 𝑅 m,q variation given by Eq. (3.11) at receding (𝑡 0.5𝑃) is also shown.

Figure 11 :

 11 Figure 11: Liquid film thickness: comparison of the experimental data of Youn et al.(2018), quasi-steady estimation with Eq. (3.12) and numerical results for water.

Figure 13 :

 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 receding with the velocity 𝑈 r to steady advance with the velocity 𝑈 a such that |𝑈 a | = 𝑈 r . The time lag between the two curves is Δt = 0.5. The inset shows the root mean square deviation from the steady advancing profile and (b) relaxation time as a function of 𝑈 a /𝑈 r . The dots are the numerical points. The line is the exponential fit.

Figure 14 :

 14 Figure 14: Sketch of the straight liquid wedges with varying angle (a) and curved wedge at evaporation (b).

  admits a solution 𝜓 ∝ 𝑟 𝜆 𝑓 (𝜑), where 𝜆 is a constant(Moffatt 1964). On the one hand, 𝜓 should scale with 𝜔 because it causes the flow. On the other, from Eq. (B 4), its dimension should be length 2 /time. The only choice is thus𝜓 = 𝜔𝑟 2 𝑓 (𝜑). (B 5)The corresponding function 𝑓 is 𝑓 (𝜑) = 𝐴 cos 2𝜑 + 𝐵 sin 2𝜑 + 𝐶𝜑 + 𝐷 (Moffatt 1964), where 𝐴, 𝐵, 𝐶, 𝐷 are constants. They can be determined with the boundary conditions. A zero velocity (𝑣 𝑟 = 𝑣 𝜑 = 0) is imposed at the liquid-solid surface 𝜑 = 0. A radial velocity 𝑣 𝜑 = 𝜔𝑟 is imposed at the liquid-vapour boundary 𝜑 = 𝜃, which is also stress free (𝜕𝑣 𝑟 /𝜕𝜑 = 0). The result is𝜓 = 𝜔𝑟 2 tan 2𝜃 (cos 2𝜑 -1) + 2𝜑sin 2𝜑 2(tan 2𝜃 -2𝜃) 𝑝 𝑙 ∼ log 𝑟.Note that for a small 𝜃, Eq. (B 7) reduces to Straight wedge flow Consider first the case with no phase change. It is described by Eq. (2.1), where 𝐽 = 0. As in the above Stokes problem (Appendix B), one can consider asymptotically (near the contact line) the straight wedge ℎ = 𝜃𝑥 with the angular velocity 𝜕𝜃/𝜕𝑡 = 𝜔. The governing equation , this expression agrees with the small 𝜃 asymptotics (B 8) 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, Eq. (4.7). For the straight wedge with the varying contact angle it becomes

  the case of varying contact angle 𝜔 ≠ 0, a solution that satisfies the condition (4.4) can be found as an asymptotic expansionΔ𝑝 = Δ𝑝 cl + 3𝜇𝜔 𝜃 3 ℓ 2 K 𝑥 2 + O (𝑥 4). Value of 𝜃 as a function of Δ𝑇 for different 𝜃 micro computed for pentane at 1 bar. The curve for 𝜃 micro = 0 corresponds to 𝜃 evp-min (Δ𝑇).

Table 1 :

 1 Dimensional reference values that are used to make the governing equations non-dimensional.

Table 2 :

 2 Fluid properties at the experimental conditions (1bar, 20 • C) of Lips et al. (2010); and key dimensionless numbers.

	parameter	notation value
	tube inner radius	𝑅	1.2 mm
	density	𝜌	625.7 kg/m 3
	surface tension	𝜎	0.0152 N/m
	shear viscosity	𝜇	2.37 • 10 -4 Pa•s
	reference velocity (velocity amplitude)	𝑈 0	0.24 m/s
	oscillation frequency		3.7Hz
	dimensionless period	P	260.1
	capillary number	𝐶𝑎 0 3.74 • 10 -3
	Reynolds number	Re 0	1521
	Weber number	𝑊 𝑒 0 5.66

Table 3 :

 3 3.12) and numerical results for water. Fluid properties at the experimental conditions (1bar, 25 • C) of Youn et al.

	parameter	notation		value
	tube inner radius	𝑅		0.5 mm
	fluid		water	ethanol
	surface tension	𝜎	7.20 • 10 -2 N/m	2.23 • 10 -2 N/m
	shear viscosity	𝜇	8.88 • 10 -4 Pa•s	1.088 • 10 -3 Pa•s
	liquid density	𝜌	997 kg/m 3	785 kg/m 3
	oscillation frequency		2 Hz	6 Hz	2 Hz	6 Hz
	oscillation amplitude	𝐴 0	19.5 mm	22.5 mm 20.69 mm 25.2 mm
	𝑈 amplitude in simulation	𝑈 0 = 2𝜋 𝐴/𝑃 0.245 m/s 0.848 m/s 0.26 m/s 0.950 m/s
	𝑈 amplitude in experiments		0.254 m/s 1.025 m/s 0.285 m/s 1.150 m/s
	dimensionless period	P	1267	966	834	659
	capillary number	𝐶𝑎 0	0.0030	0.0105	0.0127	0.0463
	Reynolds number	Re 0	275.07	952.09	187.59	685.43
	Weber number	𝑊 𝑒 0	0.8312	9.9576	2.3796	31.77

  In polar coordinates (𝑟, 𝜑), the Stokes equations for the liquid velocity ì 𝑣 = (𝑣 𝑟 , 𝑣 𝜑 ) read

	𝜕 𝑝 𝑙 𝜕𝑟	=𝜇	𝜕 𝜕𝑟	1 𝑟	𝜕 (𝑟𝑣 𝑟 ) 𝜕𝑟	+	1 𝑟 2	𝜕 2 𝑣 𝑟 𝜕𝜑 2 -	2 𝑟 2	𝜕𝑣 𝜑 𝜕𝜑	(B 1a)
	𝜕 𝑝 𝑙 𝜕𝜑	=𝑟 𝜇	𝜕 𝜕𝑟		1 𝑟	𝜕 (𝑟𝑣 𝜑 ) 𝜕𝑟		+	1 𝑟 2	𝜕 2 𝑣 𝜑 𝜕𝜑 2 +	2 𝑟 2	𝜕𝑣 𝑟 𝜕𝜑	(B 1b)
	∇ • ì 𝑣 =0.										(B 1c)

† Email address for correspondence: vadim.nikolayev@cea.fr

Acknowledgments

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 program, which has received funding from the European Union Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement No. 800945. An additional financial support of CNES awarded through GdR MFA is acknowledged. We are indebted to S. Lips for opening the access to his experimental results.

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 𝑡 rel . It is a characteristic time scale of decrease of