Spontaneous imbibition of a liquid film wetting a wall-mounted cylinder corner

Abstract

Spontaneous imbibition flows within confined geometries are commonly encountered in both natural phenomena and industrial applications. A profound knowledge of the underlying flow dynamics benefits a broad spectrum of engineering practices. Nonetheless, within this area, especially concerning complex geometries, there exists a substantial research gap. This work centers on the cylinder-plane geometry, employing a combined theoretical and numerical approach to investigate the process of a wetting film wrapping a cylinder corner. It is found that the advance of the liquid front generally follows the Lucas–Washburn kinetics, i.e., t^1/2 scaling, but it also depends on the dynamics of the liquid source. Furthermore, a theoretical estimation of the timescale associated with the imbibition process is also provided, especially the merging time as an important time length characterizing the duration of the wetting process. The timescale is highly dependent on the wettability conditions and the properties of the involved liquid. The conclusion of this work lays a theoretical foundation for comprehensively understanding the capillary phenomena in complex media and shedding light on various microfluidic applications.

---

1 Introduction

Spontaneous imbibition flows, i.e., liquids driven by capillary pressure to wet confined geometries, such as capillary tubes,¹ grooves,^2,3 porous media,^4,5etc., play a crucial role in various natural and industrial processes, such as designing microfluidic devices,^6–8 water harvesting^9,10 and surface coating and cleaning.^11,12 The pioneering research dates back to the Lucas–Washburn equation,¹³ which leads to a scaling law regarding the evolution of the liquid front h(t) in capillaries, i.e., h = Ct^1/2. This type of scaling law describes an energy balance between the capillary and viscous terms. Specifically, the wetting liquid is driven by the capillary force to spread on the surface and tends to maximize the coverage over the surface, during which the interfacial energy decreases and is consumed by viscous friction. Due to gravity, part of the released interfacial energy is transformed to the gravity potential leading to a different t^1/3 scaling. Nevertheless, the Lucas–Washburn equation was developed for circular capillaries. Once the confined geometry is complex, to what extent the t^1/2 scaling can predict the imbibition dynamics and how to estimate the scaling coefficient C remain unexplored, especially when sharp corners come up in a geometry.

There have been certain works investigating the imbibition in corners. For an open V-shaped groove, Tang and Tang² theoretically proved that the imbibition dynamics follow the t^1/2 scaling ignoring gravity or t^1/3 scaling considering gravity; Higuera et al.¹⁴ derived the same scaling law within the framework of the lubrication approximation. These scaling laws have been verified against the experimental observations.^3,14,15 Very recently, Zhou and Doi¹⁶ developed a theory model for corners with curved walls using the Onsager principle. Surprisingly, they found that the above scaling law still works while the scaling coefficient C slightly depends on the wall shape. In a closed medium, like a square or rectangular tube, if the contact angle θ < 45°, i.e., the Concus–Finn condition is satisfied,¹⁷ and the liquid can wet the interior corners and forms “finger-like” films along the corners. The imbibition flows thus become manifold, i.e., the bulk flow and corner flows, and the synergistic effect of the corner and bulk flows should be carefully considered.¹⁸ Imbibition in square tubes has been numerically and theoretically investigated.^19–21 It is found that both flows follow the Lucas–Washburn kinetics and their coupling plays an evident role.

What's more complex, in a porous medium, especially a natural one, the inter-connected angular channels randomly distribute in its solid space, where corner flows are enhanced and bulk–corner flows are expected to interplay in a more complicated manner. Cylinder-based geometries are commonly used as a surrogate model of real porous media. A number of experimental and numerical works have been reported based on these geometrical settings. Zhao et al.²² conducted a microfluidic experiment and directly visualized the process of liquid film spreading among cylinder corners in the strong imbibition regime (θ < 30°). Numerical modelling works on corner flows in cylinder-based porous media are ensued,^23–26 and the corner flow is regarded as a specific fluid–fluid displacement pattern and emerges under certain combining conditions of capillary number, viscosity ratio and wettability.

Though great progress regarding spontaneous imbibition flows within complex geometries has been made, answers to the fundamental questions posed in the beginning are still demanded because they are steps towards better engineering practices. In this work, I shall focus on the cylinder-plane geometry and investigate the process of a liquid film wrapping a cylinder theoretically and numerically.

2 Theoretical model

