Computations of spontaneous rise of a rivulet in a corner of a vertical square capillary

Accepted Manuscript Title: Computations of Spontaneous Rise of a Rivulet in a Corner of a Vertical Square Capillary Author: Vignesh Thammanna Gurumurthy Daniel Rettenmaier Ilia V. Roisman Cameron Tropea Stephen Garoff PII: DOI: Reference:

S0927-7757(18)30079-7 https://doi.org/doi:10.1016/j.colsurfa.2018.02.003 COLSUA 22255

To appear in:

Colloids and Surfaces A: Physicochem. Eng. Aspects

Received date: Revised date: Accepted date:

29-12-2017 1-2-2018 1-2-2018

Please cite this article as: Vignesh Thammanna Gurumurthy, Daniel Rettenmaier, Ilia V. Roisman, Cameron Tropea, Stephen Garoff, Computations of Spontaneous Rise of a Rivulet in a Corner of a Vertical Square Capillary, (2018), https://doi.org/10.1016/j.colsurfa.2018.02.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Computations of Spontaneous Rise of a Rivulet in a Corner of a Vertical Square Capillary

Institute of Fluid Mechanics and Aerodynamics, Alarich-Weiss-Str. 10 Technische Universit¨at Darmstadt, 64287 Darmstadt, Germany

Stephen Garoff

us

cr

Department of Physics and Center for Complex Fluids Engineering Carnegie Mellon University, Pittsburgh, PA 15213, USA

ip t

Vignesh Thammanna Gurumurthy, Daniel Rettenmaier, Ilia V. Roisman, Cameron Tropea

Abstract

ed

M

an

In this study, the spontaneous rise of a Newtonian liquid in a square capillary with completely or partially wetted walls is investigated using numerical simulations. The flow is modelled using Volume-of-Fluid method with adaptive mesh refinement to resolve the interface for high accuracy. The computations show that for contact angles smaller than 45◦ , rivulets appear in the corners of the capillary. At large times the length of the capillary growth approaches the one-third power function of time. The same asymptotic behaviour has been identified in the existing experimental observations for corners of different geometries. The computations predict the dependence of the rate of the rivulet growth on the liquid viscosity, gravity, width of the capillary and the contact angle. The flow in the rivulet is described using a long-wave approximation which considers three regions of the rivulet flow: the flow near the rivulet tip which is described by a similarity solution, intermediate region approaching a static rivulet shape, and the bulk meniscus. Finally, a scaling analysis is proposed which predicts rivulet growth rates for the given parameters.

1. Introduction

pt

Keywords: capillary flows, capillary rise, interior corner, noncircular capillaries, imbibition

Ac ce

Capillary rise in corners between any two intersecting surfaces and in capillaries and containers with corners can be found in a wide range of applications in the field of microfluidics, fluid management in space and microgravity, porous media, etc. The capillary rise in the corners can be either detrimental, such as in a microchannel, where it alters the flow by lubricating the channel and also by premature mixing of fluids [1], or beneficial, such as its use to fabricate microstructures used for force measurements by capillary micromolding [2, 3]. In the case of porous or wet granular media, where the pores are of different shapes with varying corner geometries, the water imbibition is influenced significantly by the flows in the corners of the pores [4–8]. Understanding the behavior and optimizing the performance in such

Email addresses: [email protected] (Ilia V. Roisman), [email protected] (Stephen Garoff) Preprint submitted to Colloids and Surfaces A

applications requires quantitative modelling of the dynamics of capillary flows in corners. Investigations of the liquid rise by wettable solids started as early as 1710 by Taylor [9] and Hauksbee [10] even before the phenomena of surface tension and wall wettability have been understood [11]. They report the rise of the liquid between two vertical parallel glass plates, in tubes and between the vertical plates forming a small angle. The observed hyperbolic shape of the liquid between two inclined vertical plates has been explained by the theory of Newton, published in his 31st Query [12], predicting that the height of the liquid between two parallel plates is inversely proportional to the interval between them. When the contact angle is in the range π/2 − α < θ < π/2 (where α is the half-angle between the wedge faces) the bulk meniscus shape is nearly spherical and the contact line forms a closed curve on the capillary walls. However, if the contact angle is smaller, θ + α < π/2, the near-spherical meniscus solution is impossible. In February 1, 2018

Page 1 of 13

draulic diameter approximation [24], full computation of the flow in the given cross-section [22, 23, 25, 26] and the Laplacian scaling approach [27]. Further studies [27–29] on capillary flows in the interior corner of containers has been performed under zero and microgravity conditions. Using the lubrication approximation, an asymptotic solution was derived and three regimes were identified: an initial inertial regime h ∼ t followed by a constant flow regime h ∼ t3/5 then an overlap regime and finally constant height regime h ∼ t1/2 . The capillary rise h(t) in the corner under gravity at long times [30] follows the remote asymptotic h ∼ t1/3 , if the contact angle θ < π/2 − α. A similar asymptotic [31] has been obtained for the liquid rise in a narrow wedge. The one-third asymptotic has been recently confirmed in the experiments [32] on the capillary rise in corner formed between different geometries. Numerical simulations can further enhance our understanding of the flow in corners by resolving them spatially and temporally. However, the few studies that have been undertaken only discuss the influence of surface roughness [33] on the early stages of corner rise and the static configurations of meniscus shapes [34]. This study presents three-dimensional computations of the capillary rise under gravity in a vertical square capillary. The influence of viscous and capillary forces, gravity and wall wettability are accounted in the code. The study is focused on the effect of the width of the capillary w, gravity, liquid viscosity µ and of the contact angle θ on the rate of the rivulet propagation along the corners of a vertical square capillary. The computations predict the propagation of the bulk meniscus and the emergence of the corner rivulets if the contact angle satisfies the well-known condition θ < π/4 for the square capillary (with α = π/4). The computations also show that even if the width of the capillary is much smaller than the capillary length of the fluid, the asymptotic behaviour h ∼ t1/3 remains valid. Then, the flow in the rivulet is described by consideration of the mass and momentum balance equations which are satisfied in the integral sense in the framework of the long-wave approximation. The solution for the flow consists of three main regions. A similarity solution is found for the instationary flow near the rivulet tip. This solution is matched with the intermediate region of the rivulet which at very large times approaches the static shape. The third region near the bulk meniscus has a complex three-dimensional form which cannot be easily treated analytically. Therefore the solution for the rivulet flow is not geometrically matched with the meniscus region. The shape of the bulk meniscus can be

