© N. Abbaspour et al., EDP Sciences, 2020

1 Introduction

Small scale structures added to otherwise smooth surface of the same chemical composition influence the wettability [1,2]. Adapting structure geometry may enhance water drop evaporation [3] or coalescence [4,5]. Here we concentrate our attention on structures that cause spontaneous water wicking in a privileged direction as observed on the skin of organisms living in arid areas [6–8]. The phenomenon is also the subject of the experiments of [9] that demonstrate spontaneous directed liquid wicking on artificial flat surface structured by added series of micro-metric scale parallel channels. Each of the latter can be viewed as a succession of identical unit cells composed of straight part alternating with non straight wider part. Moreover, flat surface and channel walls are made of the same material. The reference [9] points out the leading role of the capillary forces, summarized by the contact angle of the solid/liquid/air triplet.

This paper studies liquid wicking on surfaces shaped with the geometry considered in [9]. The chosen approach is three dimensional two phase flow numerical simulation that allows us studying the dynamics of liquid wicking, not addressed in previous studies. This allows us investigating time/space scales as small as desired, as well as the influence of contact angle and pattern scale. An experiment validates this approach. It is performed on the geometry studied in [9], but the material that composes the solid structure is now different: it is stainless steel in view of industrial application. In this aim, we moreover consider a scaling compatible with selective laser melting (SLM), a very promising technique that allowed us to manufacture the solid structure on which we performed our liquid wicking experiment. It turns out that the technique imposes a scaling significantly larger than in [9]. Therefore, the question of whether spontaneous liquid wicking occurs in these conditions deserves being answered since thinner patterns are known to speed up wicking and to enhance the privileged character of one specific direction [6].

Section 2 specifies the here considered channel structure. Section 3 sets the mathematical model that allows us capturing the evolution of the liquid/air interface on the structure. The simulated dynamic is discussed in Section 4. Section 5 summarizes the experiment that visualizes the hatched dynamic observed in the simulation.

2 Channel geometry

Starting from horizontal flat surface, we modify it by building vertical walls limiting micro-channels that mimic the shape of fins (or arrow heads) as the ones on which [9] observed one-directional wicking. Flat surface and channel walls are made of the same solid material. Figure 1 represents a view from above of an example of a flat surface to which such micro-channels were added. Figure 2 represents a CAD perspective view of the same. On these figures, we see that all channels are identical and have parallel axes of x direction. However, their width is not uniform: they enclose voids that are periodic (in x direction) copies of the same elementary cell, the basis of which is represented in white on Figure 3. The elementary cell is reminiscent of a pair of fins (or arrow head) followed by straight channel. Numerical simulation and experimental validation are performed on three variants A, B, and C of this geometry, according to Figure 3 and Table 1. However, we keep walls height constant, and equal to 0.1 mm.

We study liquid wicking on the modified surface by numerically simulating the evolution of the liquid/air interface, and also of the liquid velocity field. An experiment validates the simulation though boundary conditions are not exactly the same. Indeed, memory and computing time cause that the simulation can only be run in a smaller and simplified domain (a single channel complemented by a volume above it limited by a definite horizontal plane). Periodic lateral boundary conditions attempt to account for the existence of other channels. The experiment, by contrast, cannot be performed on single channel: instead it uses three stainless steel structures (A,B,C) exhibiting between 5 and 15 parallel identical channels with the parameters specified in Table 1, the same as in simulation. The view from above displayed in Figure 1 is a microscopic picture of structure A.

Fig. 1 View from above of horizontal stainless steel surface modified by super-imposing fin (or arrow head) shaped parallel channels. Channel bottom/walls appear in withe/dark. We will see that water drops preferentially wick in x direction, the opposite to the one to which the arrows point.

Fig. 2 Computer-aided design perspective view of the structure pictured on Figure 1. Observe the two directions of the faces of the channel walls.

Fig. 3 The basis of the elementary cell of a channel in white. Channel walls are in blue. In our study lengths L1 = L2 and Wi take the values specified in Table 1, and φ = 27°.

3 Mathematical model and numerical simulation

In a single channel we simulate three dimensional two phase flow according to a Volume of Fluid (VOF) method [10] including triple line modeling.

3.1 Mathematical model underlying VOF method

Volume of Fluid approach [11,12] assumes one (two phase) fluid of velocity $\overrightarrow{u}$ and volume density ρ satisfying $$\rho =\alpha {\rho }_{liq}+\left(1-\alpha \right){\rho }_{gas}$$(1)

