Universal oscillatory dynamics in capillary filling

Capillary filling action, characterized by its passive pumping ability, enjoys the focus of diverse fields of applications ranging from natural science [1–10] to engineering technology [11–19]. With recent focus on modulating reduced scale flows [20–26], surface force generated capillarity transport, thus, forms a favourite tool in cutting-edge microfluidic research [27,28]. Genesis of capillarity holds its inception in the intrinsic interfacial energies of the solid (S), fluid (L) and the vapor (V) phase. While a large number of reported studies on capillary filling dynamics aptly justify its outstanding relevance, the proposition of depicting a universal picture of the oscillatory characteristics of the capillary front, typically observed at the onset of reaching the equilibrium meniscus penetration, remains elusive. This deficit stems from the complexities in generalizing the underlying physical phenomena with the aid of unified normalizing physical parameters.

Here, we bring out a universal depiction of the oscillatory dynamics of capillary filling processes, by appealing to the relevant normalized physical parameters pertinent to Newtonian fluids. We show that the resultant data sets indeed collapse into a signature master curve expressed as a single-valued function of the significant normalized parameters, bearing significant implications towards realizing a generalization in the design principles, thereby encompassing a wide variety of applications ranging from nature to energy. Our comprehensive experiments and computational analysis provide validations to this apparently non-intuitive theoretical proposition.

The most simplistic approach towards predicting capillary filling dynamics follows the classical Lucas-Washburn approach with added inertial effect [1,29,30]. The formulation proceeds with a lumped balance of inertial (F_I), surface tension (F_S), viscous (F_V) and gravitation (F_G) forces and its dynamical analysis. Considering a representative cylindrical capillary of radius R $(r\text{-}\theta\text{-}z\,\text{coordinate system})$, with the origin located at the capillary entry, the inertial, surface tension and gravitational forces for a capillary rise of z can be written as $F_{I} =\frac{\text{d}}{\text{d}t}({\rho \pi R^{2}z\frac{\text{d}z}{\text{d}t}})$, $F_{S}=2\pi R\sigma \cos(\theta)$ and $F_{G}= \pi R^{2}z\rho g$, respectively (see schematics in fig. 1). Here ρ, σ, θ and g represent fluid density, surface tension at the $L\text{-}V$ interface, contact angle at the $S\text{-}V\text{-}L$ contact line and gravitational acceleration, respectively. The viscous force may be derived as a first estimation considering a fully developed laminar flow: $F_{V}=8\pi z\eta_{f} \frac{\text{d}z}{\text{d}t}$ (see the Supplementary Material Supplementarymaterial.pdf (SM)). The balance of these forces essentially predicts the capillary entry dynamics wherein the effect of the vapour phase is usually neglected [3,31].

Zoom In Zoom Out Reset image size

The simplified analysis as described above lacks the consideration of viscous drag $F_{\textit{VD}} =8\pi ({L-z})u_{avg} \eta_{g}$ in the vapour phase (which becomes critical for less viscous fluids and longer channels. [32]) and includes some further idealizations in capturing the details of the capillary filling process. These considerations are as follows: a) loss in the kinetic energy during the downward motion of the capillary front when it is undergoing oscillatory dynamics; b) loss at the entrance region; c) errors incurred due to lack of capturing the capillary meniscus topography; and d) errors for not considering the dynamic contact angle. The loss of kinetic energy can be accounted for through the inclusion of a parameter b (detailed below). The other factors b), c) and d) may be neglected in the model without sacrificing the essential physics of interest. For further details on accounting for the effects from gas viscosity and kinetic energy loss, reported literature can be referred [33,34]. Finally, the modified force balance is non-dimensionalized using the variables $\bar{t}=t/{\sqrt {{\pi R^{3}}/{2\sigma}}}$ and $\bar{z}=z / R$. The equivalent dimensionless governing equation for capillary rise, thus, reads (see SM for the derivations)

$$\begin{equation} \frac{\text{d}^{2}\bar{z}}{\text{d}\bar{t}^{2}}=\frac{\cos \theta }{\bar{z}}-Bo-Ca_{f} \frac{\text{d}\bar{z}}{\text{d}\bar{t}}-\frac{Ca_{g} }{\bar{z}}({\bar{L}-\bar{z}})\frac{\text{d}\bar{z}}{\text{d}\bar{t}}-\frac{b}{\bar{z}}\left( {\frac{\text{d}\bar{z}}{\text{d}\bar{t}}} \right)^{2}, \end{equation} \tag{ 1 }$$

