Capillary filling in microchannels with wall corrugations -
A comparative study of the Concus-Finn criterion by continuum, kinetic and atomistic approaches

Our CFD approach is based on the CFD-ACE+ software package (release 2008), as a multiphysics commercial tool including geometry modeling, grid generation and results visualization [34]. The spontaneous capillary filling in micro/nano channels is reproduced by means of the VOF (Volume of Fluid) scheme, based on the hypotheses of incompressibility, no-interpenetration and negligible localized relative slip of the two fluids in contact. In order to describe the transport of the volume fraction $\phi$ of one of the two fluids in each cell ($\phi$ thus ranging from 0 to 1), the Navier-Stokes equations for the fluid velocity are augmented with a scalar transport equation for the volume fraction:

where $t$ is time, $\vec{\nabla}$ is the standard spatial gradient operator, and $\bm{u}$ is the velocity vector of the fluid. The composition of the two fluids in the mixture determines for each computational cell the averaged value of physical properties such as density and viscosity. Any volume-specific quantity is evaluated in accordance with the following expression:

where $\bar{\omega}$ is the volume-averaged quantity, $\omega_{l}$ (resp. $\omega_{g}$) is the value of the property for one liquid (resp. gas). However, being $\phi$ the averaged volume fraction of fluid in each cell, the definition of the interface between the two fluids in that cell is not unique and must be dynamically reconstructed throughout the simulation. In CFD-ACE+, this is accomplished through an upwind scheme with PLIC (Piecewise Linear Interface Construction method [35, 36]), taking into account surface tension to determine accurately the interface curvature. As to boundary conditions, a zero static pressure is imposed at both inlet (liquid only) and outlet (gas only). At the initial time both fluids are at rest, with a short portion of the channel at the inlet filled with liquid. This specific configuration allows the liquid-vapour interface to assume the curvature corresponding to the hydrophilic partial wetting condition imposed at walls through the contact angle $\theta$, thus generating the capillary pressure that drives the fluid motion. The CFD simulation set-up consists of a $2D$ straight channel with a height $H=800\;$nm and overall length $L=80\;\mu$m. At the bottom wall, a squared post, of side and height $h=400\;$nm, is located. The computational domain has been discretized with a squared structured non-uniform grid consisting of $185000$ cells and $190000$ nodes.

The Lattice Boltzmann (LB) method is based on a minimal form of the Boltzmann kinetic equation describing the time evolution of discrete particle distribution function $f_{i}(\bm{x},t)$, denoting the probability of finding a particle at lattice site $\bm{x}$ and time $t$ moving along the lattice vector $\bm{c}_{i}$ ( 2), where the index $i$ labels the discrete directions of motion.

In mathematical terms [27], the LB equation reads as follows: [37]

where $\rho$ is the fluid density, ${\bm{u}}$ is the fluid velocity and $i=0-8$ labels the the nine-speed, two-dimensional $2DQ9$ model [27]. The second term on the right hand side is a model collision operator $\Omega_{i}$, describing the relaxation towards a local equilibrium on a time scale $\tau$. Finally, $F_{i}$ describes the effect of external/internal forces on the discrete particle distribution. The macroscopic density and momentum are obtained as weighted averages of discrete distributions: [27]

For the case of a two-phase fluid, the interparticle force reads as follows:

Here $\psi({\bm{x}},t)=\psi(\rho({\bm{x}},t))=(1-\exp(-\rho(\bm{x},t)))$ is the pseudo-potential functional, describing the fluid-fluid interactions triggered by inhomogeneities of the density profile, and ${\cal G}$ is the strength of the coupling (see [11, 38] for details). Finally, $F_{i}=w_{i}\;{\bm{F}}\cdot{\bm{c}}_{i}/c_{s}^{2}$, where $w_{i}=0,1/9,1/36$ are the standard weights of the nine-speed two-dimensional lattice 2DQ9 [27] and $c_{s}=1/\sqrt{3}$ is the sound speed in lattice units. Besides introducing a non-ideal excess pressure $p^{\star}=\frac{1}{2}{\cal G}c_{s}^{2}\psi^{2}$, the model also provides a surface tension $\gamma\sim{\cal G}c_{s}^{4}\frac{(\delta\psi)^{2}}{\delta_{w}}$, where $\delta\psi$ is the drop of the pseudo-potential across the interface of width $\delta_{w}$. By tuning the value of the pseudo-potential at the wall $\psi(\rho_{w})$, this approach allows to define a static contact angle, spanning the full range of values $\theta\in[0^{o}:180^{o}]$ [39].

As to the boundary conditions, the standard bounce-back rule is imposed: any flux of particles that hits a boundary simply reverses its velocity so that the average velocity at the boundary is automatically zero. This rule can be shown to yield no-slip boundary conditions up to second order in the Knudsen number in the hydrodynamic limit of single phase flows.

In this work we consider a $2D$ channel composed of two parallel plates separated by a distance $H=40\Delta$, where $\Delta$ is the space discretization unit. This two-dimensional geometry, with length $2L=3000\Delta$, is divided in two regions, as shown in 3. The left part has top and bottom periodic boundary conditions, so as to support a perfectly flat gas-liquid interface, mimicking a “infinite reservoir”. The actual capillary channel resides on the right half, of length $L$. The top and bottom boundary conditions are those of a solid wall, with a given contact angle $\theta$. Periodic boundary conditions are also imposed at the west and east sides. The squared obstacle, of side $h=H/2=20\Delta$, is placed on the lower wall.

During the last few decades Molecular Dynamics has proved to be one of the most broadly used simulation techniques, capable of reproducing many details of macroscopic fluid dynamics. It has been extensively used to study dynamic wetting problems and to shed more light on the molecular processes close to the contact line region [40, 41]. An appealing feature of MD (e.g., compared with Monte Carlo or molecular statics) is that it follows the actual dynamical evolution of the system.

We perform MD simulation of a simple generic model on a coarse-grained level. In our MD simulation, see Fig. 4, the fluid particles interact with each other through a Lennard-Jones (LJ) potential,

where $\epsilon=1.4$ and $\sigma=1.0$. The capillary walls are represented by particles forming a triangular lattice with spacing $1.0$ in units of the liquid atom diameter $\sigma$. The wall atoms may fluctuate around their equilibrium positions, subject to a finitely extensible non-linear elastic (FENE) potential

with $R_{0}=1.5$ and $\epsilon_{w}=1.0k_{B}T$ where $k_{B}$ denotes the Boltzmann constant and $T$ is the temperature of the system; $r$ is the distance between the particle and the virtual point which represents its equilibrium position in the wall structure. The FENE-potential acts like an elastic string between the wall-particles and their equilibrium positions in the lattice and keeps the structure of the wall densely-packed hexagonal. In addition, the wall particles interact with each other via a LJ potential with $\epsilon_{ww}=1.0$ and $\sigma_{ww}=0.8$. This choice of interactions guarantees no penetration of liquid through the wall, and at the same time the mobility of the wall particles corresponds to the system temperature.

Molecules are advanced in time via the velocity-Verlet algorithm with integration time step $\delta t=0.01t_{0}$ where $t_{0}=(\sigma^{2}m/48\epsilon_{LJ})^{1/2}=1/\sqrt{48}$ is the basic time-unit and we have taken particle mass $m=1$ and $k_{B}T=1$. The temperature is maintained by a Dissipative-Particle-Dynamics (DPD) thermostat, with friction parameter $\xi=0.5$, thermostat cutoff $r_{c}=2.5\sigma$ and step-function-like weight functions [42, 43]. Fluid properties and flow boundary conditions arise then as a consequence of the collision dynamics and the local friction controlled by the DPD thermostat. An advantage of the DPD method in comparison with other MD schemes is that the local momentum is conserved, so that the hydrodynamic behavior of the liquid at large scales is correctly reproduced.