I consider a wall-mounted cylinder wetted by a liquid film, as shown in Fig. 1(a). In this setting, the wetting film symmetrically spreads along the cylinder-bottom corner from a liquid source and finally merges at the other end. For describing this problem, a cylindrical coordinate (r–φ–z) is set up, where the liquid source locates at φ = 0 while the liquid front locates at φ = φ_m, as can be seen in Fig. 1(b). Here, the film wetting is along a horizontal corner rather than against the gravity direction, and I assume that the characteristic size of the liquid film is smaller than the capillary length , where ρg is the liquid gravity and γ is the surface tension. Thus, the effect of gravity can be neglected. Additionally, I assume that the liquid–gas interface on the z–r plane is an arc, as shown in Fig. 1(c). Thus, the wetting height h_w and width r_w are equal. Furthermore, since the moving velocity of the contact line remains low, I assume a constant contact angle during the whole process rather than a dynamic one.²⁷ With the contact angle θ and wetting width r_w, the film thickness h as a function of r is expressed as:where R₀ is the cylinder radius and R = r_w/(cos θ − sin θ).

Fig. 1 Schematics of the theoretical model, including a perspective view (a), a top view (b) and a sectional profile on the r−z plane of a cylindrical coordinate (c).

2.1 Time evolution equation

Using the Onsager principle, I derive the time evolution equation for the meniscus profile, which can be characterized by r_w(φ,t) for a given θ as per eqn (1). For the present problem, it is stated in this principle that the dynamics of the system can be directly determined by the minimum of the Rayleighian,²⁸

[scr R, script letter R][ṙ_w(φ,t)] = Ḟ[ṙ_w(φ,t)] + Φ[ṙ_w(φ,t)], (2)

where Ḟ is the change rate of the free energy of the film-cylinder system and Φ is the energy dissipation function.

2.1.1 The change rate of free energy. The free energy of the system is a superposition of the interfacial energies along the liquid-cylinder wetting area A_ls1(r_w), the liquid-wall wetting area A_ls2(r_w) and the liquid–gas area A_lg(r_w), and is given by

F = γ(−A_ls1 cos θ − A_ls2 cos θ + A_lg), (3)

whereand h_r and h_φ are the derivatives of h concerning r and φ, respectively. The change rate of the free energy Ḟ is thus obtained as:Here, the top dot denotes the time derivative and the prime denotes the derivative with respect to r_w. Separately, and can be directly derived as:For , since the film is much thinner than the cylinder radius, i.e., h ≪ R₀, and moreover, h_φ²/r² ≪ h_r² ≪ 1, can be given as a simplified form:An auxiliary variable a′(r_w) is defined as an integrated part of F for the convenience of following usages, i.e.,I take the volume flux Q(φ,t) of the liquid flowing across the cross-section area, as shown in Fig. 1 at φ, as an independent variable. Here, Q(φ,t) is related to ṙ_w(φ,t) by the conservation equation, which readswhere is the cross-sectional area. Using the conservation equation eqn (12), I can rewrite the change rate of free energy as a function of Q instead of ṙ_w,The definition of Ḟ as per eqn (13) suggests that Ḟ is a measurement of the power of the capillary force. Thus, the capillary pressure P_c of the film-cylinder system can be estimated as

2.1.2 Dissipation function. Assuming that the liquid imbibes slowly along a cylinder corner, the inertia effect can be neglected. The flow is almost one-dimensional since u_φ is much larger than u_r and u_z. Thus, flow dynamics can be described by the following Stokes equation:where ∂P/(R₀∂φ) is the pressure gradient along the φ-axis. Provided ∂P/(R₀∂φ), eqn (15) is solved on the domain shown in Fig. 1 with no-slip boundary conditions, i.e., u_φ = 0 at the solid walls and shear-free boundary conditions, i.e., n·∇u_φ = 0 at the gas–liquid interface, where n is the normal vector of the interface within the r−z plane. The volume flux,and according to Darcy's law,where k is the permeability of the planar meniscus with units of m². It is determined by the characteristic length of the meniscus, naturally taking r_w. Thus, k is in the form of

k = r_w²k̄(θ). (18)