where $Bo={\rho gR^{2}}/ {2\sigma }$, $Ca_{f,g} ={4\eta_{f,g} R}/{\sigma t_{0}}$ and $\eta_{f,g}$ denote the dimensionless Bond number, capillary number and viscosity, while the subscripts f and g denote the quantities for the fluid and vapour phases, respectively.

As mentioned above, during the downward motion of the capillary front, which is undergoing an oscillatory dynamics, the kinetic energy contribution of the liquid column is lost in the effective infinite fluid reservoir [33,34]. This is reflected through the parameter b which holds the value $b =1$ when ${\text{d}\bar{z}}/{\text{d}\bar{t}}>0$ (representing the upward motion of the meniscus through a positive velocity) and $b =0$ when ${\text{d}\bar{z}}/{\text{d}\bar{t}}<0$ (representing the downward motion of the meniscus through a negative velocity). Thus, with b = 0, the pure kinetic energy part of the inertial energy, that is now lost in the infinite fluid reservoir, is not accounted for in the resultant force balance while only the instantaneous acceleration of the liquid column contributes to the force balance during the downward motion. From the above governing equation, the dimensionless equilibrium height (or Jurin height, $\bar{J}$) of the capillary rise, obtained by neglecting the viscous and inertial contributions at the equilibrium state, reads $\bar{J}={\text{cos}\theta}/{Bo}$. Now, we proceed to obtain the criteria of oscillation during the rise of the capillary front. Therefore, starting from eq. (1) and setting $b =1$, we make the two consecutive transformations via $\bar{z}=\bar{J}-\bar{z}$ and $X={\bar{z}}/{\bar{J}}$ which physically connote the idea of focusing on the dynamics near a magnified region entered around the Jurin height. These transformations recast eq. (1) in the form

$$\begin{equation} \frac{\text{d}^{2}X}{\text{d}\bar{t}^{2}}+A\frac{\text{d}X}{\text{d}\bar{t}}+BX=0 \end{equation} \tag{ 2 }$$

that may be recognized as the classical second-order equation for damped oscillator where $A=Ca_{M}$ and $B={Bo}/J={Bo^{2}}/ {\cos \theta }$, where $Ca_{M} =Ca+Ca_{g} ( {{\bar{L}}/ {\bar{J}}-1})$ (details of derivation are in the SM). Equation (2) has a solution of the form $X=\exp ({m\bar{t}}\,)$ with the $m={({-A\pm \sqrt {A^{2}-4\omega^{2}}})}/2$. Since the oscillation criterion is given by $\sqrt {A^{2}-4\omega^{2}}$ where $\omega =\sqrt B$, one may quickly infer from here that the criterion for the inception of the oscillatory regime of the capillary dynamics is $Bo>{Ca_{M} \sqrt {\cos \theta } }/2$. Thus, if the system Bond number exceeds a critical Bond number $Bo_{c} ={Ca_{M} \sqrt {\cos \theta } }/ 2$, the system exhibits oscillatory dynamics. Drawing analogies from the dynamics of vibration [35], one finds that in such oscillatory regimes, the ratio of the initial peak overshoot to the equilibrium height (similar to $M={z_{p} }/J$, where z_p is the maximum height attained by the capillary front during the total oscillatory dynamical process) is given by

$$\begin{equation} M=\exp \Bigl( {\sqrt {{-\pi \xi }/{\sqrt {1-\xi^{2}} }} }\Bigr), \end{equation} \tag{ 3 }$$

where $\xi ={Ca_{M} \sqrt {\cos \theta } }/{({2Bo})}$ is the damping factor. The analysis reveals that the argument $\xi/{\sqrt {1-\xi^{2}} }$ for the function M is only a function of ${Bo}/{Bo_{c} }$ (see SM) and, therefore, a plot ${z_{p} }/J$ against ${Bo}/{Bo_{c}}$ becomes independent of Ca_M and θ. Considering that one may construct a master curve wherein systems with a wide variety of problem data converge to be represented by a single plot, this realizes an effective generalization. Although not apparent, the dynamical plots related to an annular geometry [36] or a parallel-plate geometry can also be associated with similar master curves through proper non-dimensionalization scheme, thereby, delineating the generalization of the present formulation. Comprehensive experiments and numerical investigations indeed support this apparently non-intuitive conjecture that is reached from an elusively simple theoretical formalism.