It is worth emphasizing that all contact angles $\theta$, used in the simulations, have been determined by measuring the mean static contact angle $\theta$ of a sessile meniscus in a slit, formed by two parallel atomistic walls for a few given liquid -wall interactions $\epsilon$. Specific contact angles are then chosen by interpolation between the relevant values of $\epsilon$. Therefore, the actual value of the contact angle when the fluid front is located at an obstacle edge, for instance, may incidentally deviate from the nominal value of $\theta$ . A recent MD study [44] has shown that the LW-law (Eq.(1)), holds almost quantitatively down to nanoscale tube diameters. In this study we use an obstacle shaped as a rectangular ridge which runs perpendicular to the flow direction in the slit and has the same atomic composition as that of the slit walls. The height and width of the ridge are chosen as $10\sigma$ so that the obstacle height takes half of the slit thickness.

As is well known, the Lucas-Washburn asymptotic regime sets in once inertial effects die out and steady propagation results from the sole balance between capillary drive and dissipation. During the transient regime, the dynamic contact angle, which depends on the front speed itself, changes in time until its static value is attained [46, 47].

In 4, the time evolution of the centerline ($y=H/2)$ front position is reported as a function of time, and compared with the dimensionless LW-law Eq. (2). As one can see, all three methods exhibit satisfactory agreement with the theoretical LW prediction, although on a different range of time-scales. This is due to the different values of the parameters used in each method, and particularly, to the fact that the MD capillary speed $V_{cap}=\gamma/\eta$ is $20\div 30)$ times higher than in the other methods. Indeed, it should be noted that the LW asymptotic solution sets in after a typical transient time of the order of $\tau\sim H^{2}/12\nu_{l}$, or in units of capillary times, $\tau/t_{cap}=\frac{\rho\gamma H}{12\mu^{2}}$. Based on the data in Table I, one can readily check that $\tau/t_{cap}\sim 1$ for LB and CFD, and $\tau/t_{cap}\sim 10^{-2}$ for MD. It is therefore possible to appreciate the anomalous transient in the CFD case, see 4, which is of the order $\tau/t_{cap}\sim 100$. However, since CFD and LB simulations last over $10^{3}$ capillary times, and MD simulations last about $10$ capillary times, before reaching the obstacle, the superposition with the LW regime is free from transient phenomena.

In 5, 6 and 7, snapshots of the fluid front during the surmounting of the obstacle are reported for the three simulation methods (the figures report density contours). As one can see, the front dynamics and morphology are strongly affected by the presence of the obstacle: the liquid impinging on the obstacle must adjust to a $90^{\circ}$ degrees discontinuity of the contact angle, which is clearly causing a significant deformation of the front, before the static value is recovered again on the flat-top surface of the obstacle. These changes of shape are well visible in figures 5 and 6, for both CFD and LB. The case of MD is less clear-cut, due to the absence of a well-defined interface and to molecular fluctuations. Indeed, although the fluid meniscus is clearly visible, its surface appears rough and strongly fluctuating in time. From 7 one can see that some atoms evaporate from the liquid and overcome the slit. Moreover, long before the fluid meniscus has passed the obstacle, vapor condensation and partial filling of the wedge, formed by the rear wall of the ridge and the slit wall, take place (see also 7a). Since fluid imbibition in a capillary is accompanied by the faster motion of a precursor film far ahead of the meniscus [49, 50, 48], it is also conceivable that the wedge is filled by this atomically thin precursor. Since the meniscus position at time $t$ is difficult to locate precisely, we rather measure the volume of the fluid in the capillary which is readily obtained by the total number of fluid atoms, residing at time $t$ in the slit during the filling process. For an incompressible fluid, this volume is directly proportional to the distance $z(t)$, travelled by the fluid front at time $t$. The curve giving $z(t)$ is shown in  7b for several values of the contact angle, $30^{\circ}\leq\theta\leq 64^{\circ}$. The shape changes also entail a substantial change of the front speed, as well documented in 5b, 6b and 7b , from which the front speed can be read off simply as the slope of the curves reporting the front position as a function of time. These figures show evidence of a significant acceleration of the front (centerline position $y=H/2$) in the climbing stage, followed by a deceleration in the stage where the front is approaching the rear corner, where the chance/risk of pinning is highest. Here the fate of the front becomes critical. According to the CF analysis, if the contact angle is below $45^{\circ}$, there is enough drive from the upper wall to pull the front away from the rear edge of the obstacle. Otherwise, the front stops moving.