Ac ce

pt

ed

M

an

us

cr

ip t

this case rivulets are formed in the corners near the bulk meniscus. They propagate along the wedge edge [13]. For a square capillary the half-angle is α = 45◦ . Therefore the threshold value of the contact angle, corresponding to the formation of the rivulet flow in the corners is θ = 45◦ . Further discussions on the static shape of menisci in the interior corners can be found for triangular tubes [14], polygonal containers [15] and for other complex geometries [3]. Both the bulk and corner rise have been studied extensively, but separately. The bulk rise, i.e., the rise in the central portion of the square capillary is very similar to the one observed in a circular capillary, as the forces driving and opposing the rise are the same. Hence, the model describing the capillary rise in the circular capillary can be applied also to the description of the flow in a square or rectangular capillary with minor modifications accounting for the geometry [16, 17]. It has been shown that the equilibrium height z0 reached by the bulk rise decreases only slightly due to the corner rise [18]. A comparison of several expressions for equilibrium height in square/rectangular capillaries can be found in the literature [19, 20]. The dynamics of liquid rise in an isolated corner under the assumptions of negligible inertia, slender liquid column and zero-gravity, can be modelled using a longwave approximation of the viscous flow in the rivulet. A similarity solution for the rivulet propagation is obtained in [21], which predicts the evolution of the rivulet height h proportional to the square root of time, h ∼ t1/2 , at long times. Such scaling is confirmed by the experimental data [21] in square capillaries placed horizontally, for which the effect of gravity is negligibly small. Theoretical analysis of the rivulet flows in the corner are usually based on the long-wave approximations by assuming a fully developed laminar flow, with negligible inertial effects and also, neglecting the deformation of the interface due to viscous stresses and gravity. One of the difficulties however is an accurate estimation of the viscous stresses in the rivulet cross-section. The corner flow is usually analysed independently on the propagation of the bulk meniscus or on the flow in the neighboring corners. In many studies [22–27], the liquid flow along the corner is computed numerically. The resulting solution for the viscous stresses is presented in the integral form, in terms of a dimensionless friction factor which measures the hydrodynamic resistance to the flow. The influence of the corner geometry (roundedness, corner angle) and contact angle is discussed in these studies. They differ primarily in the modelling approaches used for the friction factor which includes hy-

2

Page 2 of 13

where velocity given by  Ur is the compression ∇γ min Cγ |U|, max(|U|) |∇γ| and Cγ is the compression factor set to 1. The density and viscosity (ρ, µ) of the fluid are calculated from the volume weighted averages of the individual phase as given below

accounted for only using full numerical computations. The dimensionless rivulet profile has been obtained from the computation, using the scaling provided by the theory. This shape agrees well with the analytical solution. Finally, the dimensionless rate of the capillary growth is computed and presented as a function of the dimensionless initial p rivulet thickness, scaled by the capillary length, a ≡ σ/ρg.

φ = γφ1 + (1 − γ)φ2 ,

(5)

2. Numerical methodology For numerical simulations a modified version of the interFoam solver is used, which is implemented in the open source, computational fluid dynamics toolbox OpenFOAM [35, 36]. In interFoam, the two phase flow is modelled using the Volume-of-Fluid method, where the two incompressible, immiscible fluids are treated as one mixed fluid given by an indicator function γ. Here, γ represents the volume fraction of the liquid (dense phase) so γ = 1 represents the liquid phase and γ = 0 the gas phase and any value in between indicates the presence of the liquid-gas interface. Thus the distribution of the mixed fluid in the domain includes both phases as well as the interface between them. Both the liquid and gas phase are treated as incompressible fluids since the flow velocity is much smaller that the speed of sound. The governing equations under incompressible, isothermal conditions for a laminar flow can be written as

Fσ = σκ∇α,

cr

ip t

where, φ = ρ or µ and the suffices 1, 2 represent the liquid and gaseous phase respectively. Based on the Continuum Surface Force (CSF) approach [38], the volumetric surface tension force (Fσ ) acting on the interface can be written as (6)

ed

M

an

us

here κ represents the curvature of the liquid-gas interface given by ∇ · n and the gradient term ∇γ ensures that this force acts only on the interface region. In interFoam ∇γ , But the interface normal (n) is calculated from γ as |∇γ| here the iso-contour approach [39, 40] is used where the contour of γ = 0.5 is used for calculating the normal. At each time step the cell centered γ values are interpolated to its vertices. Then the points (x p ) where the contour cuts the cell are calculated by linear interpolation which is used to reconstruct the interface with a plane in the mixed cells. Finally the interface normal is calculated from the surface vector (S p ) to the plane as S p /|S p |. The surface wettability which governs the contact line dynamics at the solid walls is treated explicitly by imposing a static microscopic contact angle (θ) boundary condition. Numerically, the microscopic contact angle is defined as the angle between the interface normal and the normal to the wall (nw ):

∇·U =0

∂U + U · ∇U = −∇P + ∇ · τ + ρg + Fσ ∂t

(2)

n · nw = cos θ

Ac ce

ρ

pt

(1)

!