Our experiments (schematics shown in fig. 1) are performed using capillaries of diameters of the orders of millimeters, while different fluids (with their relevant properties are given in table S1a in the SM) are used in order to explore the various dimensionless numbers appearing in the theoretical formulation. In order to explicitly visualize the oscillation in the filling dynamics, a timeline plot for one experimental set is provided in fig. 2(a). For our experiments, we have used n-heptane, n-hexane, acetone and n-pentane (Merck, 99%) as the wetting fluids and atmospheric air as the non-wetting fluid. Borosilicate glass capillary tubes of 1.1 mm diameter are purchased from Spectrum Marketing (AQCAPS). A cleaning protocol is followed for the capillary tubes where they are first rinsed in isopropyl alcohol and then subjected to ultrasonication for 30 minutes to remove any entrapped moisture [37–40]. Then it is flushed and rinsed with Milli-Q ultra pure water (Millipore India Pvt. Ltd) which is followed by a further ultrasonication of 30 minutes. To enhance the accuracy and reproducibility of the experimental outcomes, each capillary tube is used only once to avoid any chance of contamination. The inner diameters of the capillary tubes are measured from microscopic images of the capillary tube cross-section using standard image processing techniques.

Zoom In Zoom Out Reset image size

The spatial movement of the capillary tube is controlled using a high-precision $x\text{-}y\text{-}z$ translation stage (Thorlabs Inc., Newton, New Jersey, United States). The tube is attached to a holder and mounted vertically on the translation stage. The fluid is filled in a glass Petri dish of large diameter to minimize the curvature effects. At the start of the experiment, the capillary tube is lowered slowly towards the Petri dish up to the point of first contact at the interface, from which the rise is measured $(t=0, h=0)$. To capture the time evolution of the capillary meniscus, a high-speed CCD camera (IMPERX: VGA 210), fitted with a macro-lens, is employed along with a diffuse light source. The images are recorded at 207 frames per second at a resolution of $640 \times 480$ pixels and the post-processing is carried out using ImageJ software (NIH, MD, USA). The fluids used for the experiments, namely n-heptane, n-hexane, acetone and n-pentane, are known to be volatile in nature, and thus, changes the environmental temperature of the experiments. Since viscosity, density and surface tension are very sensitive to temperature variations, we sought to measure the temperature variation in the experimental conditions. With the ambient temperature around 30 °C, the temperature on the volatile liquid surface on the Petri dish is measured to be 25 °C. We obtained the relevant physical properties of the testing fluids, that is, density $(\rho)$, surface tension $(\sigma)$ and viscosities $(\eta)$ at this temperature and tabulated them with the corresponding references (table S1a in the SM). Moreover, we show that previous experimental data also conforms to our claim (see fig. 3 and table S1b in the SM).

Zoom In Zoom Out Reset image size

For the computational simulation of capillary rise in the two-phase $L\text{-}V$ system, we employ the phase-field model. Phase-field formalism has been extensively used in studies relating to two-phase systems [41] and finds remarkable matching with analytical [42,43] or experimental results [42,44]. We directly report the dimensionless form for the governing phase-field equations. The dynamics of the phase-field parameter ϕ is given using

$$\begin{equation} \frac{D\phi}{Dt}=\frac{1}{Pe}\nabla^{2}\mu, \end{equation} \tag{ 4 }$$

where the chemical potential reads as

$$\begin{equation} \mu =-Cn\nabla^{2}\phi +\frac{1}{Cn}\phi ({\phi^{2}-1}), \end{equation} \tag{ 5 }$$

where Cn and Pe are the dimensionless Cahn number and Peclet number, respectively. The corresponding Navier-Stokes equation governing the two-phase fluid flow reads

$$\begin{equation} Re\rho \frac{D\mathbf{u}}{Dt}=-\nabla p+\nabla \cdot ( {\eta [{\nabla \mathbf{u}+\nabla \mathbf{u}^{\text{T}}}]})+\frac{4}{Ca}\mu \nabla \phi +\frac{8Bo}{Ca} \end{equation} \tag{ 6 }$$