The CFD and LB simulations confirm this picture, showing evidence of pinning only for the two angles above $45^{\circ}$. More precisely, they show a significant front acceleration in the climbing stage, followed by a coasting period, once the rear edge is reached. For $\theta=30^{\circ}$ and $\theta=40^{\circ}$, the coasting period exhibits a finite lifetime, after which the front regains its motion. For $\theta=50^{\circ}$ and $\theta=60^{\circ}$, the coasting period does not seem to come to an end (within the simulation time), and the front is pinned. After the obstacle, the (unpinned) front is slowed down, as one can appreciate from the slightly reduced slope of the front dynamics, see 5b (only for $\theta=40$) and 6b.

The MD simulations, on the other hand, tell a different story, in that even at $\theta=50^{\circ}$, the front proves capable of overcoming the obstacle, although with a drastically reduced speed. Note that, while both CFD and LB show evidence of a strong front acceleration as they approach the obstacle, MD simply shows a monotonic deceleration. This is due to the fact that while in CFD and LB simulations we measure the distance travelled by the interface midpoint, as already pointed out, MD measures the total volume of the fluid in the capillary.

Given the qualitatively different outcome of CFD(LB) versus MD simulations at $\theta>45^{\circ}$, we next discuss the dependence of the filling dynamics on the different contact angles, case by case.

• •

  Contact angles $30^{\circ}$ and $40^{\circ}$
  As already mentioned, in these cases, all three methods give the same qualitative outcome: in proximity of the ridge, the front deforms in response to the geometrical discontinuity, climbs up the obstacle and walks over its top. Once the rear-edge is reached, the bottom end of the front pins at the corner, and keeps moving on the top wall. Then, according to the CF criterion, the front overcomes the obstacle and completely fills up the entire channel. Manifestly, the standard $z$ versus $t$ relationship, described by the LW-law (1) is strongly violated in the vicinity of the obstacle. After overcoming the obstacle, the usual LW regime is recovered, although after a transitional period of time, which depends on the wettability $\theta$ [17], and with a reduced velocity. In particular, in all cases, the front is slowed down right after passing the obstacle, see 5b, 6b and 7b. However, at $\theta=30$, the asymptotic behaviour is quite different for the three methods: (i) In CFD simulations, 5b, after a short transient time, the velocity basically regains the initial value before passing the obstacle. (ii) LB shows a similar behaviour, with only a slight reduction of the front velocity after passing the obstacle. (iii) in MD, however, the front undergoes a substantial velocity decrease.

  This might be due to the details of the fluid-solid interaction during the obstacle surmounting (in which the contact angle is bound to undergo drastic changes), as well as to the different time-scales. A detailed analysis of this effect shall be deferred to a future study.