Hence, the interface normal must be corrected in order to satisfy (7) before the calculation of curvature at every time step. The above equations (1), (2) and (4) are discretised using a finite volume method on an unstructured grid. The pressure-velocity coupled equations resulting from the discretisation of (1), (2) are solved using the pressure implicit splitting operator algorithm (PISO). The transport equation (4) is solved using multidimensional universal limiter with explicit solution (MULES) algorithm.

where U and P represents the velocity and pressure field respectively andh g the gravity. The i viscous stress tensor T τ is given by µ (∇U) + (∇U) and Fσ represents the volumetric contribution of surface tension force acting on the interface. The distribution of γ is governed by the advection equation which is given below ∂γ + ∇ · (γU) = 0 (3) ∂t In order to avoid interface smearing due to numerical diffusion, a counter-gradient approach [37] is used. According to this method, an additional term is added to equation (3), which artificially compresses the interface, thereby maintaining its sharpness. The modified equation is given as: ∂γ + ∇ · (γU) + ∇ · (γ(1 − γ)Ur ) = 0 ∂t

(7)

2.1. Computational Setup The parameters chosen for the study are the capillary width w, the contact angle θ and the liquid viscosity µ. The widths of the capillaries used are 0.6, 1 , 3 and 4.5 mm and the corresponding Bond number Bo = ρgw2 /σ

(4) 3

Page 3 of 13

ip t cr

us

Figure 2: Meniscus shapes in 3 mm capillary at different contact angles. Images from left to right are at θ = 0◦ , 15◦ , 45◦ and 60◦ at time t = 5 s.

an

Figure 1: Computational domain

PDMS-20 930 19.2 19.8 1.47

PDMS-100 970 97 20.1 1.45

M

Properties Density(kg/m3 ) Viscosity (mPa·s) Surface tension (mN/m) Capillary length (mm)

for the liquid flow and the top face is open to the atmosphere. The front and left faces of the column have the constant contact angle and no-slip boundary condition. This study primarily focuses on the rise of rivulets in the corners and the capillary number based on its speed is very small. Under such conditions, the effect of the dynamic contact angle is minor and therefore a static contact angle at the boundaries can be assumed. Such an assumption is the same as in [27]. In the Volume of Fluid method the simulations on contact line dynamics can be mesh dependent [41–43]. Therefore the dependence of computed results on the mesh refinement has been examined and the difference in the predictions of the rivulet growth when the chosen mesh size is reduced by half is less than 5%. We characterize the capillary rise in the bulk and in the corner by measuring the liquid height in the center and in the corner of the capillary. When describing the corner flow in detail we adopt the convention [44] of measuring the corner rise from the center of the bulk meniscus, which will be described in the section 4.

ed

Table 1: Material properties of the liquids used in the simulations.

Ac ce

pt

based on the capillary sizes are 0.17, 0.46, 4.15 and 9.33 respectively. The contact angles chosen for the study are 0◦ , 15◦ , 30◦ . The physical properties of the liquids used are given in the Table 1 and are set to values for a common fluid used in wettability experiments, polydimethyl siloxane (PDMS). Since the square capillary is symmetric along its cross-section, only one quarter of its cross section is simulated. The computational domain consists of a square column, as shown in the Figure 1, of size w/2 × w/2 × L. The length of the computational domain L is set to 30 mm for all the capillary sizes to ensure that the rising bulk meniscus and the rivulets are located in the computational domain. The base mesh size for each capillary is w/10 and an adaptive mesh refinement is used with load balancing along the liquid-gas interface for higher accuracy. Based on grid convergence tests, the refinement levels are kept at 2, 3, 4 and 4 for capillary sizes 0.6, 1, 3, and 4.5 mm respectively. Figure 1 shows also the typical mesh employed in the simulations. The boundary conditions are symmetry at the back and right faces of the column, inlet at the bottom

3. Computations of the liquid rise in a square capillary: bulk and corner flows Figure 2 shows the computed meniscus shapes in a square capillary of width w = 3 mm for various contact angles (0◦ , 15◦ , 45◦ , 60◦ ) at t = 5 s. For θ ≥ 45◦ the liquid in the corner reaches equilibrium thus approaching a long-time stationary situation. For θ < 45◦ , the computations show an emergence of continuous rivulet flow in the corners. These results agree with the known theoretical predictions [13] on the threshold value of θ < 45◦ for the rivulet flow in the corners of a square capillary. 4

Page 4 of 13

ip t cr us an

(a) w = 1 mm

(b) w = 3 mm

M

Figure 3: Instantaneous snapshots of capillary rise in square capillaries. The snapshots are taken at t = 0.01, 0.1, 1, 10 and 15 s for θ = 30◦ .

20 10

ed

h [mm]

5 2 1

pt

0.5

corner center

0.2 0.1 0.01 0.02

0.05

0.1

0.2

0.5

1

2

5

10

Ac ce

t [s]

(a) w = 1 mm

20 10

h [mm]

5 2 1 0.5 0.2

0.1 0.01 0.02

corner center

0.05

0.1

0.2

0.5

1

2

5

10

t [s]

(b) w = 3 mm Figure 4: Temporal variation of rise height in the center and in the corner for θ = 0◦ . The dashed line indicates the time t90% , when the bulk rise in the center has reached 90% of its equilibrium height.

Figures 3a and 3b show the instantaneous snapshots of capillary rise observed in w = 1 mm and w = 3 mm capillary at θ = 30◦ . As the contact angle satisfies both the general condition for the capillary rise in the bulk (θ < 90◦ ) and in the corner (θ < 45◦ ), we observe both bulk and corner rise (compound capillary flow). Eventually the bulk rise reaches equilibrium, where it attains a constant height and thereafter only the corner rivulet rises. Figures 4a and 4b show the progression of meniscus height measured in the center and in the corner of the capillary over time for w = 1 mm and w = 3 mm respectively. The corner rise always precedes the bulk, unlike the capillary rise observed in wedges where the maximum height slowly approaches the corner from the bulk [32, 45]. Thus the capillary rise in a square capillary has two different regimes defined in [44]: compound capillary flow at the early stages and corner flow at the later stages, when the bulk rise has reached equilibrium. In the early stages of the liquid propagation in a square capillary the bulk rise and corner rise are coupled. This coupling is reflected in the rise, where we observe the corner rise following the same trend as the bulk rise, as shown in the Figures 4a and 4b. During its rising stage, the bulk flow undergoes several regimes [46] in which its speed continuously decreases due to the increasing friction and gravitational pull, eventually reaching equilibrium. As the bulk rise approaches the