Here, k̄(θ), as a function of wettability, describes the effect of the meniscus shape and is obtained numerically, see Appendix A for details. The dissipation function is then expressed as:Ḟ and Φ are expressed with respect to Q, and the Rayleighian is given as:The governing equation is derived from the Onsager variational principle, δ[scr R, script letter R]/δQ = 0,Using the conservation equation eqn (12) again, I express the governing equation concerning r_w,Substituting h, A_l and k in eqn (22), a dimensionless form of the governing equation is obtained,Its length is scaled by R₀ and its time is scaled by a characteristic time t*

2.2 Theoretical analysis

The time evolution equation eqn (23) suggests a scaling relationship,and thus it admits a self-similar solution in the form ofwhere H(·) is a function to be determined. Substituting eqn (26) in eqn (23), it gives an ordinary differential equation,

2HH′ + (8H′² + 4H′′H + H′)χ = 0, (27)

where the prime represents the derivative regarding χ. When χ = 0, it corresponds to the boundary condition at the liquid source (φ = 0), i.e., H(0) = r_w|_φ=0 > 0, and from eqn (27) it leads to:

H′(0) = 0 (28)

Another boundary condition is at the liquid front where H(χ) approaches zero at a certain value χ = χ₀, i.e.,

H(χ₀) = 0. (29)

Substituting eqn (29) in eqn (27), I obtainTo satisfy eqn (29) and (30), H(χ) is assumed to be in the form of:where parameters n_i and a_i are to be determined. According to eqn (28), I obtain

I consider a situation with a fixed r_w at the liquid source (φ = 0) i.e., r_w|_φ=0 = r⁰_w, and it leads to

Anticipating n_i > 1, the upper and lower bounds of are determined from eqn (32),where n^max_i and n^min_i are the maximum and minimum values of n_i. Thus, can be estimated as:where n^max_i ≤ n̄ ≤ n^max_i. Furthermore, substituting it into eqn (33), an asymptotic solution of the liquid front φ_m is obtained,This suggests φ_m ∼ t^1/2 which aligns with the liquid imbibition in a capillary tube or homogeneous porous media described by the Lucas–Washburn equation.²⁹ Furthermore, the merging time t_merge, at which two liquid fronts from both sides touch each other, can be estimated. Here, I only consider the contribution of the linear term in eqn (31), and by letting φ_m = π,

2.3 Numerical solution

I now numerically solve the time-evolution equation eqn (23) for validating my proposed law φ_m ∼ t^1/2. Besides the boundary conditions at the liquid source (φ = 0), the one at the merging point (φ = π) is set as r_w|_φ=π = r^min_w. Then, the capillary pressure is calculated as per eqn (14),Since the liquid front is regarded as a point, r^min_w should be zero. However, the capillary pressure would be an infinite value if r^min_w = 0 as per eqn (38), resulting in a convergence issue. Therefore, I take a finitely small value as r^min_w, and r_w is initialized with r^min_w, i.e., r_w|_t=0 = r^min_w. eqn (23) with the boundary conditions is solved on a domain φ ∈ [0, π] using the finite element method.

I first investigate the effect of r^min_w. As shown in Fig. 2, cases with r^min_w ranging from 3 × 10⁻⁵ to 1 × 10⁻³ are almost overlapped regarding the time evolution of the liquid front position in Fig. 2(a) and the r_w profiles in Fig. 2(b). A difference is observed in the zoomed-in plot around the liquid front in Fig. 2(b), suggesting that the value of r^min_w only influences the local region in the vicinity of the liquid front. More importantly, the measured log–log slope of curves φ_mvs. t, as shown in Fig. 2(a), confirms φ_m ∼ t^1/2 at the later stage times.

Fig. 2 (a) The evolution of φ_m for the group of cases with r^w₀ = 0.03 and various r^w_min ∈ [1 × 10⁻³, 3 × 10⁻⁴, 1 × 10⁻⁴, 3 × 10⁻⁵]. (b) The corresponding r_w profiles at different times which are marked by black triangles in (a), and the inset is a zoomed-in plot of liquid fronts.

Another scaling law that , suggested by eqn (36), is rationalized and verified. From eqn (38), it suggests that the larger r⁰_w is, the stronger P_c is and thus the faster the wetting film spreads along the corner. Furthermore, as shown in Fig. 3(a), cases with various r⁰_w ranging from 0.03 to 0.12 collapsed as one line on the space. In addition, the merging time t_merge for each case is directly measured from the numerical result and compared against the theoretical estimation from eqn (37). Fig. 3(b) shows that both numerical solutions and theoretical estimations have the same trend, but eqn (37) underestimates t_merge as per the comparison. This inconsistency should be attributed to the transition period at the early time, as can be seen in Fig. 3(a). During the transition period, the interfacial profile is relaxed and self-adjusted to progressively follow the law φ_m ∼ t^1/2. Nevertheless, predicting the transition period is out of the scope of the theoretical model.