• •

  Contact angle $50^{\circ}$
  In this case, while results from CFD and LB confirm the CF criterion, MD simulations show a different behaviour for the front motion: interestingly, the Concus-Finn (Gibbs) criterion for contact line pinning at the edge of the ridge is found to break down. Indeed, for $\theta=50^{\circ}$, the fluid front overcomes the obstacle in manifest violation of the Concus-Finn criterion. In order to vividly explain this feature, we show in 8 the snapshots and in video 1 and 2 the movie of the front dynamics, respectively for the cases $\theta=40^{\circ},~50^{\circ}$. These observations suggest that, on the nanoscale, the overcoming of topological obstacles is strongly affected by interface fluctuations, thus undermining the deterministic nature of the imbibition process. Indeed, both the CFD and LB method work in absence of fluctuations, and this would explain the difference with MD simulations. The problem of contact line pinning during capillary imbibition acquires thus a stochastic character and is most probably governed by the size of the obstacles around the CF critical point. In fact, depending on the height of the ridge obstacle, a coalescence of the pinned meniscus with the molecules ahead of the obstacle, in the vicinity of the edge, may occur at later times.

• •

  Contact angle $60^{\circ}$
  Again, in this case all the three methods give the same result: the front deforms, climbs the obstacle, walks on its top, but pins at the back edge and definitely stops moving.

  

  4.3 Further discussion

  In order to gain a better understanding of the previous results, and particularly of the violation of the CF criterion for mildly super-critical angles in MD simulations, additional runs have been performed. More precisely, we have run MD simulations at $\theta=50^{\circ}$ with a taller obstacle, $h=15\sigma$, in order to inspect whether the front is still capable of surmounting the obstacle in violation of the CF criterion (in this case, the total slit height was correspondingly increased, so as to leave the same clearance above the obstacle as in the previous simulations). Moreover, we have also run additional CFD and MD simulations at $\theta=50^{\circ}$ with shorter obstacles, $h=H/4$, in order to inspect whether even CFD(LB) would show violations of the CF criterion upon reducing the obstacle size. The main outcome is as follows: MD simulations with $h=15\sigma$ do show front pinning, indicating that violations of the CF criterion disappear once the obstacle is made sufficiently tall (see Fig. 9). This corroborates our previous conjecture of a (Arrehnius-like?) dependence of the CF criterion on the obstacle height, in the presence of molecular fluctuations. At the same time, CFD and LB simulations with $h=H/4$ at $\theta=50^{\circ}$ keep showing evidence of pinning, thereby lending further support to the idea that the CF criterion remains insensitive to obstacle size, so long as molecular fluctuations can be neglected.

  

  

  Before concluding, a few words of caution should be spent on the fact that, being diffuse-interface methods, both CFD and LB propagate a finite-width interface. As a result, one may wonder whether and to what extent finite-width effects could interfere with the physics of post surmounting. Indeed, as shown in previous studies by these and other authors ([45, 13]), finite interface width eventually leads to mild deviations from the LW law. However, they do not affect the qualitative outcome of a pinning/no-pinning test, unless the interface width becomes comparable to the characteristic obstacle width. Both CFD and LB simulations are visibly far from this critical limit, which is why we are confident that the qualitative conclusions of the present work are not significantly affected by finite-width effects. One may wonder whether such deviations might be paralleled to the effect of molecular fluctuations. To this regard, it is worth underlining that such a parallel has to be taken very cautiously, since molecular fluctuations stem from the microscopic physics of the problem, while finite-width effects are due to a lack of numerical resolution. Sometimes a mapping between the two can be established, but this can by no means be taken as a general rule. Finally, another potential source of discrepancy between CFD(LB) and MD is the fact that the latter allows slip motion while CFD and LB (in this set up) do not. Although any solid statement on the effect of slip flow on the dynamics of post surmounting must necessarily be deferred to a detailed quantitative analysis, we observe that due to the weakness of inertial effects, hydrodynamic boundary conditions have little or no effect, on the dynamics/energetics of the obstacle surmounting.

  This is confirmed by the fact that slip flow effects should manifest through visible deviations from the LW regime, of which we have no evidence, at least in the parameter regime investigated in this work. This situation can drastically change in the presence of superhydrophobic effects, although we shall not be concerned with these problems in the present work.