in which α and ρ_liq are liquid phase volume fraction and volume density while the volume density of the gaseous phase is ρ_gas. Moreover, the VOF approach assumes two phase fluid viscosity μ related to liquid and gaseous phase viscosities μ_liq and μ_gas according to $$\mu =\alpha {\mu }_{liq}+\left(1-\alpha \right){\mu }_{gas}.$$(2)

In this model, fluid velocity $\overrightarrow{u}$ and liquid volume fraction α evolve according to $$\begin{array}{l}\frac{\partial (\rho \overrightarrow{u})}{\partial t}+\nabla .\left(\rho \overrightarrow{u}\otimes \overrightarrow{u}-\mu \left(\nabla \overrightarrow{u}+\nabla {\overrightarrow{u}}^T\right)\right)\hfill \\ =-\nabla p+{\overrightarrow{f}}_{\sigma }+\rho \overrightarrow{g}\hfill \end{array}$$(3)

and $$\frac{\partial \alpha }{\partial t}+\nabla .\left(\alpha \overrightarrow{u}\right)=0.$$(4)

The momentum equation (3) includes the continuous surface force ${\overrightarrow{f}}_{\sigma }$ introduced by [11]. The variable α is assumed to continuously vary between 0 and 1 in a neighborhood of the liquid gas interface. There, the continuous surface force is given by $${\overrightarrow{f}}_{\sigma }=\sigma \frac{\rho \kappa \nabla \alpha }{\frac{1}{2}\left({\rho }_{liq}+{\rho }_{gas}\right)},$$(5)

where κ is a continuous function representing the free surface curvature. It is defined in the interface vicinity by $$\kappa =\nabla .\overrightarrow{n}.$$(6)

In this equation $\overrightarrow{n}$ is continuous also and represents the unit vector normal to the interface. It satisfies $$\overrightarrow{n}=\frac{\nabla \alpha }{\Vert \nabla \alpha \Vert }.$$(7)

Moreover, the wettability prescribes $\overrightarrow{n}$ at the tripleline that is the intersection of the solid surface with the liquid/gas interface.

3.2 Triple line modeling

At solid boundaries we consider wall adhesion. At the triple air/liquid/solid interface there is the so-called triple line. At equilibrium, the static contact angle θ characterizes the wettability. Out of equilibrium, the macroscopic contact angle in general depends on the front dynamics. In principle, it differs from the static contact angle [13]. However, at the mesoscopic scale, it is still relevant to consider the latter contact angle θ [14]. Therefore, following [11] we assume that the normal $\overrightarrow{n}$ to the interface satisfies $$\overrightarrow{n}=\cos \theta {\overrightarrow{n}}_w+\sin \theta {\overrightarrow{n}}_t$$(8)

in which θ is the equilibrium contact angle. Vector ${\overrightarrow{n}}_w$ is the unit vector normal to the wall, and ${\overrightarrow{n}}_t$ is tangent to the wall and normal to the contact line. At the triple line, equations (7) and (8) imply $$\cos \theta {\overrightarrow{n}}_w+\sin \theta {\overrightarrow{n}}_t=\frac{\nabla \alpha }{\Vert \nabla \alpha \Vert }.$$(9)

The code deduces κ and ${\overrightarrow{f}}_{\sigma }$ from this equation, in view of (5)–(7).

3.3 Liquid inlet, and other specifications

We consider two possibilities for liquid inlet. The first one is equivalent to initial condition: at time instant t = 0 we put on the channel a small cuboid. It is 0.3 mm high and its width is W2. Cuboid length is between one and two elementary cells. The volume in the channel below the cuboid is filled up as well. In order to increase the liquid volume, a second inlet condition imposes a localized liquid flux (at channel bottom).

3.4 Dimensional analysis

For a given geometry, the fluid evolution is governed by liquid density, viscosity and surface tension (ρ, μ, σ) and by the static contact angle θ that quantifies fluid/solid interactions. This dynamics is ruled by Reynolds number $Re=\frac{\rho VW1}{\mu }$ (V characteristic velocity, W1 defined by Fig. 3), Capillary number $Ca=\frac{\mu V}{\sigma }$, Weber number $We= \frac{\rho V^{2}W1}{\sigma }$ and Bond number $Bo=\frac{W_1}{{\mathcal{l}}_c}$ with ℓ_c = 2.7 mm the capillary length of water. All simulations reveal a maximum velocity V of 1 m/s, hence a Reynolds number about 100 showing that inertia plays a significant role, and that the flow is far from equilibrium. Inertia force competes with surface tension since the Weber number is about 1. Moreover, Capillary and Bond numbers much smaller than 1 suggest surface tension dominating over viscous and gravity forces.