Fig. 3 (a) The evolution of scaled φ_m for the group of cases with r^w_min = 3 × 10⁻⁴ and various r^w₀ ∈ [0.03, 0.06, 0.09, 0.12]. (b) The comparison of t_merge obtained from the numerical solution and theoretical estimation (eqn (37)) under various r^w₀.

3 Volume-of-Fluid simulation

Given that my theoretical model is developed on the foundational assumption of the “arc-shape interface”, it is necessary to gauge the practical applicability of my theoretical model and further test the proposed scaling law. In this section, I will conduct numerical simulations using the volume-of-fluid (VoF) method. Not only for the verification, I also investigate the film wrapping problems under diverse conditions.

3.1 Governing equations

I consider the imbibition as a laminar, incompressible, and immiscible two-phase flow, which is governed by the Navier–Stokes equations,

∇·v = 0, (39)

ρ∂v/∂t + ρ∇·(vv) = −∇p + μ∇²v + F_γ, (40)

where v denotes the velocity vector; p, ρ, and μ are respectively the fluid pressure, density and viscosity; F_γ is the surface tension force per unit volume. The interface between two phases is tracked by the volume-of-fluid (VoF) method, wherein a scalar transport equation regarding the volume fraction α is introduced,

∂α/∂t + ∇·(vα) = 0. (41)

The interface is reconstructed based on α-field and related geometric features including interface normal n_α and curvature κ are obtained. Then, F_γ is calculated as³⁰

F_γ = γκ∇α. (42)

Wetting conditions are implemented by correcting n_α in the vicinity of the solid walls,³¹

n_α = n_s cos θ + t_s sin θ, (43)

where n_s and t_s are the unit normal and tangent vectors to solid walls, respectively, see Patel et al.³² for more details on determining t_s. eqn (39)–(41) with the following boundary conditions are solved using OpenFOAM.^33,34

3.2 Numerical model

I build up a three-dimensional numerical model, as shown in Fig. 4(a). Considering this problem is a symmetric one, and a half-cylinder zone is adopted as the computation domain. The symmetry plane, as marked by dash-dot lines in Fig. 4(b), is divided by the cylinder wall into two face boundaries, i.e., the left and right faces. At the right face, where the liquid fronts from both sides will touch, symmetric boundary conditions are imposed for the flow field and the α field. At the left face, I control the α field to simulate different types of the liquid sources, including the “fixed boundary” mimicking the situation where r_w is fixed at the liquid source and the “free boundary” where r_w can freely grow at the liquid source as described in detail in the following. Wetting wall boundary conditions are set on the cylinder wall and the bottom wall, as marked in Fig. 4(a), following eqn (43). Other boundaries connect to the environment and thus a zero-pressure condition and a zero-gradient α field are imposed.

Fig. 4 Geometrical settings of the numerical model in a perspective (a), front (b), and top view (c).

The radius of the domain is 3R₀ and its height is 2R₀. The upper limit of the mesh size is set as R₀/100, which has passed the mesh-sensitive test. I set the viscosity ratio as 100 which is large enough to represent a gas–liquid situation. The quantities including r_w, h_w, and φ_m are directly measured from the reconstructed interface. For the convenience of comparing with the theoretical model, all lengths and times presented in the following have been scaled by R₀ and t*, respectively.

3.2.1 Fixed boundary. I first simulate the situation with the fixed α field at the left face, which is expected to agree with the theoretical predictions in Section 2.2. Specifically, provided r⁰_w and θ, the interface position at the left face is calculated as per eqn (1), and then the liquid and gas phases separated by the interface are mapped on the α field at the left face.

I conduct simulations over a range of r⁰_w ∈ [0.3,0.5] and θ ∈ [15°, 20°, 25°, 30°]. Fig. 5(a) shows the evolution of φ_m scaled by in the log–log space. For the group of cases with the same r⁰_w, simulation results from various-θ cases are overlapped, suggesting that the effect of wettability is well considered in t*. Moreover, the scaling law is also verified to a good extent, since the two groups are significantly close to each other and almost collapse as one line, though a small gap is observed. To better gain insights into the evolution of φ_m, I calculate the secant slopes of φ_m–t curves in the log–log space, as defined as