5

Page 5 of 13

equilibrium height, at some point of time (t = t s ) it enters the visco-gravitational regime, during which it exponentially relaxes [47] to equilibrium. We consider the duration of this stage to last until t = t s . Here, we take t s as the time to reach 90% of its equilibrium height (t s = t90% ).

10

2

4. Flow in the rivulet in the corner: large times At times larger than t s , the height of the bulk meniscus quickly approaches a plateau value while the rivulet in the corner continues to rise (see Figure 4). This stationary bulk meniscus situation is described by the solution in the corner flow limit [44]. Denote h(T ) as the total height of the rivulet, which is the distance from the central, lowest point of the main bulk meniscus to the rivulet tip. The time T is measured from the instant when the meniscus height reaches 90% of its maximum height, T = t − t s . Figure 5 shows the temporal evolution of liquid height in the corner at θ = 15◦ . A few seconds after the bulk rise reaches equilibrium, we find the corner rise to follow the one-third asymptotic h ∼ T 1/3 . This result is in the agreement with the existing theoretical predictions for an open corner [44] and the experimental observations [30, 32, 45].

1 0.01 0.02

0.05 0.1

0.2

ip t

h [mm]

5

0.5

1

2

5

10

20

2

5

10

20

2

5

10

20

2

5

10

20

cr

T [s]

(a) w = 0.6 mm

us

5

an

h [mm]

10

2

M

1 0.01 0.02

4.1. Scaling of the flow in the rising rivulet The flow in the corner is governed by surface tension, gravity and viscous stresses. The velocity of the rivulet rise in the computations is of order 1 - 5 millimeter per second at the larger times after the bulk meniscus has reached its stationary position. The rivulet thickness is 0.1 - 0.5 millimeter. The corresponding Reynolds number associated with this flow is Re ≈ 0.01 − 0.1. Therefore the inertial effects in the flow in the rivulet can be neglected since they are much smaller than the viscous forces. The shape of the rivulet in the corner is three dimensional. An exact analysis of such flow is therefore rather complicated and so is better left to the simulation results such as shown in Figures 2 and 3. However, far from the bulk meniscus the rivulet is thin and a long-wave approximation can be applied for the estimation of the thickness distribution and the evolution of the rivulet height far from the bulk meniscus. In this case a theoretical description of the rivulet flow using the long-wave, lubrication approximation is straightforward [30, 44, 45]. Below the conventional derivations for the rivulet flow in the framework of the long-wave approximation are presented and examined. Consider a coordinate system with the z-coordinate along the corner with the origin of coordinates at the

0.05 0.1

0.2

0.5

1

T [s]

(a) w = 1 mm 20

h [mm]

5

2

pt

ed

10

1 0.01 0.02

0.05 0.1

0.2

0.5

1

T [s]

Ac ce

(a) w = 3 mm 20

h [mm]

10

5

2

1 0.01 0.02

0.05 0.1

0.2

0.5

1

T [s]

(a) w = 4.5 mm Figure 5: Temporal variation of rise height h in the corner for θ = 30◦ . The dashed line is the one-third asymptotic least squares fit h = c1 T 1/3 + c2 , where c1 and c2 are fitting constants.

6

Page 6 of 13

Remote stationary shape of the rivulet. At large times a portion of the rivulet shape near the bulk meniscus approaches a static shape. This shape can be estimated by equating the expression for pressure (10) with the hydrostatic pressure, leading to the equation for the distribution of the rivulet thickness in the form √ 1 − 2 cos θ + p0 + ρgz = 0, (11) σ δstat (z)

ip t

lowest point of the bulk meniscus. The rivulet thickness is denoted δ(z, T ). See the definition of δ in Figure 6 (left graph). In the long-wave approximation, ∂δ/∂z  1, the curvature of the free interface can be roughly estimated as √ 1 − 2 cos θ ∂2 δ − 2, (8) κ= δ(z, t) ∂z where the first term in the right-hand side of the equation is the curvature of the rivulet horizontal crosssection, obtained from pure geometrical considerations, and the second term is associated with the curvature of the rivulet in the vertical symmetry plane which includes the corner. This second term near the bulk meniscus is determined by the local meniscus radius Ro and can be roughly estimated in the form ∂2 δ/∂z2 ∼ 1/Ro . At large distance from the bulk meniscus the thickness of the rivulet is much smaller than the meniscus radius δ  Ro and the second term can be therefore neglected. We consider only the flow at the distances from the bulk meniscus much larger than the rivulet thickness. The expression for the curvature κ in this remote asymptotic solution can therefore be further simplified for the region where δ  Ro √ 1 − 2 cos θ . (9) κ≈ δ(z, t)

us

cr

where −p0 = 2σ/Ro is the known pressure at the bulk meniscus, Ro is its radius at z = 0. An approximate solution for the static meniscus is Aσ , (12) δstat = ρg(z0 + z) √ 2σ A = 2 cos θ − 1, z0 = , (13) ρgRo

M

an

where z0 is the height of the bulk meniscus. The solution for δstat is positive only if θ < π/4 which is in the agreement with the existing threshold for the rivulet growth [13]. The cross-section area S r of the rivulet in a corner can be estimated from geometrical considerations

ed

Ac ce