4 The dynamics of drop wicking

Simulating the evolution of the liquid/air interface reveals strong differences in liquid progress in the two directions x and − x. We outline the most significant trend observed on the liquid/air interface dynamic started from an initial amount of liquid set on single channel of A type, the solid/liquid/air contact angle θ being of 50°. Then, assuming the same initial condition, we discuss the influence of θ and of the elementary channel cell size. Second liquid inlet condition is considered at the end of the section.

4.1 Liquid/air interface motion in x and −x directions, at 50° contact angle

An amount of liquid initially set on a single channel (e.g. A type) spreads differently in x and − x directions.

The simulation reveals more rapid liquid/air interface motion in x direction. We see it on Figure 4 composed of four sub-figures. Each one displays snapshots of liquid volume and velocity field (at bottom an top, respectively). The figure also reveals quite uneven liquid/air interface velocity: interfaces proceed faster in straight channels, more rapidly in x than in − x direction. Flooding the fins takes more time, even more for the interface that moves in − x direction (i.e. to the right). Moreover, with the initial liquid amount representing the initial condition of the simulation that returned Figure 4, the right front stops after having flooded one elementary cell plus one straight channel. In the VOF simulation it definitively stops at next fin entrance. The left front, oppositely, continues to the left. Except regarding right front pinning, these findings agree with the visualizations described in Section 5.

Deeper inspection of the simulation shows that even the left front takes pauses near each next straight channel entrance while an eddy dissipates kinetic energy in fin tails. The snapshots of the velocity field at t = 1.56 ms and t = 1.94 ms in Figure 4 reveal the eddy. The successive pauses taken by the left front have increasingly longer durations (not represented on the figure). However, the left front floods the left part of the domain represented on Figure 4. Upon increasing simulation time and channel length we observe a left front that stops before reaching channel end. We conjecture that smallertotal liquid volume causes fronts to stop earlier. We discuss this point in Section 4.4.

Fig. 4 Simulated liquid velocity and free surface height during the wicking of a parallelepiped of liquid over a single fin shaped channel, using the VOF. The top of each sub-figure documents the velocity, and just below we see liquid height distribution at z = 0.05 mm above channel bottom. Contact angle is 50°. The simulation domain consists of 6 elementary cells. The initial liquid parallelepiped is one elementary cell long. Dashed lines mark its limits. Just after t = 0 the liquid parallelepiped deforms quickly to smoother shape and wicks on channel bottom in the two directions x and − x. At t = 0.34 ms it has invaded the straight channel at the right, and is still progressing in the wider part at the left. However, at t = 0.64 ms the left front has compensated its delay. At t = 1.56 ms the left front has clearly covered more distance than the right front. At t = 1.94 ms the right front spreads a little further than $E_r^3$ before immediately receding because of surface tension that tends to decrease the free surface. Hereafter the right front remains pinned. On the other side, differently, we see that the flow suddenly accelerates in the fins since the time the left front touches $E_l^3$ at t = 1.56 ms. Between time instants t = 1.56 ms and t = 1.94 ms, the acceleration is caused by the tapered geometry of the fin that increases the wet solid surface over liquid volume ratio. At t = 1.94 ms, the liquid completely wets the fin surface but the kinetic energy is only partly dissipated. This results into reverse flow in the fins, converging to the middle of the channel. The volume increase in the center of the channel increases the free surface area: this transforms a fraction of the kinetic energy into surface energy. Relaxation then causes water to spread again, and to flood the next straight channel which acts as capillary.

4.2 Influence of contact angle

Since capillary forces are smaller at larger contact angle slower wicking is expected in this case. This is corroborated by the simulation.

Starting the VOF simulation from the cuboid mentioned in Section 4.1, we observe that increasing the static contact angle θ decreases the mean front velocity in x direction. No water spreads through the channel beyond θ_lim ≥ 70°. The liquid volume converges rapidly to an equilibrium state reminiscent of the steady ridge observed on a flat substrate [15] with periodic lateral boundary conditions.

Conversely, decreasing θ from 50° speeds up the wicking front progression in x direction. However, this affects the front pinning in the −x direction. Indeed, at θ = 30° the fluid spreads over a larger area before retracting and definitively stopping. This dynamic results froma competition between inertia and surface tension. There exists a critical value, noted θ_c ≤ 30°, below which the two fronts just take pauses before proceeding further so that water spreads in both directions. Even in this case the wicking is slower in −x than in x direction. Though the front is no longer pinned, the edges constrain the flow to spread out by forming a thin film. Since the viscous dissipation is high for the spreading of thin films, the kinetic energy almost totally dissipated near the right front explains the slower wicking in −x direction.