where dt is the scaled time interval. As shown in Fig. 5(b), each case has a transition period at the early time, during which the secant slope sharply decreases from a large value and then becomes flattened. The length of such a transition period depends on r⁰_w and θ, but it generally takes around 2t* before the evolution reaches the steady state. The steady slope, though floating over a range of [0.47,0.55], is close to 0.5, indicating that the proposed law φ_m ∼ t^1/2 effectively governs the imbibition dynamics. Besides, t_merge measured from simulation results is compared with the theoretical estimation from eqn (37), as presented in Fig. 5(c). The scaled t_merge seems a function of θ, while it should be independent of θ according to the theoretical model where the impacts of θ have been considered in t*. This is owing to the transition period which is θ-dependent and involved in the measured t_merge. Although deviations between predicted and measured t_merge are observed, the theoretical model provides a reasonable lower-bound estimation of t_merge.

Fig. 5 Simulation results of the fixed-boundary situation with r⁰_w ∈ [0.3,0.5] and θ ∈ [15°, 20°, 25°, 30°]. The evolution of (a) the scaled φ_m and (b) the corresponding secant slope. (c) The comparison of t_merge against the theoretical prediction. (d) The wetting height h_wvs. the wetting width r_w at φ = π/2.

What's more, to further confirm the practicability of my theoretical model, I test the foundational assumption that the interface on the r–z plane maintains an arc-shape. Fig. 5(d) shows the evolution of h_w − r_w at φ = π/2 of each case. With imbibition ongoing, the wetting film expands within the r–z plane and h_w should increase at the same rate with r_w as per the assumption, i.e., r_w = h_w as marked by the dashed line in Fig. 5(d). It is observed that the measured r_w − h_w aligns well with the assumption, especially at the early time when r_w is small. With the increasing r_w, though a slight deviation occurs, i.e., h_w becomes smaller than r_w, the assumption is still acceptable. Noteworthily, this deviation is only determined by the relative size of the wetting film to the cylinder radius. In my theoretical model, only the curvature within the r–z plane is considered for calculating the capillary pressure. However, with the wetting film expanding and r_w increasing to close to 1, the contribution of the other principle curvature to the capillary pressure may not be neglected. Thus, the effective scope of my theoretical model should be limited to the “small-film-size” regime. Additionally, the deviation from the “arc-shaped interface” assumption could be another source of the failure in precisely predicting t_merge.

3.2.2 Free boundary. I then extend my focus to another situation where the size of the wetting film at the liquid source can freely grow. Correspondingly, the zero-gradient boundary condition for the α field is imposed at the left face.

The simulation cases cover various θ ∈ [15°, 20°, 25°, 30°]. Initially, a small arc-shaped patch (around 0.05R₀) is set as a liquid phase at the corner of the left face. It relaxes and evolves to form a meniscus after one recording time step dt. I regard the size of such formed meniscus as an initial value r⁰_w at the liquid source, which depends on θ, as shown in Fig. 6(a). However, since growth curves under various θ are observed parallel, the growths of r_w at the liquid source are in a similar track, approximately following a power law. The average power is measured as 0.23, which is marked in Fig. 6(b). Equivalently, for the theoretical model, the boundary at the liquid source r_w|_φ=0 is time-dependent, i.e.,

r_w|_φ=0 ≈ r⁰_wt^0.23. (45)

The analysis in Section 2.2 remains effective but an adaption is needed. Considering the transient formation of r_w|_φ=0, eqn (36) is modified asThus, I obtain an approximate scaling law φ_m ∼ t^0.615 adapted to the free-boundary situation. The liquid front position is measured from my simulation results, and its evolution and secant slopes are demonstrated in Fig. 6(c) and (d). Similarly, after a transition period, liquid front advancing reaches a steady state. The steady slope of each case tends to be around 0.6, as marked in Fig. 6(d), which is comparable to the theoretically predicted value of 0.615. Moreover, based on eqn (46), I can estimate t_merge asThis estimation still serves as a lower bound of t_merge, as observed in Fig. 6(e). Again, I test the foundational assumption of the “arc-shaped interface” in the free-boundary situation using h_w − r_w on the φ = π/2 plane. As shown in Fig. 6(f), the deviation is linearly enlarged with r_w, and the relative error (r_w − h_w)/r_w is larger than 10% when r_w = 0.8, probably suggesting that the contribution of the secondary principle curvature has to be considered if r_w further increases.