along with the continuity of both the phases $\nabla \cdot\mathbf{u}=0$. For full details of the formulation, dimensionless numbers and symbols, refer to the SM. For the system domain and corresponding boundary conditions, refer to fig. 1(b) above. We have considered a 2D axisymmetric domain for the 3D simulation which is shown in the SM for visual clarity. The red portion is the cylindrical tube where initially the rarer phase (air) is considered to be present. The Petri-dish–like structure, with blue color, initially houses the denser phase (water) which fills the capillary tube due to the capillarity phenomena. The inlet has the natural pressure gradient due to gravity while the outlet is open to the atmosphere. The other walls (c) have a wetted boundary condition with a specified contact angle. The governing equations along with the boundary relations are solved using the commercial software COMSOL Multiphysics.

Comparison of the experimental and numerical results is pictorially depicted in figs. 2(a) and (b). In fig. 2(a), we delineate the experimental timeline snapshot of the capillary front evolution with time. The figure clearly demonstrates the oscillatory nature of the contact line which eventually attains equilibrium height at large times. Such oscillations are also demonstrated via the phase-field modelling as shown in fig. 2(b) which shows the evolution of the $L\text{-}V$ interface. Since the numerical solutions are solved in axisymmetric configuration, we represent this axisymmetric domain with the idea that it can be easily rotated about its symmetric axis to showcase the 3D viewing of the system. The details of the physical properties used for the experiments and simulation are provided in the caption of figs. 2(a), (b).

Here we begin by describing the physical origin of the oscillations. Although the equilibrium position of the meniscus represents the balance of the gravitational and surface tension forces, the moment when the meniscus approaches the equilibrium position for the first time (when no oscillation has initiated yet), there may be residual inertia forces (due to the initial high kinetic entry energy which is not yet damped by the hindering viscous losses) which propels the meniscus to overshoot the equilibrium position to attain a position that, by virtue of force balance at that section, is unstable. Therefore, the meniscus motion ceases at the point where gravitational forces start to exceed the inertial and surface tension effects and from thereon the meniscus starts to drop due to the action of the gravity. Here the viscous force, by virtue of its nature, opposes the downward motion. Again, the meniscus may have picked up enough inertial force from the gravitational forces that it crosses the equilibrium position. Thus, the phenomenon of oscillation about the Jurin height is observed until the viscous forces quickly dampen out the inertial energy near the equilibrium position. In retrospect, oscillations are not visible for fluids with high viscosity or in channels with very low hydraulic diameter where the viscous effects dominate over inertial effects that are inherently responsible for the occurrence of oscillations. These oscillations are observed in fig. 2(a) (experimental results) and fig. 2(b) (numerical results). An enlarged view around the Jurin height for experimental results and a corresponding video (see supplementary movie Movie_S1.avi) will clearly showcase the above physical expression. For the timeline depiction for the capillary rise in the complete channel, refer to fig. S2b in the SM. Figure 2(b) depicts the phase-field simulation results that have also captured the oscillatory dynamics of the system. The qualitative nature of the oscillations seems similar, which is about the Jurin height. Although these figures present an interesting qualitative understanding of the capillary front oscillations, the intricacies in the underlying physics need detailed attention through proper analysis of the phenomenon and quantification of the same. These will help us to build the base for comprehensive illustration of the characteristic universality that subtly exists across various conditions of the phenomenon.

In fig. 3(a) we attempt to reproduce the experimental observation by employing our analytical framework. The experimental data (depicted using markers) agree reasonably well with the results obtained from the corresponding analytical formalism (depicted by solid lines; the oval structure groups the experimental plot with the corresponding analytical plot). The plots in the main figure reflect in-house performed experimental analysis while those in the inset are data obtained from previous studies. Finally, the plot of group 5 demonstrates the matching between the numerical data and the analytical formulation. All these data points depict a clear impression of the oscillatory nature of the capillary front.

With reasonable good agreement among the generalized analytical approach, comprehensive numerical analysis and experimental data, we next explain the generalization in the master curve. In this curve, plots with blue, green and red markers describe the variation in ${\bar{z}_{p} } / {\bar{J}}$ as a function of ${Bo}/{Bo_{c}}$. With a given choice of this combination of the coordinate axes, all the data points collapse into a single unique characteristic curve. However, for any other combinations of the coordinate axes, such a master curve will not be observed (as detailed above and explicitly demonstrated in the SM).