4.3 Influence of geometry

Simulating wicking on structures B and C reveals that the geometry does not qualitatively change the above described scenario, within the here considered limits. Indeed, the three structures exhibit similar trend for the contact angle dependency of the average velocity, showing in all cases a superior limit θ_lim for the wicking in the two directions and a critical value θ_c below which the pinning effect vanishes in −x direction.

However, pattern size influences wicking velocity, found larger on smaller patterns, as observed in [6]. Since the transport is driven by the capillary force related to the surface energy of the solid, down-scaling decreases the capillary force to inertia ratio, larger in smaller cell.

4.4 Continuous flow

To increase the total volume of liquid in simulation, we replace the initial condition by a continuous volume rate of 1.06 mm³∕s at channel bottom.

With a contact angle of 50°, the simulation returns a dynamic that follows the trends described in Section 4.1, except that now the spread in the positive direction is no longer limited. This confirms the conjecture claimed at the end of Section 4.1, regarding the influence of the total amount of liquid. However, the VOF simulation still exhibits definitively stopped right front, except for θ < 30°. This may be an issue of the front diffusion caused by the interface modeling included in VOF method that excessively increases the transition region at the interface and undervalues the capillary force. This may cause more damage in the direction that involves more severe obstacles and drains less liquid. Indeed, the main drawback of the VOF simulation is exaggerating the numerical interface diffusion, though we use Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM), a better choice than Geometric Reconstruction Scheme (GRS) and High Resolution Interface Capturing (HRIC).

The color gradient Lattice Boltzmann Method (LBM) described in [16] restricts interface diffusion, and does not exhibit the right front pinning observed in the VOF simulation: the right front still stops at the same sharp corner at fins entrance. However, a small amount of water is observed in the tip of the fin. The liquid volume in this region growths andprogressively fills the fins. Then, the water continues to the next channel. Except regarding right front progression, the LBM reveals the same trends observed with the VOF method and described in the two previous sections.

5 Experimental validation

We checked the qualitative features revealed by the simulation upon visualizing drop wicking on stainless steel structures exhibiting series of parallel identical channels as described in Section 2, and built on the elementary cells A, B and C specified in Figure 3 and Table 1. We see the three structures in Figure 5.

The stainless structures were manufactured by SLM, a technique that allowed us to build channel walls out of stainless steel powder on a flat plate of the same material. Photos and videos captured single distilled water drop evolution on each of the three structures, with the help of Keyence VHX-600 microscope equipped with VH-Z20R and VH-Z100R lens. The steel plate with the structures was horizontal and the microscope was positioned above. Photos and videos were captured using Samsung Galaxy Note 8. In all tests we observed contact angles of about 50 degrees for the distilled water/air/stainless steel triplet.

Figure 5 displays four snapshots representing the different steps (t₀, t₁, t₂, t₃) of the evolution of single drop set on the middle of structure B, the one with the widest channels. Drop wicking on other structure only differs by the time instants at which the liquid completely wets the structure in x or − x direction. On each structure wicking was found to start as soon as the drop was set on it at initial time t₀. We quantify the mean velocity in x or − x direction by recording  x₀, x_f and t_f that represent initial drop meniscus position and structure end in the considered direction, and the time at which the fluid reaches structure end. With these notations, the average velocity in each direction is $$\overline{V}=\frac{x_f-x_0}{t_f-t_0}.$$(10)

Not surprisingly, we observe larger such velocities in both directions x and − x on structure A whose smaller pattern enhance capillary forces that drive front progression. However, on each structure the absolute value of the average velocity in − x direction is nearly half the value observed in x direction. This is closely related to the details of the wicking dynamics. Indeed, Figure 5 shows fast front progression, but slowing down the videos suggests quite uneven dynamics, with much slower steps between rapid ones, especially in − x direction. This is reminiscent of the pauses taken by the simulated liquid fronts according to Section 4 at cross sections exhibiting a wedge with an obtuse angle ($E_r^{2i+1}$ and $E_l^{2i+1}$ on Fig. 4). There, the front adapts its shape to satisfy the contact angle requirement beyond the wedge: the convex front that proceeds to the right experiences a more drastic evolution that takes more time, during which it stays pinned by the wedge. LBM simulations at θ = 50° reveal thin films that slowly bypass the wedge and help the right front to escape. These films resemble the ones observed on VOF simulations at θ < 30° at the end of Section 4.2. The latter discussion and the comparison with the quite different fin filling process observed near $E_l^{2i+1}$ on the concave left front suggest that more drastic shape transformation and kinetic energy dissipation stronger in thin films cause the assymmetry between average velocities of fronts progressing in x end − x directions.