pt

In figures 6 and 7 the shapes of interfaces at various heights z are shown for w = 0.6 and w = 4.5 mm, respectively. For both cases the interface shapes are scaled by the corresponding values of δ (shown in right graphs, Figures 6b and 7b). For smaller capillary width the shapes quickly approach the circular shape of the curvature predicted by the approximate expression (9). For the larger capillary shown in Figure 7, w = 4.5 mm, the shapes are also nearly circular, but the deviation from the predicted curvature (9) can be noted. This deviation is explained by the fact that the width of the capillary is comparable with the capillary length. The forces associated with gravity lead to the interface deformation, especially in the regions near the bulk meniscus where the gradient ∂δ/∂z is not negligibly small. Such deviation from the theory is not observed for the capillary width much smaller than the capillary length, shown in Figure 6b. Next, since the velocity component in the direction normal to the solid wall of the capillary is negligibly small, the pressure is uniform in the cross-section normal to the z-axis and can be estimated using the YoungLaplace equation p = σκ(z).

δ2 [4θ − π + 4 cos θ(cos θ − sin θ)] . (14) 4A2 total mass of the rivulet for long times (Mr ≡ RThe ∞ ρ 0 S r dz) is roughly estimated with the help of (14) as Sr =

Mr =

σ2 (2 − π + 4θ + 2 cos θ − 2 sin 2θ) (15) 4g2 ρz0

The force balance equation in the capillary accounts the total weight of the liquid in the capillary and the total capillary force applied to the entire contact line leads to the following equation 4Mr + w2 ρgz0 = 4wσ cos θ.

(16)

The height of the bulk meniscus for the contact angles θ < π/4 is the root of equation (16) obtained in [18] with the help of (15)  √ σ  (17) z0 = 2 cos θ + π − 4θ + 2 sin 2θ . ρgw The average radius Ro of the meniscus is therefore Ro =

2σ 2w = , √ ρgz0 2 cos θ + π − 4θ + 2 sin 2θ

(18)

The value of R0 is uniform at the meniscus surface only in the limit w  a. When the size of the meniscus is comparable or larger than the capillary length a the curvature of the meniscus is not uniform and depends on the local position at the interface.

(10) 7

Page 7 of 13

3.0

0.15

z

= 0.1 mm

z

= 0.2 mm

z

= 0.3 mm

2.5

mm

= 0.1 mm

z

= 0.2 mm

z

= 0.3 mm

theory

2.0

1.5

cr

0.10

z

1.0 0.05

)