Fig. 6 Simulation results of the free-boundary situation with θ ∈ [15°, 20°, 25°, 30°]. The evolution of (a) r_w at the liquid source and (b) the corresponding secant slope. The evolution of (c) φ_m and (d) the corresponding secant slope. (e) The comparison of t_merge against the theoretical prediction. (f) The wetting height h_wvs. the wetting width r_w at φ = π/2.

I now shift my focus to imbibition dynamics after merging. Though post-merging behaviours are beyond the scope of the theoretical model, my simulation results provide insights into them. After the two fronts merge at the right face, the film continues to expand under the free-boundary conditions. I show the evolution of r_w at the right face in Fig. 7(a) and the secant slopes in Fig. 7(b). The expanding rate of r_w decreases at the beginning and gradually tends to be a constant value, i.e., 1.11 as marked in Fig. 7(b). In another word, r_w increases with time approximately in a linear mode, which is significantly faster compared to the one at the liquid source, see Fig. 6(b).

Fig. 7 The post-merging dynamics of the free-boundary situation including (a) the evolution of r_w at the right face and (b) the corresponding secant slope.

4 Conclusions

In this work, I have theoretically and numerically investigated the spontaneous imbibition of a liquid wetting a cylinder corner. Using the Onsager variational principle, a time evolution equation for the meniscus profile was derived. Based on the time evolution equation, I derived an asymptotic solution of the liquid front . This suggests that the advance of the liquid front follows the Lucas–Washburn kinetics, i.e., the t^1/2 scaling, if the boundary r⁰_w is time-independent; otherwise, the effect of the dynamic boundary should be included and the scaling accordingly changes. Then, the imbibition process was numerically simulated using the VoF method, and the simulation results can be well rationalized by my proposed scaling law to a large extent. Furthermore, I provide a theoretical prediction of t_merge, which is demonstrated as a lower bound of the real one.

My theoretical model is extensible. More complex geometries, such as tapered, ellipse, or even any arbitrary-shape symmetric cylinders, can be modelled by modifying the expression of free energy. I can expect that the scaling coefficient C and characteristic time t* vary with the geometry while the scaling t^1/2 remains effective. Moreover, another demanding aspect for future works is to investigate the imbibition flows in a cylinder group, and model how the liquid front spreads among neighboring cylinders.

The developed framework is expected to provide theoretical support for deep understanding of the capillary-driven flows in complex media, e.g., fibre-based fabrics and artificial functional surfaces. The obtained conclusions build up the scientific foundation for optimal designs of microfluidic devices.

Conflicts of interest

There are no conflicts to declare.

Appendix A. Determination of k̄(θ)

I determine the relative permeability k̄(θ) using numerical experiments. eqn (15) is solved on an axisymmetric meniscus domain, as shown in Fig. 8, whose geometry is dependent on θ and r_w. I sweep the parameter combinations of θ ∈ [15°, 20°, 25°, 30°] and r_w ∈ [0.10, 0.15, 0.20, 0.25, 0.30], and calculate the permeability k according to eqn (17). Fig. 9 shows that the permeability k is proportional to r_w² for any θ. Thus, the relative permeability k̄(θ) can be obtained by measuring the slope of k – r_w² lines, which are summarized in Table 1.

Fig. 8 The computation model for determining the relative permeability k̄(θ).

Fig. 9 Permeability k vs. r_w² for various θ.

Table 1 Relative permeability k̄(θ)

15°	20°	25°	30°
0.01772	0.02032	0.02305	0.02591

---

Acknowledgements

This research was undertaken with the financial support of the Swedish Foundation for Strategic Research under SSF-FFL15-001. Access to the supercomputing resources of the PDC Center for High Performance Computing and National Supercomputer Centre (NSC) used for this work were provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS). I acknowledge the early discussion with Prof. Shervin Bagheri, the casual discussion with Mr Boxue Zheng at the dinner table, and suggestions from Dr Zhaojun Zhang.

References