Fig. 5 Time evolution of actual drop wicking on stainless steel structure B (middle structure). The bottom of the structure appears as white and the top is dark. Each structure is 24 mm long. Water is colored for better visibility. Times are in seconds. The first picture shows the pipette that places a water drop on structure B at time t0. At time t1 the drop has already begun to wick through the channels. At time t2, the structure has been completely flooded in x directionwhile wicking is still progressing in the opposite direction. Transport in − x direction is completed at time t3 only: slower wicking is observed in this direction.

6 Conclusion

Numerical simulation and experiment reveal that single drop wicking on flat metal surface equipped of fin (or arrow head) shaped parallel channels as in [9] induces a net wicking in the direction opposite to the one pointed by the arrow heads. They confirm more rapid motion on thinner channels, already revealed by [6]. We demonstrate that efficient directional wicking actually occurs on micro-sized steel patterns manufactured by SLM, paving the way to applications in industrial framework.

Numerical results retrieve the uneven dynamics with pulsing process. Moreover, the simulation provides access to the dynamics description at small space and time scale. In particular, it points out that the full non-linear Navier-Stokes equation has to be considered in contrast to the quasi-static model assumed in the literature. Nevertheless, the VOF simulation apparently exaggerates interface diffusion and underestimates the capillary forces that drive wicking on channel bottom: its issues agree with the experimental right front progression provided we insert smaller θ value. The LBM does not need such trick to give an idea of the causes of here described non-symmetric front dynamics. This motivates further efforts with this method, to be implemented in a parallel code in view of initial conditions and flow domains more similarto realistic experimental conditions. An experimental device allowing better time resolution is also desirable to asses such numerical attempt.

Author contribution statement

All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.

References

All Tables

All Figures

	Fig. 1 View from above of horizontal stainless steel surface modified by super-imposing fin (or arrow head) shaped parallel channels. Channel bottom/walls appear in withe/dark. We will see that water drops preferentially wick in x direction, the opposite to the one to which the arrows point.
In the text

	Fig. 2 Computer-aided design perspective view of the structure pictured on Figure 1. Observe the two directions of the faces of the channel walls.
In the text

	Fig. 3 The basis of the elementary cell of a channel in white. Channel walls are in blue. In our study lengths L1 = L2 and Wi take the values specified in Table 1, and φ = 27°.
In the text

Fig. 4 Simulated liquid velocity and free surface height during the wicking of a parallelepiped of liquid over a single fin shaped channel, using the VOF. The top of each sub-figure documents the velocity, and just below we see liquid height distribution at z = 0.05 mm above channel bottom. Contact angle is 50°. The simulation domain consists of 6 elementary cells. The initial liquid parallelepiped is one elementary cell long. Dashed lines mark its limits. Just after t = 0 the liquid parallelepiped deforms quickly to smoother shape and wicks on channel bottom in the two directions x and − x. At t = 0.34 ms it has invaded the straight channel at the right, and is still progressing in the wider part at the left. However, at t = 0.64 ms the left front has compensated its delay. At t = 1.56 ms the left front has clearly covered more distance than the right front. At t = 1.94 ms the right front spreads a little further than $E_r^3$ before immediately receding because of surface tension that tends to decrease the free surface. Hereafter the right front remains pinned. On the other side, differently, we see that the flow suddenly accelerates in the fins since the time the left front touches $E_l^3$ at t = 1.56 ms. Between time instants t = 1.56 ms and t = 1.94 ms, the acceleration is caused by the tapered geometry of the fin that increases the wet solid surface over liquid volume ratio. At t = 1.94 ms, the liquid completely wets the fin surface but the kinetic energy is only partly dissipated. This results into reverse flow in the fins, converging to the middle of the channel. The volume increase in the center of the channel increases the free surface area: this transforms a fraction of the kinetic energy into surface energy. Relaxation then causes water to spread again, and to flood the next straight channel which acts as capillary.

In the text

Fig. 5 Time evolution of actual drop wicking on stainless steel structure B (middle structure). The bottom of the structure appears as white and the top is dark. Each structure is 24 mm long. Water is colored for better visibility. Times are in seconds. The first picture shows the pipette that places a water drop on structure B at time t0. At time t1 the drop has already begun to wick through the channels. At time t2, the structure has been completely flooded in x directionwhile wicking is still progressing in the opposite direction. Transport in − x direction is completed at time t3 only: slower wicking is observed in this direction.

In the text