One subtle aspect of the master curve shows that as long as ${Bo}/ {Bo_{c} }\le 1$, ${\bar{z}_{p} }/{\bar{J}}$ remains equal to unity. It shows that until the Bond number exceeds the critical Bond number, the maximum peak attained by the capillary front remains equal to the Jurin height; or, in other words, no oscillation of the capillary front is observed. As soon as ${Bo}/ {Bo_{c} }>1$, ${\bar{z}_{p} }/{\bar{J}}$ becomes greater than unity showing the oscillatory nature of the curve. At large ${Bo}/{Bo_{c}}$ ratios, the oscillation of the capillary front becomes very prominent. The capillary number and contact angle independence of the master curve is not a mere artefact of the analytical formalism but is also reflected in the realistic situations. This is delineated from the experimental and numerical data in a sense that, although performed in completely different dimensional scales, the resultant characteristics closely coincides with the master curve. In the literature, a criterion of oscillation was obtained numerically from numerous simulation results [34] which showed that above a critical ω, expressed by $\omega =({{\rho^{3}g^{2}R^{5}}/{128\eta_{f}^{2} \sigma \cos \theta }})$, oscillations start occurring. They found that $\omega_{cr}=0.25$. Now, referring back to the present work, we find that the oscillation criterion is $Bo_{c} ={Ca_{M} \sqrt {\cos \theta } }/2$ where the critical Bo is dependent on the capillary number of the flow. Rearranging in terms of the physical parameters $Bo={\rho gR^{2}}/{2\sigma }$, $Ca_{M} =Ca_{f} ={4\eta_{f} R}/{\sigma t_{sc} }$, $t_{sc} =\sqrt {{\rho R^{3}}/{2\sigma }}$, we find $({{\rho ^{3}g^{2}R^{5}}/{128\cos \theta \sigma \eta_{f}^{2} }})=1/4$. Now, the expression on the left-hand side is identical to the parameter ω defined by Masoodi et al. [34]. Thus, $\omega_{cr} =0.25$ (since this form came from the criterion equation, we have replaced ω by critical $\omega~(\omega_{cr})$).

The slight discrepancy between the analytical and experimental results may be attributed to the fact that phenomena such as entrance effects and capillary interface geometry are not accounted for in the simplistic theoretical analysis [32]. The entrance effect usually acts in a way to dampen the initial inertial energy that the liquid column may gain from the actuating surface tension forces. This is reflected in the first oscillatory dynamics where the maximum attained Jurin height is reduced as compared to that anticipated from the present formalism. Further, such an effect will cause the capillary entry rate to reduce, thereby decreasing the slope of the z vs. t plot (details in the SM). Interestingly, this is exactly the value of critical ω that was obtained numerically, and, thus, it establishes the consistency and universality of our formulation. However, it must be noted that their work has neglected the effect of the viscous drag, which may need to be accounted for while predicting the experimental or numerical trends.

Thus, the actual peak overshoot observed remains slightly below the one predicted analytically (fig. 3(a)). Nevertheless, we comment here that including such higher-order complexities in the analytical solution will compromise on the elegant nature of the theoretical paradigm without adding much to the essential physics governing the oscillatory regimes about the equilibrium (Jurin) height.

As mentioned above, with the choices of the axes ${\bar{z}_{p}}/{\bar{J}}$ and ${Bo}/{Bo_{c} }$, the plot portrays a universal depiction of the oscillatory dynamics of the capillary filling phenomenon by collapsing all data, as aptly justified by eq. (3). Our analytical derivations of the dimensionless system characteristics justify that other ways to group the dimensionless number may fail to generate the master curve under the present configuration. To follow up the logic, we proceed with an example to plot with different groupings of the dimensionless numbers. As seen in fig. 3(c), a plot with $\bar{z}_{p}$ against ${Bo}/{Bo_{c}}$ fails to reproduce the universality nature of the problem since the curves show to vary with Ca_M. Again in the inset of fig. 3(c), we consider the other choices of the dimensional groups, namely, ${\bar{z}_{p}}/{\bar{J}}$ plotted against Bo. Yet again, we find that the above choice fails to collapse the different Ca_M data in a single universal curve.

To summarize, we have discovered a unique universal characteristic of the oscillatory regime observed in capillary filling processes of Newtonian liquids. With a suitable normalization of the relevant parameters having inception from a simple analytically based theory, we show that the resultant generalization can collapse the problem data in the form of a single master curve that remains invariant across experiments. This may turn out to be of profound importance in developing a universal design paradigm for engineered capillary systems having a wide variety of applications ranging from oil recovery to medical diagnostics.