1. J. Cai, T. Jin, J. Kou, S. Zou, J. Xiao and Q. Meng, Langmuir, 2021, 37, 1623–1636.
2. L.-H. Tang and Y. Tang, J. Phys. II, 1994, 4, 881–890.
3. D. Deng, Y. Tang, J. Zeng, S. Yang and H. Shao, Int. J. Heat Mass Transfer, 2014, 77, 311–320.
4. S. Suo, M. Liu and Y. Gan, Transp. Porous Media, 2019, 126, 615–631.
5. J. Ha, J. Kim, Y. Jung, G. Yun, D.-N. Kim and H.-Y. Kim, Sci. Adv., 2018, 4, eaao7051.
6. Y.-H. Lai, J.-T. Yang and D.-B. Shieh, Lab Chip, 2010, 10, 499–504.
7. M. A. Burns, B. N. Johnson, S. N. Brahmasandra, K. Handique, J. R. Webster, M. Krishnan, T. S. Sammarco, P. M. Man, D. Jones and D. Heldsinger, et al. , Science, 1998, 282, 484–487.
8. S. Suo and Y. Gan, Langmuir, 2020, 36, 13642–13648.
9. J. Li and Z. Guo, Nanoscale, 2018, 10, 13814–13831.
10. H. Chen, P. Zhang, L. Zhang, H. Liu, Y. Jiang, D. Zhang, Z. Han and L. Jiang, Nature, 2016, 532, 85–89.
11. M. Heinz, P. Stephan and T. Gambaryan-Roisman, Colloids Surf., A, 2021, 631, 127450.
12. S. Suo, H. Zhao, S. Bagheri, P. Yu and Y. Gan, Chem. Eng. Sci., 2022, 264, 118134.
13. E. W. Washburn, Phys. Rev., 1921, 17, 273.
14. F. Higuera, A. Medina and A. Linan, Phys. Fluids, 2008, 20, 102102.
15. R. Rye, F. Yost and J. Mann, Langmuir, 1996, 12, 4625–4627.
16. J. Zhou and M. Doi, J. Fluid Mech., 2020, 900, A29.
17. P. Concus and R. Finn, Proc. Natl. Acad. Sci. U. S. A., 1969, 63, 292–299.
18. M. M. Weislogel, J. Fluid Mech., 2012, 709, 622–647.
19. T. Yu, J. Zhou and M. Doi, Soft Matter, 2018, 14, 9263–9270.
20. J. Zhao, F. Qin, R. Fischer, Q. Kang, D. Derome and J. Carmeliet, Water Resour. Res., 2021, 57, e2020WR029190.
21. V. T. Gurumurthy, D. Rettenmaier, I. V. Roisman, C. Tropea and S. Garoff, Colloids Surf., A, 2018, 544, 118–126.
22. B. Zhao, C. W. MacMinn and R. Juanes, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 10251–10256.
23. B. K. Primkulov, A. A. Pahlavan, X. Fu, B. Zhao, C. W. MacMinn and R. Juanes, J. Fluid Mech., 2021, 923, A34.
24. B. Zhao, C. W. MacMinn, B. K. Primkulov, Y. Chen, A. J. Valocchi, J. Zhao, Q. Kang, K. Bruning, J. E. McClure and C. T. Miller, et al. , Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 13799–13806.
25. S. Cox, A. Davarpanah and W. Rossen, Transp. Porous Media, 2023, 147, 197–216.
26. R. Hu, J. Wan, Z. Yang, Y.-F. Chen and T. Tokunaga, Geophys. Res. Lett., 2018, 45, 3077–3086.
27. Y. Zhang, M. Guo, D. Seveno and J. De Coninck, Adv. Colloid Interface Sci., 2023, 313, 102861.
28. M. Doi, Soft Matter Physics, Oxford University Press, USA, 2013.
29. J. Cai, Y. Chen, Y. Liu, S. Li and C. Sun, Adv. Colloid Interface Sci., 2022, 304, 102654.
30. J. U. Brackbill, D. B. Kothe and C. Zemach, J. Comput. Phys., 1992, 100, 335–354.
31. A. A. Saha and S. K. Mitra, J. Colloid Interface Sci., 2009, 339, 461–480.
32. H. Patel, S. Das, J. Kuipers, J. Padding and E. Peters, Chem. Eng. Sci., 2017, 166, 28–41.
33. J. Roenby, H. Bredmose and H. Jasak, R. Soc. Open Sci., 2016, 3, 160405.
34. H. Scheufler and J. Roenby, J. Comput. Phys., 2019, 383, 1–23.

---