,t (z

0.00

0.0 0.05

0.10

0.15

0.20

0.0

0.5

1.0

an

mm

us

0.5

0.00

ip t

0.20

(a)

1.5

2.0

2.5

3.0

(b)

ed

M

Figure 6: Shapes of the rivulet horizontal cross-section at various distances from the bulk meniscus for w = 0.6 mm, θ = 15◦ . (a) Dimensional shapes and (b) the shapes are scaled by the corresponding value of δ and compared with the theoretically predicted circular interface of constant curvature, determined in equation (9).

2.0

z

= 0.2 mm

z

3.0

z

= 0.1 mm

z

= 0.2 mm

= 0.3 mm

z

= 0.3 mm

z

= 0.4 mm

z

= 0.4 mm

z

= 0.5 mm

z

= 0.5 mm

Ac ce

mm

1.2

= 0.1 mm

pt

1.6

z

2.5

2.0

theory

1.5

0.8

1.0

0.4

0.5

0.0

0.0

0.4

0.8

1.2

1.6

0.0 2.0

0.0

0.5

1.0

1.5

2.0

2.5

3.0

mm

(a)

(b)

Figure 7: Shapes of the rivulet horizontal cross-section at various distances from the bulk meniscus for w = 4.5 mm, θ = 15◦ . (a) Dimensional shapes and (b) the shapes are scaled by the corresponding value of δ and compared with the theoretically predicted circular interface of constant curvature, determined in equation (9).

8

Page 8 of 13

subject to the boundary conditions

Dynamics of the flow in the rivulet. Consider now the flow in a thin rivulet. The momentum balance equation of the flow in the rivulet can be obtained neglecting the inertial effects ! ∂p ∂2 u ∂2 u + = 0, (19) − − ρg + µ ∂z ∂x2 ∂y2

F(0) = 0,

ip t

us

(20)

M

an

where B is a dimensionless coefficient of order unity, which depends only on the contact angle θ. The value of B is associated with the friction factor [44], which is often introduced to account for the shape of the crosssection of the rivulet. U(z, T ) is the characteristic longitudinal velocity of the liquid in a rivulet cross-section. This velocity can now be estimated from (20) with the help of (9) and (10) gρh2 σAδz U=− − Bµ Bµ

4.2. Rivulet profile δ(z)

(21)

ed

The mass balance equation of the rivulet yields ∂δ2 ∂δ2 U + = 0. ∂T ∂z

(22)

Ac ce

pt

The equation for the evolution of the rivulet thickness is obtained by substitution of (21) in (22)   4gρδz δ2 − 2BµδT + Aσ δzz δ + 2δ2z = 0 (23) This equation has a similarity solution far from the bulk meniscus of the capillary: h i δ(z, T ) = LF W − G(z + z0 )T −1/3 T −1/3 . (24) L

=

G

=

A1/3 B1/3 µ1/3 σ1/3 , 2 × 31/3 g2/3 ρ2/3 2B1/3 g1/3 µ1/3 ρ1/3 , 31/3 A2/3 σ2/3

(25) (26)

where W is an unknown constant and F(ξ) is a dimensionless function of a similarity variable ξ = W − G(z + z0 )T −1/3 .

(29)

cr

U(z, T ) ∂p − ρg − Bµ = 0, ∂z δ(z, T )2

4 , W −ξ

since ξ = 0 corresponds to the rivulet tip and the region ξ → W corresponds to the static solution (12). It is clear from the form of similarity solution that this instationary flow cannot be matched directly with the bulk meniscus. Therefore the similarity can appear only if the rivulet flow is subdivided into three main regions: (i) the similar flow (24) near the rivulet tip, (ii) the intermediate region which at large times approaches the static solution (12) and (iii) the complex three-dimensional shape of the rivulet near the bulk meniscus. The regions (i) and (ii) are matched by the second condition in (29). It is therefore difficult to obtain exact values for W from this analysis since the equations are derived on the long-wave approximation which is not valid near the meniscus. This is why full three-dimensional computations are necessary for a reliable solution of this problem.

where u is the streamwise velocity component. Here only the major viscous terms are taken into account, as in the lubrication approximation, valid for a long rivulet. Equation (19) can be simplified assuming the similarity of the flow in a rivulet cross-section. In this case the momentum equation is reduced to −

F(ξ → W) =

The dimensionless rivulet profile along the corner, F(ξ), is estimated from (24) using the values for δ(T ) from the full numerical computations. The unknown dimensionless parameter B(θ) is fitted using the condition that F(ξ) does not depend on θ. The results of computations of F(ξ) are shown in Figure 8. The predicted shapes F(ξ) coincide on the segment ξ < 1. At larger values of ξ the different profiles diverge, indicating the region where the similarity solution is not valid. This is not surprising since the similarity solution is not able to satisfy the matching conditions at the bulk meniscus. The dimensionless rivulet shape for various liquid viscosities and capillary widths are shown in Figure 9. In all the cases the similarity solution is valid near the rivulet tip (for small values of ξ) and breaks up near the static bulk meniscus. The effect of the bulk meniscus is especially clear in the case of different capillary widths. This effect leads to the dependence of the rivulet growth on the value of w. 4.3. The rate of the rivulet propagation The long time asymptotics for the position of the tip of the rivulet at ξ = 0 can now be determined from (24) as

(27)

Function F(ξ) can be obtained by integration of the ordinary differential equation

h=K

−Fξ F 2 + Fξξ F + F + 2Fξ2 − WFξ + ξFξ = 0. (28)

A2/3 σ2/3 B1/3 g1/3 µ1/3

T 1/3 − z0 ,

(30)

9

Page 9 of 13

5

12.5

4

10

3

θ = 15◦

H/a

F (ξ)

7.5

2

5 2.5

µ = 19.2 mPas µ = 97 mPas

0 0

0.5

1

1.5

2

2.5

0

ξ

5

10

(a)

ip t

θ = 0◦ θ = 15◦ θ = 30◦ theory

1

θ = 30◦

15

20

25

(σT /µa)1/3

cr

10

Figure 10: Effect of viscosity and contact angle on the dimensionless corner rise h/a for the width of the capillary w = 3 mm.

8

us

B

6

3 4

2.5 2

an

2 0 10

20

30

40

K

0

θ [◦ ]

(b)

1.5 1

Figure 8: Dimensionless rivulet profile F(ξ) defined in equation (24) for w = 3 mm, µ = 19.2 mPa·s in (a), obtained using the fitted values of the parameter B(θ), shown in (b). The numerical predictions are compared with the theory (26) computed for W = 3.4

θ = 0◦ θ = 15◦ θ = 30◦

M

0.5 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

δo /a

ed

5

4

F (ξ)

3

pt

2

1

where K is a dimensionless rate of the rivulet growth. This relation is similar to the corresponding expression obtained in [44].

0

Ac ce

µ = 19.2 mPas µ = 97 mPas

0

0.5

1

1.5

2

2.5

Equation (30) can be rewritten in the form h/a ∼ (σT/µa)1/3 . This relation, proposed in [32] describes well the evolution of the rivulet height in a corner formed by two vertical cylinders. Figure 10 presents the temporal variation of the dimensionless rivulet height in a corner of the 3 mm capillary at two values of the contact angle: θ = 15◦ and 30◦ . The curves for the same contact angle collapse into one curve for both viscosities.

ξ

(a) θ = 15 , w = 3 mm. ◦

5

4

F (ξ)

3

2

w w w w

1

= 0.6 mm = 1 mm = 3 mm = 4.5 mm

The predicted dimensionless rate of the rivulet growth K, defined in (30), is of order unity and slightly depends on the width w of the capillary through the influence of the shape of the bulk meniscus. We assume that the influence of w is introduced through the dependence of the initial rivulet thickness δo near the bulk meniscus, which can be estimated from geometri-

0 0

0.5

1

1.5

2

Figure 11: Computed values for the dimensionless rate K of the rivulet rise in a corner, defined in (30). Dependence on the wall contact angle and the dimensionless initial rivulet thickness δo /a.

2.5

ξ

(b) θ = 15◦ , µ = 19.2 mPa·s. Figure 9: Dimensionless rivulet profile F(ξ) defined in equation (24) for various viscosities (9a) and capillary widths (9b).

10

Page 10 of 13

cal considerations √ 2 δo ≈ w − Ro . (31) 2 The influence of the value of δo is not reflected in the scaling (30) since the effect of the second derivative ∂2 δ/∂z2 is neglected in the long-wave analysis, which make it impossible to match the solution to the bulk meniscus region. The dependence of the computed values of K on the ratio δo /a is shown in Figure 11. The dependence of K on δo /a is indeed relatively weak, but the rate of the rivulet growth slightly increases with δo /a. It is remarkable that the obtained values of K almost do not depend on the contact angle.

by the capillary length a is obtained from the full numerical simulations of the flow. It is shown that the function K(δo /a) is almost independent of the contact angle, liquid viscosity or surface tension.

ip t

Acknowledgement

cr

This work has been performed in the framework of the Marie Curie Initial Training Network ”Complex Wetting Phenomena” (CoWet), Grant Agreement no. 607861. The authors gratefully acknowledge the computing time granted on the Lichtenberg-High Performance Computer at Technische Universit¨at Darmstadt.

us

[1] S. Girardo, R. Cingolani, S. Chibbaro, F. Diotallevi, S. Succi, D. Pisignano, Corner liquid imbibition during capillary penetration in lithographically made microchannels, Appl. Phys. Lett. 94 (17) (2009) 2007–2010. doi:10.1063/1.3123804. [2] T. W. Sowers, R. Sarkar, S. E. Prameela, E. Izadi, J. Rajagopalan, Capillary driven flow of polydimethylsiloxane in open rectangular microchannels, Soft Matter 12 (12) (2016) 5818–5823. [3] J. Feng, J. P. Rothstein, Simulations of novel nanostructures formed by capillary effects in lithography, J. Colloid Interface Sci. 354 (1) (2011) 386–395. doi:10.1016/j.jcis.2010.10.030. [4] N. R. Morrow, G. Mason, Recovery of oil by spontaneous imbibition, Current Opinion in Colloid & Interface Science 6 (4) (2001) 321–337. [5] M. Piri, M. J. Blunt, Three-dimensional mixed-wet random pore-scale network modeling of two-and three-phase flow in porous media. i. model description, Physical Review E 71 (2) (2005) 026301. [6] T. K. Tokunaga, J. Wan, Water film flow along fracture surfaces of porous rock, Water Resources Research 33 (6) (1997) 1287– 1295. [7] A. Yiotis, A. Boudouvis, A. Stubos, I. Tsimpanogiannis, Y. Yortsos, Effect of liquid films on the drying of porous media, AIChE Journal 50 (11) (2004) 2721–2737. [8] J. Schoelkopf, C. J. Ridgway, P. A. Gane, G. P. Matthews, D. C. Spielmann, Measurement and network modeling of liquid permeation into compacted mineral blocks, Journal of colloid and interface science 227 (1) (2000) 119–131. [9] B. Taylor, Concerning the ascent of water between two glass plates, Philos. Trans. 27 (27) (1710) 538. [10] F. Hauksbee, An Account of an Experiment Touching the Ascent of Water between Two Glass Planes, in an Hyperbolick, Philos. Trans. (27) (1710) 539–540. doi:10.1098/rstl.1763.0053. [11] Report of the Fourth Meeting of the British Association for the Advancement of Science, J. Murray., London, 1835. [12] I. Newton, Opticks: or, A Treatise of the Reflexions, Refractions, Inflexions and Colours of Light, William Innys, London, 1730. [13] P. Concus, R. Finn, On the behavior of a capillary surface in a wedge, Appl. Math. Sci. 63 (2) (1969) 292–299. doi:10.1007/BF02786711. [14] G. Mason, Capillary Behaviour of a Perfectly Wetting Liquid in Irregular Triangular Tubes, J. Colloid Sci. 141 (1). [15] H. Wong, S. Morris, C. J. Radke, Three-dimensional menisci in polygonal capillaries, J. Colloid Interface Sci. 148 (2) (1992) 317–336. doi:10.1016/0021-9797(92)90171-H. [16] N. Ichikawa, K. Hosokawa, R. Maeda, Interface motion of capillary-driven flow in rectangular microchan-

5. Conclusions

Ac ce

pt

ed

M

an

In this study, the capillary rise of a wetting liquid in a square capillary is investigated using numerical simulations. The early stage of capillary rise is dominated by the bulk flow with the corner rise preceding the bulk rise for θ < π/4. In this stage, the amount of liquid rising in the corner is smaller than at later times; hence, its influence on the bulk flow can be neglected as the bulk meniscus evolves. When the bulk rise starts relaxing towards a static height during the visco-gravitational regime, the corner rise starts its transition towards pure corner flow. At longer times, the capillary rise in the corner is found to rise according to the one-third asymptotic: h ∼ T 1/3 which is observed in the previous experimental studies and predicted in several theoretical studies [30, 32, 44, 45]. The flow in the rivulet is described using the longwave approximation, like in many previous theoretical studies of this or similar problems [30, 44, 45]. We have shown that the rivulet flow can be subdivided into three regions. The flow near the rivulet tip is described using a similarity solution. This similarity solution is matched with the intermediate region which at large times approaches the static solution. The static solution is connected with the bulk meniscus region. However the influence of the meniscus geometry is not considered in the analysis since the long-wave approximation is not valid in this region. The full numerical computations of the rivulet profile δ(z, T ) agree very well with the theoretical predictions in the regions far from the bulk meniscus. Finally a scaling for prediction the rate of capillary rise is proposed. The dependence of the dimensionless growth rate K on the initial rivulet thickness δo , scaled

11

Page 11 of 13

[23] [24]

[25]

[26]

[27] [28]

[29]

[30] [31]

[32]

[33]

[34]

[35]

[42]

[43]

ip t

cr

[41]

us

[22]

[40]

an

[21]

[39]

[44] [45]

M

[20]

[38]

[46] [47]

ed

[19]

[37]

pt

[18]

[36]

Methods in Numerical Dynamics, Vol. 1000, IUC Dubrovnik, Croatia, Dubrovnik, Croatia, 2007, pp. 1–20. Open∇Foam. the open source cfd toolbox, https://www.openfoam.com/. H. Weller, A new approach to vof-based interface capturing methods for incompressible and compressible flow, Tech. Rep. 04, OpenCFD Ltd (2008). J. Brackbill, D. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335 – 354. C. Kunkelmann, Numerical Modeling and Investigation of Boiling Phenomena, Ph.D. thesis, Technische Universit¨at Darmstadt (2011). S. Batzdorf, Heat transfer and evaporation during single drop impingement onto a superheated wall, Ph.D. thesis, Technische Universit¨at Darmstadt (2015). S. Afkhami, S. Zaleski, M. Bussmann, A mesh-dependent model for applying dynamic contact angles to VOF simulations, J. Comput. Phys. 228 (15) (2009) 5370–5389. J. B. Dupont, D. Legendre, Numerical simulation of static and sliding drop with contact angle hysteresis, J. Comput. Phys. 229 (7) (2010) 2453–2478. doi:10.1016/j.jcp.2009.07.034. D. Legendre, M. Maglio, Comparison between numerical models for the simulation of moving contact lines, Comput. Fluids 113 (2014) 2–13. doi:10.1016/j.compfluid.2014.09.018. M. M. Weislogel, Compound capillary rise, J. Fluid Mech. 709 (2012) 622–647. doi:10.1017/jfm.2012.357. F. J. Higuera, A. Medina, A. Li˜na´ n, Capillary rise of a liquid between two vertical plates making a small angle, Phys. Fluids 20 (10) (2008) 102102. doi:10.1063/1.3000425. N. Fries, M. Dreyer, The transition from inertial to viscous flow in capillary rise., J. Colloid Interface Sci. 327 (1) (2008) 125–8. N. Fries, M. Dreyer, An analytic solution of capillary rise restrained by gravity., J. Colloid Interface Sci. 320 (1) (2008) 259– 63.

Ac ce

[17]

nel, J. Colloid Interface Sci. 280 (1) (2004) 155–164. doi:10.1016/j.jcis.2004.07.017. F. F. Ouali, G. McHale, H. Javed, C. Trabi, N. J. Shirtcliffe, M. I. Newton, Wetting considerations in capillary rise and imbibition in closed square tubes and open rectangular crosssection channels, Microfluid. Nanofluidics 15 (3) (2013) 309– 326. doi:10.1007/s10404-013-1145-5. J. Bico, D. Qu´er´e, Rise of liquids and bubbles in angular capillary tubes., J. Colloid Interface Sci. 247 (1) (2002) 162–166. doi:10.1006/jcis.2001.8106. K. S. Birdi, D. T. Vu, A. Winter, A. Nørregård, Capillary rise of liquids in rectangular tubings, Colloid Polym. Sci. 266 (5) (1988) 470–474. doi:10.1007/BF01457265. P. Wu, H. Zhang, A. Nikolov, D. Wasan, Rise of the main meniscus in rectangular capillaries: Experiments and modeling, J. Colloid Interface Sci. 461 (2016) 195–202. doi:10.1016/j.jcis.2015.08.071. M. Dong, I. Chatzis, The Imbibtion and Flow of a Wetting Liquid along the Corners of a Square Capillary Tube, J. Colloid Interface Sci. 172 (1995) 278–288. A. Singhal, W. Somerton, Two-phase flow through a noncircular capillary at low Reynolds numbers, J. Can. Pet. Technol. P. S. Ayyaswamy, I. Catton, D. K. Edwards, Capillary Flow in Triangular Grooves, J. Appl. Mech. (1974) 332–336. R. Lenormand, C. Zarcone, Role of roughness and edges during imbibition in square capillaries, in: SPE annual technical conference and exhibition, no. SPE-13264-MS, Society of Petroleum Engineers, Houston, Texas, 1984. doi:10.2118/13264-MS. H. B. Ma, G. P. Peterson, X. J. Lu, The influence of vaporliquid interactions on the liquid pressure drop in triangular microgrooves, Int. J. Heat Mass Transf. 37 (15) (1994) 2211–2219. doi:10.1016/0017-9310(94)90364-6. T. Ranshoff, C. Radke, Laminar Flow of a Wetting Liquid along the Corners of a Predominantly Gas-Occupied Noncircular Pore, J. Colloid Interface Sci. 121. M. M. Weislogel, S. Lichter, Capillary flow in an interior corner, J. Fluid Mech. 373 (1998) 349–378. Y. Chen, M. M. Weislogel, C. L. Nardin, Capillary-driven flows along rounded interior corners, J. Fluid Mech. 566 (2006) 235. doi:10.1017/S0022112006001996. M. M. Weislogel, J. A. Baker, R. M. Jenson, Quasi-steady capillarity-driven flows in slender containers with interior edges, J. Fluid Mech. 685 (2011) 271–305. doi:10.1017/jfm.2011.314. L.-H. Tang, Y. Tang, Capillary rise in tubes with sharp grooves To cite this version :, J. Phys. II 4 (1994) 881–890. G. J. Gutierrez, A. Medina, F. J. Higuera, Equilibrium height of a liquid in a gap between corrugated walls under spontaneous capillary penetration, J. Colloid Interface Sci. 338 (2) (2009) 519–22. doi:10.1016/j.jcis.2009.06.056. A. Ponomarenko, D. Qu´er´e, C. Clanet, A universal law for capillary rise in corners, J. Fluid Mech. 666 (2011) 146–154. doi:10.1017/S0022112010005276. S. Girardo, S. Palpacelli, A. De Maio, R. Cingolani, S. Succi, D. Pisignano, Interplay between Shape and Roughness in EarlyStage Microcapillary Imbibition, Langmuir 28 (5) (2012) 2596– 2603. S. Son, L. Chen, Q. Kang, D. Derome, J. Carmeliet, Contact Angle Effects on Pore and Corner Arc Menisci in Polygonal Capillary Tubes Studied with the Pseudopotential Multiphase Lattice Boltzmann Model, Computation 4 (1) (2016) 12. doi:10.3390/computation4010012. H. Jasak, A. Jemcov, Z. Tukovic, Openfoam: A C++ library for complex physics simulations, in: Int. Workshop on Coupled

12

Page 12 of 13

ce

g

Ac

Rivulet flow

pt

ed

Spontaneous rise in a square capillary

bulk meniscus

Page 13 of 13

θ = 0◦

θ = 15◦

θ = 45◦

θ = 60◦