# Directional Water Wicking on a Metal Surface Patterned by Microchannels

Nima Abbaspour, Philippe Beltrame, Marie-Christine Néel, Volker Schulz

► To cite this version:

Nima Abbaspour, Philippe Beltrame, Marie-Christine Néel, Volker Schulz. Directional Water Wicking on a Metal Surface Patterned by Microchannels. Materials, 2021, 14 (3), pp.490. 10.3390/ma14030490. [hal-03769541](https://hal.inrae.fr/hal-03769541)

HAL Id: hal-03769541

https://hal.inrae.fr/hal-03769541

Submitted on 5 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Directional Water Wicking on a Metal Surface Patterned by Microchannels

Nima Abbaspour <fnref n="1" />*, Philippe Beltrame <fnref n="1" />*, Marie-Christine Néel <fnref n="1" /> and Volker P. Schulz <fnref n="2" />

**Abstract:** This work focuses on the simulation and experimental study of directional wicking of water on a surface structured by open microchannels. Stainless steel was chosen as the material for the structure motivated by industrial applications as fuel cells. Inspired by nature and literature, we designed a fin type structure. Using Selective Laser Melting (SLM) the fin type structure was manufactured additively with a resolution down to about 30 µm. The geometry was manufactured with three different scalings and both the experiments and the simulation show that the efficiency of the water transport depends on dimensionless numbers such as Reynolds and Capillary numbers. Full 3D numerical simulations of the multiphase Navier-Stokes equations using Volume of Fluid (VOF) and Lattice-Boltzmann (LBM) methods reproduce qualitatively the experimental results and provide new insight into the details of dynamics at small space and time scales. The influence of the static contact angle on the directional wicking was also studied. The simulation enabled estimation of the contact angle threshold beyond which transport vanishes in addition to the optimal contact angle for transport.

**Keywords:** directional wicking; patterned surface; wetting dynamics; capillarity, 3D simulation; Volume of Fluid; Lattice Boltzmann Method; Selective Laser Melting Manufacturing; microstructure

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Abbaspour, N.; Beltrame, B.; Néel, M.-C.; Schulz, V. Directional Water Wicking on a Metal Surface Patterned by Microchannels.
*Materials* **2021**, *14*, 490. https://doi.org/10.3390/ma14030490

Academic Editor: Hansang Kwon
Received: 4 November 2020
Accepted: 14 January 2021
Published: 20 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

In microfluidics, water forms droplets which remain pinned on solid surfaces. The transport or spread of these droplets is a challenging task in an industrial framework as lubrication, chemical reaction, heat exchanger, and water management in fuel cells. In most microfluidic devices, liquid motion is provided by a pump or an external energy input such as piezo elements. Nevertheless, surface heterogeneity as a gradient of wettability, i.e., the solid-liquid contact angle, allows droplet spreading without any external energy input [1,2]. It is a relaxation dynamics of the drop/surface system which decreases its free energy. Such liquid movement is usually called *passive transport* when the liquid flow has a privileged direction. Passive water transport at microscale offers many applications [3,4].

Such a directional transport also appears in the living world. Some species such as insects, lizards and plants living in deserts use it to harvest water [5–8]. Heterogeneous wettability driven liquid motion was also observed on wharf roach skin [9], and its frequent use by living organisms inspired a variety of industrial applications [10–12]. In all these examples, it is the surface heterogeneity that causes liquid drops to move. Heterogeneity can result from spatial variations in the wettability [13–15] but can also be caused by small scale (micrometric) topographic structures on the smooth solid surface acting as an asymmetric capillary [5–7,16–18]. The water spreading through these microchannels is called water wicking in the literature [16].

The initial motivation of this work is the use of such a passive transport to remove the liquid water excess in the Proton-Exchange Membrane (PEM) fuel cell [19]. Indeed,

<footnote n="1">¹ UMR1114 EMMAH INRAE—Avignon Université, F-84914 Avignon, France;
marie-christine.neel@univ-avignon.fr</footnote>

<footnote n="2">² Department of Mechanical Engineering, Baden-Württemberg Cooperative State University Mannheim,
Coblitzallee 1-9, D-68163 Mannheim, Germany; volker.schulz@dhbw-mannheim.de</footnote>

<footnote>* Correspondence: eng.abbaspour@gmail.com (N.A.); philippe.beltrame@univ-avignon.fr (P.B.)</footnote>
water management is a crucial issue as evidenced by several patents and literature on this topic over the last decade. Part of researches has focused on the design of the channels containing the reactant gas flow and water in the bipolar plate [20–22]. Another possibility is to implement additional micro-channels to drive the water by capillarity in the bipolar plate [23] or in the membrane [24]. A third possibility is the improvement of the shape of the channel cross section [25,26]. Passive transport by modifying the surface properties of the channel is an innovative approach that has not been studied in the framework of PEM Fuel Cell (PEMFC). Transport using gradient wettability requires a physicochemical treatment applied to the metal surface [27–29]. The drawback of this treatment is the possible modification of the metal properties as conductivity and the durability of this coating. Thus, we turned to the water wicking due to microstructure on the surfaces as in [30].

The realisation of such structures in the PEMFC context requires that chosen materials are compatible with low temperature hydrogen fuel cell constraints, i.e., robust and electrical conduction. These constraints raise open questions that have not been addressed in the literature. Indeed, experiments having reported passive directional water motion driven by microstructures [5–7,9] used other liquids than water and other materials than metal. Then, there is no evidence that reproducing the pattern found in literature will lead to water wicking because water on metal is less wettable than other couples liquid/substrate used in the previous references. Moreover, if the technology progress of additive manufacturing allows complex geometry designs in metal, especially in the field of fuel cells [31–33], none have been made, to our knowledge, at the scale of a few tens of microns.

The main goal of the present paper is to manufacture micro-sized pattern on a metal surface in order to obtain spontaneous directional water drop wicking. It is known that open microchannels are acting as micro-capillary [34] and spreading a drop even if the liquid is poorly wettable on the substrate. The directional spread is obtained only if the structure is asymmetric. We considered open capillaries constituted by repeated fin-shaped patterns such as those proposed by [16]. Drawing on the latter reference, we present an experiment that shows the directional spreading of a water drop on the patterned surface. The experiment was prepared and interpreted by numerical simulation of the evolution of the drop. Several numerical simulations have analyzed the evolution of velocity field, water content, interfaces and pressure in the latter case [1, 2, 13–15, 27–29, 35]. Theoretical approaches to drop wicking in channels of complex geometry [16,17] merely use semi-empirical formulas in the framework of quasi-static flow. The main idea of the directional behavior lies in the hysteresis phenomenon of a front liquid arriving at an edge: a threshold of the advancing contact angle should be overcome to further advance according to the Gibb's criterion [36]. The relevance of this criterion in our case is analysed in Section 5. However, none of these references have attempted to simulate the time evolution of liquid drops far from quasi-static regime in complex geometries such as the fin-shaped structures considered by [16] nor a periodic distribution of tilted pillars [18]. A presentation of the experiment and a first numerical analysis reported in [37] suggested that the fluid inertia cannot be neglected during the wicking. Simulating a free flow surface wetting dynamics with the complete Navier-Sotkes in a complex geometry is a challenging task. The simulation is performed using two different approaches: the Volume of Fluid method and the Lattice Boltzmann Method [38]. The advantages and the limits of each method is discussed in the paper.

The paper is organized as follows: in the next section, we specify the constraint due to the PEM fuel cell framework and the manufacturing (Section 2). Then, the results of the drop wicking experiments are reported in Section 3. After presenting the models used for the numerical simulations we analyze the simulated dynamics in Section 5. Finally, Section 6 summarises our main results and the future outlook.
## 2. Problem and Method

We applied numerical and experimental approaches to predict and evaluate spontaneous liquid displacement on manufactured structures made by micro laser sintering with stainless steel.

### 2.1. *Structure Geometry*

The structure of the micro-channels is based on the fin-shaped one on which [16] observed one-directional wicking. Figure 1 represents, left panel, the micro-channels considered here and, right panel, the original version of [16]. As can be seen on Figure 1, the studied channels are about four times wider and two times deeper than the original version of [16]. At larger scale the one-directional wicking is not a priori guaranteed since the experiments conducted in [17] demonstrated that pattern scale influences spontaneous directional wicking, which [16] found to be more rapid on thinner channels.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Fin shaped structures added to a flat surface. Top views of the fin-shaped structure studied here and manufactured by Selective Laser Melting (SLM) in a picture from above (**left**), CAD in an isometric view (**middle**), and original Polydimethylsiloxane (PDMS) part (**right**). Notice the different scales and compare the rough metallic structure manufactured by SLM (**left**) with the smooth surface considered in [16] made of resin (**right**).</figcaption></figure>

The fuel cell framework also imposes the use of a material that significantly differs from those used in [39]. These changes are discussed below.

### 2.2. *Material Constraints*

A robust material that conducts electric current and resists corrosion and heat is highly desirable for our purpose. Stainless steel satisfies these requirements. Since [16] found smaller patterns more efficient, we manufacture our structures by the Selective Laser Melting (SLM) which is adapted to steel and achieves high resolution. Nevertheless, sizes as small as in [16] are out of reach. Moreover, the SLM technique cannot avoid slight deviations between CAD design and manufactured component, which we also expect to influence wicking.

Finally, the reaction product in a PEM-fuel cell is pure water, which returns a water/steel contact angle larger than in the drop transport experiment [16,39]. The latter was performed on Polymethylmethacrylate (PMMA) and Polydimethylsiloxane (PDMS) substrates, which of course do not satisfy our robustness and electric conductivity requirements. Moreover, reference [39] associated these substrates with IPA and soapy water. Since it seems that drop wicking on a given structure is decreased when the contact angle is increased, we expect to find the optimal conditions upon exploring the 0–70 degrees range. Numerical simulation is the most flexible tool to study the influence of contact angle. However, it needs to be validated by experiments. Our method consists in combining these two approaches.

### 2.3. *Method*

A numerical approach using either interface tracking by Volume of Fluid (VOF) or the Lattice Boltzmann Method (LBM), captures drop motion in a channel limited by the patterns described in Figure 1. We will see in Section 5 that interface tracking suggests directional drop motion at contact angles smaller than 70 degrees. The numerical approach was
validated by an experimental set-up on steel structures based on the geometry introduced by [16]. Liquid transport was observed immediately after a distilled water drop had been placed on each one of three steel structures manufactured in-house and exhibiting the fin-shaped geometry of Figure 2 with slightly varied scales. It turns out that channel width influences spontaneous drop displacement. Manufacturing and measurements are detailed in Section 3. After a brief description of the underlying physical model in Section 4, numerical issues are detailed in Section 5 which also compares the simulation results with the experiment.

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Three adjacent structures exhibiting fin shaped channels added to the steel base by SLM. (left) Actual A, B and C structures. (right) CAD design documenting the total length and width of each structure, in millimeters.</figcaption></figure>

**3. Experimental Study of Liquid Wicking on Fin-Shaped Patterns**

The manufacturing process of adding fin-shaped patterns of various widths on the steel surface is the first step of our experiment. The next step consists in observing drop evolution on the structure thus obtained.

3.1. *Substrate Manufacturing and Physical Properties*

As in [16] the substrates used here to study drop wicking are flat surfaces with fin
patterns forming channels. One of them is shown in Figure 1. The platform and patterns in
the present study were made of 1.4112 stainless steel, and three such structures denoted A,
B and C were built on a single platform (Figure 2). For manufacturing first the CAD data
is provided (Figure 1 middle panel). In second step, the CAD data is converted to layers.
In the third step of manufacturing by Selective Laser Melting (SLM) the machine deposits
an uniform thin layer of metal powder on the surface and melts each layer successively
onx the platform by a laser beam. The past two steps of spreading the uniform powder
layer and sintering are repeated until the channel wall height reaches the desired value
which here is 100 [µm] (shown in Figure 2 right panel). Then, the powder residue must be
cleaned and the manufacturing is completed. A single cell of the geometries from structure
A and B are depicted in Figure 3 left and right respectively. The height of the structure
is visualized by a 3D microscopic scan. The deviation from CAD design of the structure
A is relatively larger than for the C structure. For the structure A, the roughness is not
negligible compared to the size of the structure.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The pictures represent the basic element of structure manufactured by SLM. (left) Structure A and (right) structure C of the experimental set-up (Figure 2) corresponding to the smallest and largest lengths specified in Figure 4. Colours indicate height measured from the steel platform.</figcaption></figure>

This process results in a series of parallel capillaries. Each capillary is a periodic repetition in the longitudinal x direction of an elementary pattern shown Figure 4 left panel. The blue region in the rectangle is the melted metal of a height 100 µm while the white region is an empty cavity bounded by the walls of the blue region. The cavity constitutes the elementary cell of the capillary. Note, that the parallel capillaries are separated by the walls of a width at least larger than W3. Thus, water can only flood from one capillary to the other by passing over the wall. The geometric characteristics of the fin shape is given in the table of Figure 4 right panel. In the elementary cell, we distinguish the straight channel, called *channel*, from the wider region constituted by the fins, called *V* region. Figure 2 right panel shows the top view of the CAD drawing of the three sets of capillaries called structures A, B and C. Their respective widths are 5, 5 and 6 mm. Diameter and height of the platform are 50 and 10 mm respectively. Left panel of Figure 2 is a photo of the round platform and three adjacent structures manufactured by SLM at the center.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** A single cell of the pattern from which is deduced the capillaries. The blue colored part is the fused melted metal with a height of 0.1 mm. The bottom of the capillary colored white. Both parts create an elementary cell which is included in the rectangle 'abcd'. Adding copies of this rectangle adjacent to its parallel to y sides defines one channel and its walls. Then, adding adjacent copies of the ensemble forms the capillaries for structure A, B or C.</figcaption></figure>

<table><thead><tr><th>Variable</th><th>A</th><th>B</th><th>C</th></tr></thead><tbody><tr><td>φ</td><td>27°</td><td>27°</td><td>27°</td></tr><tr><td>L1=L2[mm]</td><td>0.2</td><td>0.64</td><td>0.4</td></tr><tr><td>W1[mm]</td><td>0.1</td><td>0.16</td><td>0.1</td></tr><tr><td>W2[mm]</td><td>0.3</td><td>0.8</td><td>0.5</td></tr><tr><td>W3[mm]</td><td>0.1</td><td>0.16</td><td>0.1</td></tr><tr><td>depth[mm]</td><td>0.1</td><td>0.1</td><td>0.1</td></tr></tbody></table>

Geometric specification is not sufficient to determine the functionality of the structure to promote directional wicking since initial surface properties and powder grain size influence wetting dynamics of each structure. A slight initial platform roughness is necessary to enhance powder adhesion. Here, the initial platform roughness was 207 RMS. However, increasing roughness may block the wetting front. Another influential parameter is the surface roughness of the laser sintered part. On one hand, we used steel powder with a grain diameter of about 5 µm. On another hand, the resolution of the final product, it is also
strongly influenced by the laser beam focus. Here the final resolution was 30 µm, and it can be seen on Figures 1 and 3 that the channel walls deviate slightly from the isometric CAD drawing of Figure 1 (middle) and exhibit small scale defects. In addition, the steel pattern surface (Figure 1 (left)) is rough compared to the resin surface (Figure 1 (right)) of [16].

We thus obtained three manufactured structures A, B and C positioned side by side on the same steel plate. Each structure is wide enough for a drop to be placed on it without touching its boundaries. A picture and a schematic of the resulting device are shown Figure 2.

## 3.2. Experiment

The spontaneous drop motion on A, B and C is captured by videos of three liquid drops successively placed on these structures. For this purpose, we used a Keyence VHX-600 microscope equipped with a VH-Z20R and VH-Z100R lens. The platform was placed horizontally under microscope during the experiment. In addition, photos and videos were captured by Samsung Galaxy Note 8 and Olympus Stylus SH-50. In all the experiments the distilled water has a contact angle of about 50 degrees.

Qualitatively, the wicking dynamics is similar on all structures. Especially, wicking starts as soon as the drop is placed on it at initial time $t_0$. Figure 5 displays four snapshots representing the different steps ($t_0, t_1, t_2, t_3$) of the evolution of a single drop placed in the middle of structure B. The snapshots were extracted from the video. The video can be found in the experiment subsection of Supplementary Materials. At the beginning, the water fills only the capillaries where the drop is initially setting, i.e., the three center capillaries. The mostleft and mostright capillaries are progressively flooded in a second step.

The wicking is faster in the $+x$ than in the $-x$ direction. The mean velocity of wicking in each direction is estimated by recording the initial positions of the front ($x_i^+$) and the back ($x_i^-$) of the drop at $t_0$ and by measuring the time taken to reach the ends ($x_f^+$) of the structure (Figure 5). With these notations, the average velocity in each direction is

$$ \langle V^{\pm} \rangle = \frac{x_f^{\pm} - x_i^{\pm}}{t_f^{\pm} - t_0}. \qquad (1) $$

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Time evolution of directional drop wicking on structure B. The bottom of the channel appears as white and the top of the structure is dark. Water is colored for better visibility. The first picture at time $t_0$ present the pipette that places a water drop on B structure. At time $t_1$ the drop has begun to wick through the channels. At time $t_2$, the structure is completely flooded in the positive direction while wicking progress is slower in the negative direction. Finally at time $t_3$ transport in the negative direction is completed.</figcaption></figure>
Figure 6 documents absolute values observed in these two directions on the three structures A, B, and C. On each of these three structures the average velocity in the negative direction is nearly half the value observed in the positive direction. According to the capillary sizes shown in the table of Figure 4, we observe larger velocities in both directions x and −x on thinner structures.

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Average velocity absolute values of water wicking in x and −x directions for the three different channels A, B and C.</figcaption></figure>

To better understand, the wicking dynamics in both directions, in the following paragraph, snapshots of the front evolution at the cell scale are displayed.

In the negative direction, we were able to capture slow down of the front dynamics during the experiments of Figure 5. In the first snapshot ($t = 0$ [ms]) of Figure 7 the front is blocked at the right ($-x$ direction) of the first cell in the thin straight channel. However, the second cell is not completely dry: a small amount of water can be observed in the tip of the fin (see red arrows in the first panel of Figure 7). The fin volume is progressively filled up by water, first slowly till $t = 12$ [ms] and then faster until $t = 24$ [ms]. At $t = 30$ [ms] the water continues into the next cell in a same way as for the snapshot at $t = 0$ [ms]. Therefore, the front remains blocked at the entry of the second cell but the water finds another path to advance to the second cell. The reason of the existence of the thin water film connecting the blocked front is unknown. One possibility is that water is flowing on the upper surface of the pattern. However, the video cannot corroborate or infirm this assumption.

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Evolution of water transport in channel type *B* in the negative direction. Experimental conditions are the same as in Figure 5. Water is colored in blue for better visibility and the triple line is highlighted by a bold blue line. Both red arrows for the panel at $t = 0$ [ms] indicates the presence of water in the tip of the fin.</figcaption></figure>

Due to fast progress of the liquid in positive direction, we were not able to have similar snapshots for the front evolution. However, during the first experimental tests, the front progress was much slower than reported in Figure 6. We believe that chemical impurities on the surface are responsible for pinning effect as reported in [40]. Figure 8 shows the snapshots front evolution over two cells of the fine structure A in this context. The triple
line displays a concave line as described in the literature. Between the time $t = 1.6$ [s] and 3.8 [s], the front is asymmetric and remains almost static. We argue that the chemical and topographic defects of the melted metal surfaces visible in the first panel are responsible for this pinning effect. Once the front is depinned, the curve retrieves the $y \leftrightarrow -y$ parity symmetry. At $t = 9.2$ [s] the front is near touching both acute edges highlighted with red dots. Below the camera time resolution, the front passes through the straight channel and the water floods the second cell. A similar scenario occurs for the second cell. Note that this time the front breaks the symmetry in $y \leftrightarrow -y$ before it touches both acute edges however, that does not affect the sudden front burst. Therefore, there is a high velocity contrast between the flooding of the V region and the fast evolution by touching the acute edges.

The numerical simulation in Section 5 is performed by neglecting the surface defects but reveals an analogous front burst due to the acute edges. Therefore, the snapshots of Figure 8 can help to understand the liquid progress in a clean surface as well.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Evolution of water transport in channel type A in the positive direction. Experimental conditions are the same as in Figure 5. Many chemical defects slow down the dynamics. Water is delimited by a blue line for better visibility. Red dots of the snapshot at $t = 9.2$ [ms] highlight the acute edges discussed in the text.</figcaption></figure>

We discuss the different steps of the dynamics in Section 5, with the help of the numerical simulation that returns a much smaller time resolution.

## 4. Physical Model and Numerical Methods

The wicking dynamics is governed by the liquid properties $\rho, \mu, \sigma$ (density, viscosity and surface tension) and by the interaction between fluid and metal surface through the static contact angle $\theta$. We employ two methods to numerically simulate the two-phase flow: the liquid and the gas. Most of simulations are performed using the Volume of Fluid (VOF) method [41–43]. The temporal evolution of the velocity field in each fluid phase is governed by the Navier-Stokes equation in viscous and incompressible fluid. At solid boundaries we consider wall adhesion. The liquid/gas interface is modelled by a region of non-zero thickness where one passes continuously from one phase of the fluid to the other. Finally, at the triple air/liquid/solid interface there is the so-called triple line. At equilibrium, the static contact angle $\theta$ characterises the wettability. Out of equilibrium, the macroscopic contact angle depends on the front dynamics and it differs from the static contact angle [44]. However, at the mesoscopic scale, it is still relevant to consider the contact angle $\theta$ [45]. Therefore, following [42] we assume that the contact angle during the front dynamics is equal to the static contact angle $\theta$. Equations of the VOF method are detailed in the present framework in [37].
However, the main drawback of this method is the interface diffusion which may underestimate the driving forces at the interface. Then, when the interface diffusion becomes large, we employ the color gradient Lattice Boltzmann Method (LBM). The equations of the LBM code is described in [46]. The relevance of each code depending on the simulation is discussed in Section 5. Usually, generating structured mesh for VOF would be more time consuming than for the LBM method. However, for compensation we have applied symmetry to tetrahedral mesh for 2D case to increase the computational speed and avoid asymmetric behaviors due to unstructured mesh. We have not used the symmetry for 3D simulation. Then, the 3D simulation requires more CPU time than LBM method. The type of mesh and the number of cells is summarized in Table 1.

<table><caption>**Table 1.** Mesh properties of the Volume of Fluid (VOF) and Lattice-Boltzmann (LBM) methods used for the simulations of Section 5.</caption><thead><tr><th></th><th>VOF 2D</th><th>VOF 3D</th><th>LBM 2D</th><th>LBM 3D</th></tr></thead><tbody><tr><td>Mesh type</td><td>Tetrahedral</td><td>Hex dominant</td><td>Pure Hex</td><td>Pure Hex</td></tr><tr><td>Cells number</td><td>3904</td><td>71,899–758,198</td><td>26,180</td><td>551,368</td></tr></tbody></table>

In order to determine the main parameter of the dynamics, we estimate the following dimensionless numbers: Reynolds number (Re), Capillary number (Ca), Weber number (We) and Bond number (Bo) (Table 2).

<table><caption>**Table 2.** The dimensionless numbers and their magnitude. V characteristic velocity, W1 minimum channel width and $l_c = 2.7$ mm the capillary length of water.</caption><thead><tr><th>Dimensionless Number</th><th></th><th>Magnitude</th></tr></thead><tbody><tr><td>Reynolds number</td><td>Re = <sup>ρ<sub>V</sub>W1</sup>&frasl;<sub>μ</sub></td><td>0–100</td></tr><tr><td>Capillary number</td><td>Ca = <sup>μV</sup>&frasl;<sub>σ</sub></td><td>≈1</td></tr><tr><td>Weber number</td><td>We = <sup>ρV<sup>2</sup>W1</sup>&frasl;<sub>σ</sub></td><td>0–1</td></tr><tr><td>Bond number</td><td>Bo = <sup>W1</sup>&frasl;<sub>l<sub>c</sub></sub></td><td>≈1</td></tr></tbody></table>

In contrast to the quasi-static assumptions found in literature, the simulation shows that the flow velocity is about 1 [m/s] and the Reynolds number is about 100 during fast phases of the transport. Then, the flow is far from equilibrium and the inertia plays an important role. During the fast flow phase, the surface tension is in competition with the inertia force (We ≈ 1). Otherwise, the surface dominates over the inertia force (We ≪ 1). Moreover, since Ca ≪ 1 and Bo ≪ 1 the surface tension dominates over the viscous force and the gravity force. Hence, we neglect the gravity force as in [2].

Because of the computational cost, the numerical simulation is performed in a single capillary composed of 5 to 8 elementary cells repeated in the x direction of the geometries A, B and C. In a first step, we consider a 2D problem. The 2D simulation is obtained by considering the vertical z direction as the invariant direction. The most significant results are in 3D, but 2D simulation helps to understand the key role of the edges exhibited by the fin-shaped geometry.

In all simulations, we assume fixed and constant atmospheric pressure of the gas phase.

The water is introduced in the structure in two different manners: First one considers an initial volume of liquid and the second one considers a constant liquid flux throughout the simulation. In most 3D simulations, as for the experiment a small liquid volume is placed on the channel at the time instant $t = 0$. The initial liquid volume is a cuboid 0.3 mm high and the same W2 width as the channel and the volume in the channel below the cuboid is also filled up. This total volume is much smaller than the droplet volume in the experience. To avoid the limited liquid volume, a liquid flux in middle of the channel at a constant rate of 1.06 mm³/s is imposed all along the simulation. This condition is applied to 2D simulation too. In this case, the liquid is injected at two inlets
positioned in the middle of the channel boundary of one elementary cell as indicated in Figure 9 with the same flux rate.

<figure><img src="image_12.png" /><figcaption>**Figure 9.** VOF two-dimensional simulation with $\theta = 50^\circ$ and the geometry of the horizontal section of channel A.</figcaption></figure>

The simulation provides a description of the wicking dynamics at smaller space and time resolution than the experiment. Moreover, simulation allows to explore a wide range of physical parameters and providing the optimal parameters that enhance the directional transport.

## 5. Water Wicking Simulation

When compared with our experimental conditions, the 2D simulation amplifies the role of the edges of the fin-shaped channel walls and retrieves (see Section 5.1) the selective pinning predicted by the Gibbs criterion near obstacles exhibiting edges [47]. In three-dimensional geometry, however, the role of the edges is mitigated by the channel bottom and the VOF results come closer to our experimental observations.

### 5.1. 2D Simulation

The 2D variant of the numerical code described in Section 4 represents a three-dimensional two-phase flow invariant in the $z$ direction, between infinite walls based on the boundary in the $(x, y)$ plane of the elementary cell defined in Figure 4. Such a flow is not expected to resemble our experiment closely, and especially its initial condition that depends on $z$. Instead, we impose a liquid flux (10.6 [mm$^2$/s]) at inlets located on both parallel sides of the domain boundary as on Figure 9. These conditions are reminiscent of reference [47] that describes a ($z$ dependent) liquid flooding a flat surface decorated with obstacles that exhibit wedges.

In our simulation we observe that the liquid, driven by the capillary force, begins by completely filling the straight narrow channel that includes the inlet (Figure 10a). The left and right fronts are pinned at edges called $E_r$ and $E_\ell$, respectively located at the right and at the left of the inlet. To move again, a driving force is required. In our context, it is the pressure due to the continuous flow at the inlet. Since $E_r$ is much sharper than $E_\ell$, the left front rapidly depins whereas at $E_r$ the interface bulges gradually while the triple line remains pinned (see Figure 10a-c). The right $E_r$ edge is thus responsible for front pinning while the edge at $E_\ell$ is not: this is the key of the directional transport in 2D geometry. Note that the acute wedges at the end of the fins do not pin the front moving to the left since they are wetted by water on both sides (Figure 10c,d). References [2,28] already suggested this quasi-static explanation for the directional transport, based on contact angle hysteresis. At $t = 3.62$ [ms] (Figure 10), the left front is similar to that of the experiment at $t = 9.2$ [s] (Figure 8). However, the 2D simulation does not accurately describe the flow in the fins, as an air volume is trapped in the fins, blocking water progression (see Figure 10d), a phenomenon that we do not observe either in the experiment or in the 3D simulation.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Snapshots of the asymmetric transport using the 2D VOF simulation. Parameters as in Figure 9.</figcaption></figure>

## 5.2. 3D Simulation at 50° Contact Angle

In the 3D simulations the geometry of one cell is identical to that used in the experiment. However, one focuses on a unique series of lined up cells and periodic boundary conditions are imposed in the lateral y direction. An initial condition similar to the experiment consists of a rectangular ‘droplet’ placed on the structure (Figure 11). The case of water introduced by an inlet at the bottom is studied at the end of Section 5.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** 3D simulation of liquid wicking driven by capillary forces on channel A. Blue color represents liquid volume. Contact angle is 50°. The simulation domain is eight elementary cells long. The parallelepiped size is $l_x = 0.7$ mm, $l_y = 0.38$ mm, and $l_z = 0.4$ mm.</figcaption></figure>

In these two cases, the 3D simulated flow remains symmetric regarding the $y \leftrightarrow -y$ in contrast to the experiment (Figure 8). This corroborates the hypothesis made in the Section 3 that asymmetry comes from local defects and not from a spontaneous symmetry-breaking of the flow. Moreover, 3D simulation mimics the intermittent behavior suggested by the videos at the end of Section 3, alternating fast and quasi-static steps. The smaller resolution of the numerical approach uncovers the details of these different phases which we now describe in the case of a 50° contact angle.

Just after having initially been placed on the structure (at $t = 0$) the liquid parallelepiped deforms quickly to a smoother shape and begins to wick on the channel bottom in the two directions $x$ and $-x$. This can be seen on Figure 11 that shows four snapshots of the simulated liquid volume: at time $t = 2.5$ [ms] it has flooded one cell to the right of its initial position, and two cells to the left. However, the right front remains pinned until $t = 4.4$ [ms] while the left front floods one cell further and proceeds in the $x$ direction. After $t = 4.4$ [ms] the right front floods one cell more to the right before stopping. Meanwhile, the left front continues proceeding to the left until it reaches, at $t = 8.6$ [ms], the end of the computational domain.

Figure 12 investigates the details of the above described pinning and fin filling processes in the $-x$ and $x$ directions. Figure 12 shows snapshots of a simulation started from an initial parallelepiped shorter than in Figure 11, and slightly longer than the elementary cell. The figure is divided into four panels captured at different times. Each panel documents liquid height below, and the velocity field at $z = 0.05$ [mm] above. Since the two-dimensional approach and the discussion of Figure 11 pointed out the crucial role of convex edges, symbols $E_t^i$ and $E_l^i$ ($i = 1, 2, ...$) mark the channel cross-sections that include such elements to the right and to the left of the initial liquid volume. At $t = 0.34$ [ms], shortly after the initial condition, the first cell is filled and the water spreads to the left
and to the right. The view from above of the front evolution (especially the concave left front) is reminiscent of the two-dimensional simulation sketched on Figure 10. However, this time the convex front at the right stays pinned on $E_r^1$ for a finite time after which the liquid spreads out further to the right into the wider V part of the channel (Figure 12 at $t = 0.64$ [ms]). Two main reasons explain this discrepancy with Section 5.1. The first one is that a part of the triple line now lies on the channel bottom while contact angle hysteresis only occurs at the edges. Advancing by spreading on the bottom surface and circumventing edges is easier, and it is what we observe at $t = 0.64$ [ms]. The second reason is the velocity field at $t = 0.34$ [ms] which suggests that inertia should not be neglected, differently from the quasi-static approach of [16].

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Simulated liquid velocity and free surface height during the wicking of a parallelepiped of liquid over a single fin shaped channel. For each sub-figure, the color scale displays the norm of the velocity field in the upper panel and the liquid height distribution for the bottom panel. Contact angle is 50°. The simulation domain consists of 6 elementary cells. The initial liquid parallelepiped is one elementary cell long. Figure credit from [37] Figure 4.</figcaption></figure>

Once the fins between $E_r^1$ and $E_r^2$ are filled (shortly after $t = 0.64$ [ms]) the right front continues to advance in the straight narrow channel. At $t = 1.56$ [ms] it is at $E_r^3$. However, this time the flow velocity is about two times slower. Therefore, the right front spreads a little beyond $E_r^3$ before immediately receding ($t = 1.94$ [ms] in Figure 12). Surface tension causes the back flow by decreasing the free surface, and thereafter the front remains pinned. Kinetic energy is now dissipated.

On the opposite side, the left front approaches $E_l^2$ at $t = 0.64$ [ms]. When it touches the acute wedge $E_l^3$ at $t = 1.56$ [ms], the flow suddenly accelerates as observed for the experiment (Figure 8). At the time instant $t = 1.94$ [ms] the velocity is about 1 [ms] in the fin volume. The main explanation of this acceleration is that the tapered geometry of the fin increases the ratio of the wet solid surface to the liquid volume in the fin implying the increase of the ratio of the driving capillary force over the inertia force. During the fin flooding, the kinetic energy is only partly dissipated when the liquid completely wets the fin surface. This results in a reverse flow in the fins, converging to the middle of the channel between $E_l^2$ and $E_l^3$ (Figure 12 at $t = 1.94$ [ms]). This back flow increases the free surface area of the liquid. Thus a fraction of the kinetic energy is transformed into surface
energy. This potential energy is added to the capillary effect due to the straight channel between $E_l^3$ and $E_l^4$. A fast flooding then occurs in this channel.

The same trend is observed on Figure 13 that details the fin filling process and represents the free surface and its velocity from a simulation in which the water input was different, however the figure specifies the behavior of the front near a cross section that it reaches at $t = 1.3$ [ms], and plays the role of $E_l^3$. Observe the velocity increase that develops in the fins before returning to the median region. Meanwhile the front advances slightly in the straight channel before receding during a pause that lasts over 0.8 [ms]. However, after this pause at $t = 2.1$ [ms] the velocity field is non-zero in the median region. The relaxation of the free surface is responsible for this significant velocity field. Then, this kinetic energy helps the front to flood the next straight channel to the left. A video presenting the evolution presented here can be found in the simulation subsection of Supplementary Materials in the experimental part.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** VOF simulation detailing fin filling in the x direction. The continuous color represents the liquid free surface. Arrows display the velocity on the liquid free surface with a color code for the amplitude. However, instead of being started from an initial condition as for Figures 11 and 12, the simulation is issued from the continuous fluid inlet condition described in Section 5.4. The time is set to zero when the liquid enters the V region.</figcaption></figure>

Numerical simulation provides a better understanding of the dynamics during the fast phase that video of the experiment cannot capture. In particular, the interplay between inertia and surface tension explains that the dynamics remains fast after the filling of the fin region.

The simulation differs from experiment during the slower phase when the front advances in the wider part of V. The simulated front dynamics is faster, regular and symmetric ($y \leftrightarrow -y$) contrary to the stick-slip motion and to the asymmetric of the shape front described in Section 3. This proves the role of surface defects which are responsible for the pinning effect that slows significantly down the dynamics. In addition, in the simulation the left directional transport is limited by the amount of liquid deposited and we expect that the meniscus will stop before reaching the end of the channel in a longer device with the initial condition considered here. Likewise, the simulated meniscus on the right of the fluid volume is pinned at $E_l^3$ whereas in the experiments it continues unevenly to the right. A possible cause is the initial amount of liquid, which is about 50 mm³ in the
experiment whereas it is only 0.1 mm³ in the simulation. An attempt to circumvent the problem of initial volume is studied in Section 5.4.

The number of cells filled by the water also depends on the static contact angle between liquid and the structure and on pattern size A, B or C. The next section investigates the influence of the contact angle and the pattern size.

### 5.3. Influence of the Contact Angle and Geometry

In the first part, only the static contact angle $θ$ varies but all other parameters remain the same and the pattern size is that of structure A. In other words, we simulate a change of substrate wettability.

We observe that increasing the static contact angle $θ$ decreases the mean front velocity in the positive direction $x$ and the right front remains pinned. This behaviour can be explained by the decrease of capillary forces as the contact angle increases. Beyond $θ_{lim} \simeq 70°$, the wicking through the capillaries does no longer occur. The rectangular water volume converges rapidly to an equilibrium state which is similar to the steady-ridge on a flat homogeneous substrate described in [40] and with periodic boundary conditions in $y$ direction.

Conversely, decreasing $θ$ from $50°$ speeds up the wicking front progression in the $x$ direction (Figure 14). However, this affects the front pinning at $E_r^3$ in the $-x$ direction. For instance, at $θ = 30°$, the fluid spreads over a larger area before retracting and remaining pinned at $E_r^3$ (Figure 15). This dynamics results from a competition between inertia and surface tension. There is a critical value, noted $θ_c < 30°$, below which the front is no longer pinned at $E_r^3$ and water spreads in both directions. Even in this case the wicking is slower in the $-x$ than in the $x$ direction. Though the front is no longer pinned, the $E_r^3$ edges constrain the flow to spread out in a thin film beyond $E_r^3$. Since viscous dissipation is strongly related to the height of liquid films, the kinetic energy is almost totally dissipated passing through $E_r^3$ in the $-x$ direction as a thin film. Thus, the flooding in the $V$-region is slow and this explains the slower wicking than in the $x$ direction.

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Average velocity of the front as a function of the static contact angle $θ$. Blue crosses indicate front velocity in the $+x$ direction while red dots indicate the front velocity $-x$. Beyond $θ \ge 30°$ the front in the $-x$ direction is pinned and then velocity is not indicated. The geometry corresponds to the smallest pattern A of Figure 4.</figcaption></figure>

To quantitatively study the transport efficiency, we estimated the transport speed averaged on the few cells that were filled in the $x$ direction (Figure 14). The mean velocity in both directions decreases as the contact angle increases. The decrease in the $x$ direction is not linear. Between 0 and $10°$ the velocity is divided by 1.6 while in the larger range [$10°, 50°$] it is divided by 2. For $θ = 0°$ or $10°$ the asymmetry of the wicking is weak, as indicated by the similar values of velocity in the $x$ and $-x$ directions. Thus, the contact angle $θ_c \simeq 30°$
represents an optimal value of the directional wicking since it is the maximum velocity of
the front in the x direction for which the front is still blocked in –x direction.

<figure><img src="image_18.png" /><figcaption>**Figure 15.** Spreading of water beyond the pinning point (t = 1.17 ms) and backflow (t = 1.70 ms) for a static contact angle θ = 30°. Velocities and liquid height are documented at the top and bottom of each figure, respectively. The initial condition is a two elementary cells long liquid parallelepiped. The channel geometry corresponds to the pattern A.</figcaption></figure>

We also simulated the influence of the other geometries B and C of the experiment. The geometry does not qualitatively change the scenario described above. The three structures exhibit a similar trend for the contact angle dependency of the average velocity, showing in all cases an upper limit $\theta_{lim}$ for the wicking in the positive direction and a critical value $\theta_c$ below which the pinning effect vanishes in the negative direction.

The velocity of fluid wicking is influenced by pattern size (Figure 16). This is especially visible on structures A and C, for which the velocity difference is maximum for θ = 30°. This result is qualitatively very similar to the experiment (see Figure 6). Indeed the capillary force over the inertia ratio is larger in smaller cells: dividing the sizes in the (x, y) plane by factor f divides the water mass by $f^2$ whereas the capillary force, proportional to the contact line length, is divided by f.

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Influence of the pattern geometry on the average velocity of the front in the positive x direction. The contact angle is θ = 30°.</figcaption></figure>

Despite the qualitative agreement with experiment, the average velocity of the liquid spreading through the capillary is about ten times that of the experiment. Surface defects as shown in Figure 3 are known to slow down the wetting dynamics [40]. Even considering a smooth surface, the low volume of water in the simulation compared to the experiment is another possible reason for this discrepancy. Firstly, the low volume does not make it possible to average out over a large number of cells implying a large uncertainty. Secondly the Laplace pressure in a small droplet is higher than in a larger one because the Laplace
pressure is proportional to curvature. In the next section, we propose a continuous flow to
explore the dynamics without the constraint of a limited volume.

5.4. *Continuous Flow*

To circumvent the problem due to the low volume of liquid in the simulation, we impose a continuous inlet condition at the bottom of the straight channel in a similar way to the 2D simulation. An interesting feature is that this boundary condition is reminiscent of the context of the bipolar plates of a fuel cell as the reaction at the cathode produces a continuous water flow. We simulated the water production by constant flux (1.06 [mm³/s]) through a rectangular surface positioned at the bottom of the channel (Figure 17).

The simulated dynamics displays an uneven wicking progression in the positive
direction as described in the previous sections. This time, the spread in the positive
direction is no longer limited. In contrast, the front is still definitively pinned in the –x
direction at the $E_r^3$ edge. That proves out that the cell geometry is effective in obtaining a
directional wicking.

<figure><img src="image_20.png" /><figcaption>**Figure 17.** VOF simulation with an inlet at the bottom of channel A. Inlet flux is 1.06 mm³/s and $\theta = 50^\circ$,</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 18.** (a) Diffusion at the water/air interface in VOF simulation with compressive scheme. Parameters are as in Figure 11. (b) LBM simulation performs clear phase separation without excess diffusion at the interface with an equal number of cell elements. Red arrows indicate the location of the inlet. Contact angle $\theta = 50^\circ$.</figcaption></figure>

Nevertheless, the long duration of VOF simulation increases the numerical interface
diffusion. Figure 18a displays the large diffusion of the interface between water and air.
The overestimate transition region at the water-air interface impacts the modeling of the
driving Laplace force and the capillary force at the triple line. To decrease this numerical
effect, we used Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM)
that proved to be a better choice than the Geometric Reconstruction Scheme (GRS) or
High Resolution Interface Capturing (HRIC). To reduce diffusion at the water/air interface
local mesh refinement can be used which adds to the total computation time by mesh
generating for each specified timestep. Indeed, the main drawback of the VOF simulation
is that it exagerates the numerical interface diffusion [48]. We believe that the interface
diffusion may underestimate the capillary force and then mislead the wicking dynamics simulation. The color gradient Lattice Boltzmann Method (LBM) used in [46] provides an alternative to model interface dynamics. The interface diffusion is drastically reduced by LBM simulation compared to VOF for a similar mesh refinement (Figure 18). Using similar inlet conditions, the LBM simulation mimics the above described fast/slow alternating steps without definitively pinning the right front: the latter stops at the entrance of the fins but continues its progression after some time (Figure 19). The steps of the front progression in the negative direction (Figure 7) are retrieved in the LBM simulation. Even if the front is pinned at the sharp corner $E_r^i$, a small amount of water in the tip of the fin is present. The volume in this region grows and progressively fills the volume of fins. Then, the water continues to the next channel. Depinning does not take place but the water bypasses the $E_r^i$ edges. The details of the LBM simulation suggest that a thin film develops on the surface of the structure and opens a path to the fin. Nevertheless, according to [49], it is possibly an artifact of the LBM simulation due to instability at the interface between the liquid and gas phases.

<figure><img src="image_22.png" /><figcaption>**Figure 19.** LBM simulation with an inlet (red arrows) at the bottom of channel A. Contact angle $\theta = 50^\circ$.</figcaption></figure>

However, the similarity with the experiment leads us to infer that such thin liquid films actually exist. An explanation for these thin films could be the presence of surface irregularities at the microscale made by SLM manufacturing (Figure 3). The interstices between defects act as capillaries into which water may infiltrate and bypass the corner $E_r^i$. Therefore, the defects decrease the dissymmetry of water wicking by, on one hand, allowing a wicking in the negative direction and another hand, by slowing down the front (stick-slip motion) in the positive direction.

Finally, according to the VOF and LBM simulations the wicking in the negative direction occurs only if the water finds a path over the structure or a path due to the defects. Thus, the depth of the structure should adapted to the liquid volume to transport in order to avoid overflooding above the structures.

## 6. Conclusions

Motivated by the problem of water drainage in a fuel cell, our study demonstrated the feasibility of the directional spreading of a water drop on a structured steel substrate. The fin geometry of the micro-channels was based on literature, in particular [16]. However, we employed a conductive substrate in contrast to the non-conductive PDMS polymer of [16] and to organic substrates [5,7,8,10,39]. Because of the low wettability of water on metal, there was no evidence that the devices reported in the literature could apply in the current framework. Moreover, the manufacturing of channels with complex small-scale geometry was a challenging problem. To our knowledge, it is the first time that this type of structure has been manufactured additively with a resolution down to about 30 µm. Experiments using three different pattern sizes provide a wicking in a preferential direction. Its efficiency is correlated to the structure scale and this highlights the importance of making small structures.
For an in-depth understanding of the mechanisms involved in the liquid wicking we used a three-dimensional numerical simulation of multiphase Navier-Stokes equations taking into account free surface, contact angle and inertia. This approach contrasts with the quasi-static description in the literature [5,16–18]. We used two complementary methods: Volume of Fluids (VOF) and Lattice Boltzmann Method (LBM) to capture the complex dynamics. The VOF and LBM simulation retrieved the jerky dynamics of the wicking in the desired direction. The jerky dynamics brings into play inertia during burst phases. The fluid inertia has an influence on the onset of depinning and has therefore to be taken into account in the design of the channel shapes. This fact is reminiscent of the imbibition in porous media for which inertial forces may affect fluid front displacement dynamics as proved by [50]. In the negative direction, i.e., unwanted direction, the simulation reveals how water may bypass the edges responsible to the pining effect.

Numerical simulation provides a limit angle around 70° beyond which the liquid will no longer spread through the micro-channels, whereas the optimal contact angle of the asymmetric spreading is about 30°. This value is significantly below metal/water values which are typically superior to 50°. Consequently, the parameter range for designing efficient microchannels is narrow, pointing out the key role of numerical simulation in the design process.

LBM simulation tracks the gas/liquid interface better than VOF. Nevertheless, the BGK code in LBM that we used can be unstable at realistic Reynolds numbers. Other classes of LBM deserve to be tried in the future, as proposed in [51,52] for instance. Thus, a first outlook is related to the improvement of the numerical tools for the numerical simulation. Moreover, taking into account the surface defects in the simulation is a crucial outlook since our study showed that defects impact the directional wicking. An economic way to perform this task is to replace topographic defects by wettability defects. It is known to mimic the effect of surface irregularities [40].

Even if the fin shape proves to be efficient for the directional wicking, we aim at applying the numerical simulation to perform a parametric study and an optimization of the geometry structure.

In the PEM fuel cell framework, a path for improving the efficiency of the directional wicking is the surface treatment, e.g. [53], in order to achieve a lower contact angle between steel and water closer to the optimal one. In addition, in the parallel flow field, the gas shear flow drags the liquid. It is known that in confined geometry the liquid may block the gas flow [54]. However, the interaction of shear flow and wicking has thus far not been studied. Such a study will be of crucial interest in applying the structure in the parallel flow field of PEM Fuel Cells.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/1996-194
4/14/3/490/s1, Video S1: Videos of experiments and simulations.

**Author Contributions:** Conceptualization, N.A. and V.P.S.; Methodology, N.A.; Data curation, N.A.; Formal Analysis, N.A. and P.B.; Investigation, N.A.; Software, N.A. and V.P.S.; Validation, N.A. and P.B. and M.-C.N.; Visualization, N.A.; Writing—original draft preparation, N.A.; Supervision, P.B. and V.P.S. and M.-C.N.; Writing—review-editing, N.A. and P.B. and M.-C.N.; Project administration, V.P.S. and P.B. and M.-C.N.; Funding acquisition V.P.S. and P.B. and M.-C.N.; Resources V.P.S. and P.B. and M.-C.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Federal Ministry for Economic Affairs and Energy in Germany as a part of ZIM (Zentrales Innovationsprogramm Mittelstand) with the project number ZF4040702CL6.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding authors.
**Acknowledgments:** N.A. and P.B. would like to thank Stéphane Ruy, UMR1114 EMMAH, for his support throughout the project.

**Conflicts of Interest:** The authors declare no conflicts of interest.

**Abbreviations**

The following abbreviations are used in this manuscript:

SLM Selective Laser Melting

IPA Isopropyl alcohol

CAD Computer Aided Design

LBM Lattice Boltzmann Method

VOF Volume of Fluid

BGK Bhatnagar-Gross-Krook

PEMFC Proton-exchange membrane fuel cells

**References**

1. Mettu, S.; Chaudhury, M.K. Motion of Drops on a Surface Induced by Thermal Gradient and Vibration. *Langmuir* **2008**, *24*, 10833–10837, doi:10.1021/la801380s.

2. Ahmadlouydarab, M.; Feng, J.J. Motion and coalescence of sessile drops driven by substrate wetting gradient and external flow. *J. Fluid Mech.* **2014**, *746*, 214–235, doi:10.1017/jfm.2014.133.

3. Ballerini, D.R.; Li, X.; Shen, W. An inexpensive thread-based system for simple and rapid blood grouping. *Anal. Bioanal. Chem.* **2011**, *399*, 1869–1875, doi:10.1007/s00216-010-4588-5.

4. Martinez, A.W.; Phillips, S.T.; Nie, Z.; Cheng, C.M.; Carrilho, E.; Wiley, B.J.; Whitesides, G.M. Programmable diagnostic devices made from paper and tape. *Lab Chip* **2010**, *10*, 2499, doi:10.1039/c0lc00021c.

5. Parker, A.R.; Lawrence, C.R. Water capture by a desert beetle. *Nature* **2001**, *414*, 33–34, doi:10.1038/35102108.

6. Guadarrama-Cetina, J.; Mongruel, A.; Medici, M.G.; Baquero, E.; Parker, A.R.; Milimouk-Melnytchuk, I.; González-Viñas, W.; Beysens, D. Dew condensation on desert beetle skin. *Eur. Phys. J. E* **2014**, *37*, doi:10.1140/epje/i2014-14109-y.

7. Gürsoy, M.; Harris, M.; Carletto, A.; Yaprak, A.; Karaman, M.; Badyal, J. Bioinspired asymmetric-anisotropic (directional) fog harvesting based on the arid climate plant Eremopyrum orientale. *Colloids Surf. A Physicochem. Eng. Asp.* **2017**, *529*, 959–965, doi:10.1016/j.colsurfa.2017.06.065.

8. Aliabadi, M.; Zarkesh, A.; Mahdavinejad, M. Bioware fog collectors: The Texas horned lizard as a model for a biomimetic fog-harvesting. *Mater. Res. Express* **2018**, *5*, 115502, doi:10.1088/2053-1591/aadab4.

9. Muto, K.; Ito, S.; Ishii, D. Liquid Transport in Bio-Inspired Capillary-Driven Open-Air Channels. *MRS Adv.* **2017**, *2*, 1111–1116, doi:10.1557/adv.2017.77.

10. Brown, P.S.; Bhushan, B. Bioinspired materials for water supply and management: Water collection, water purification and separation of water from oil. *Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.* **2016**, *374*, 20160135, doi:10.1098/rsta.2016.0135.

11. Yao, X.; Song, Y.; Jiang, L. Applications of Bio-Inspired Special Wettable Surfaces. *Adv. Mater.* **2010**, *23*, 719–734, doi:10.1002/adma.201002689.

12. Li, J.; Guo, Z. Spontaneous directional transportations of water droplets on surfaces driven by gradient structures. *Nanoscale* **2018**, *10*, 13814–13831, doi:10.1039/c8nr04354j.

13. Liu, Q.; Xu, B. A unified mechanics model of wettability gradient-driven motion of water droplet on solid surfaces. *Extrem. Mech. Lett.* **2016**, *9*, 304–309, doi:10.1016/j.eml.2016.09.004.

14. Huang, J.J.; Huang, H.; Wang, X. Numerical study of drop motion on a surface with stepwise wettability gradient and contact angle hysteresis. *Phys. Fluids* **2014**, *26*, 062101, doi:10.1063/1.4880656.

15. Chowdhury, I.U.; Mahapatra, P.S.; Sen, A.K. Self-driven droplet transport: Effect of wettability gradient and confinement. *Phys. Fluids* **2019**, *31*, 042111, doi:10.1063/1.5088562.

16. Feng, J.; Rothstein, J.P. One-way wicking in open micro-channels controlled by channel topography. *J. Colloid Interface Sci.* **2013**, *404*, 169–178, doi:10.1016/j.jcis.2013.02.052.

17. Buchberger, G.; Hischen, F.; Comanns, P.; Baumgartner, R.; Kogler, A.; Buchsbaum, A.; Bauer, S.; Baumgartner, W. Bio-inspired Microfluidic Devices for Passive, Directional Liquid Transport: Model-based Adaption for Different Materials. *Procedia Eng.* **2015**, *120*, 106–111, doi:10.1016/j.proeng.2015.08.576.

18. Muto, K.; Ishii, D. Effects of anisotropic liquid spreading on liquid transport in arrow-like micropillar arrays. *Colloids Surf. A Physicochem. Eng. Asp.* **2018**, *544*, 86–90, doi:10.1016/j.colsurfa.2018.02.023.

19. Abbaspour, N. Numerical and Experimental Approach for Bipolar Plates in PEM Fuel Cells: Novel Designs of Fluid Domain. Ph.D. Thesis, Avignon Université, Avignon, France, 2020.

20. Badduri, S.R.; Srinivasulu, G.N.; Rao, S.S. Influence of bio-inspired flow channel designs on the performance of a PEM fuel cell. *Chin. J. Chem. Eng.* **2020**, *28*, 824–831, doi:10.1016/j.cjche.2019.07.010.
21. Iranzo, A.; Arredondo, C.; Kannan, A.; Rosa, F. Biomimetic flow fields for proton exchange membrane fuel cells: A review of design trends. *Energy* **2020**, *190*, 116435, doi:10.1016/j.energy.2019.116435.

22. Wang, J.; Wang, H. Flow-Field Designs of Bipolar Plates in PEM Fuel Cells: Theory and Applications. *Fuel Cells* **2012**, *12*, 989–1003, doi:10.1002/fuce.201200074.

23. Chen, T.; Tang, M.N.; Liu, S.H.; Zhang, J.W.; Yang, L. Flow Field Plate Structure for Proton-Exchange Membrane Fuel Cell. *CN 107946610 B*, 22 November 2017.

24. Mossman, A.; Wells, B.; Barton, R.; Voss, H. Passive Recovery of Liquid Water Produced by Fuel Cells. U.S. Patent 829,871,9B2, 30 October 2012.

25. Ahmadi, N.; Rezazadeh, S.; Dadvand, A.; Mirzaee, I. Study of the Effect of Gas Channels Geometry on the Performance of Polymer Electrolyte Membrane Fuel Cell. *Periodica Polytech. Chem. Eng.* **2017**, doi:10.3311/ppch.9369.

26. Salah, Y.B.; Tabe, Y.; Chikahisa, T. Gas Channel Optimisation for PEM Fuel Cell using the Lattice Boltzmann Method. *Energy Procedia* **2012**, *28*, 125–133, doi:10.1016/j.egypro.2012.08.046.

27. Moumen, N.; Subramanian, R.S.; McLaughlin, J.B. Experiments on the Motion of Drops on a Horizontal Solid Surface Due to a Wettability Gradient. *Langmuir* **2006**, *22*, 2682–2690, doi:10.1021/la053060x.

28. Chaudhury, M.K.; Chakrabarti, A.; Daniel, S. Generation of Motion of Drops with Interfacial Contact. *Langmuir* **2015**, *31*, 9266–9281, doi:10.1021/la504925u.

29. Chaudhury, M.K.; Whitesides, G.M. How to Make Water Run Uphill. *Science* **1992**, *256*, 1539–1541, doi:10.1126/science.256.5063.1539.

30. Modroukas, D.; Modi, V.; Fréchette, L.G. Micromachined silicon structures for free-convection PEM fuel cells. *J. Micromech. Microeng.* **2005**, *15*, S193–S201, doi:10.1088/0960-1317/15/9/s04.

31. Piri, H.; Bi, X.T.; Li, H.; Wang, H. 3D-printed fuel-cell bipolar plates for evaluating flow-field performance. *Clean Energy* **2020**, *4*, 142–157, doi:10.1093/ce/zkaa007.

32. Yang, G.; Yu, S.; Mo, J.; Kang, Z.; Dohrmann, Y.; List, F.A.; Green, J.B.; Babu, S.S.; Zhang, F.Y. Bipolar plate development with additive manufacturing and protective coating for durable and high-efficiency hydrogen production. *J. Power Sources* **2018**, *396*, 590–598, doi:10.1016/j.jpowsour.2018.06.078.

33. Jayakumar, A.; Singamneni, S.; Ramos, M.; Al-Jumaily, A.; Pethaiah, S. Manufacturing the Gas Diffusion Layer for PEM Fuel Cell Using a Novel 3D Printing Technique and Critical Assessment of the Challenges Encountered. *Materials* **2017**, *10*, 796, doi:10.3390/ma10070796.

34. Yang, D.; Krasowska, M.; Priest, C.; Popescu, M.N.; Ralston, J. Dynamics of Capillary-Driven Flow in Open Microchannels. *J. Phys. Chem. C* **2011**, *115*, 18761–18769. doi:10.1021/jp2065826.

35. Bliznyuk, O.; Jansen, H.P.; Kooij, E.S.; Zandvliet, H.J.W.; Poelsema, B. Smart Design of Stripe-Patterned Gradient Surfaces to Control Droplet Motion. *Langmuir* **2011**, *27*, 11238–11245, doi:10.1021/la201671w.

36. The Collected Works of J. Willard Gibbs. *Nature* **1929**, *124*, 119–120, doi:10.1038/124119a0.

37. Abbaspour, N.; Beltrame, P.; Néel, M.-C.; Schulz, V. Geometrically driven liquid wicking: Numerical study and experimental validation. *Eur. Phys. J. Appl. Phys.* **2020**, *91*, 31101, doi:10.1051/epjap/2020200193.

38. Prosperetti, A.; Tryggvason, G. *Computational Methods for Multiphase Flow*; OCLC: 437186667; Cambridge University Press: Leiden, The Netherlands, 2007.

39. Comanns, P.; Buchberger, G.; Buchsbaum, A.; Baumgartner, R.; Kogler, A.; Bauer, S.; Baumgartner, W. Directional, passive liquid transport: The Texas horned lizard as a model for a biomimetic 'liquid diode'. *J. R. Soc. Interface* **2015**, *12*, 20150415, doi:10.1098/rsif.2015.0415.

40. Beltrame, P.; Knobloch, E.; Hänggi, P.; Thiele, U. Rayleigh and depinning instabilities of forced liquid ridges on heterogeneous substrates. *Phys. Rev. E* **2011**, *83*, doi:10.1103/physreve.83.016305.

41. Hirt, C.; Nichols, B. Volume of fluid (VOF) method for the dynamics of free boundaries. *J. Comput. Phys.* **1981**, *39*, 201–225.

42. Brackbill, J.; Kothe, D.; Zemach, C. A continuum method for modeling surface tension. *J. Comput. Phys.* **1992**, *100*, 335–354, doi:10.1016/0021-9991(92)90240-y.

43. Malgarinos, I.; Nikolopoulos, N.; Marengo, M.; Antonini, C.; Gavaises, M. VOF simulations of the contact angle dynamics during the drop spreading: Standard models and a new wetting force model. *Adv. Colloid Interface Sci.* **2014**, *212*, 1–20, doi:10.1016/j.cis.2014.07.004.

44. Bonn, D.; Eggers, J.; Indekeu, J.; Meunier, J.; Rolley, E. Wetting and spreading. *Rev. Mod. Phys.* **2009**, *81*, 739–805, doi:10.1103/revmodphys.81.739.

45. Cox, R.G. The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow. *J. Fluid Mech.* **1986**, *168*, 169, doi:10.1017/s0022112086000332.

46. Schulz, V.P.; Abbaspour, N.; Baumeister, T.; Röder, T. Lattice-Boltzmann Simulation and Experimental Validation of a Microfluidic T-Junction for Slug Flow Generation. *ChemEngineering* **2019**, *3*, 48, doi:10.3390/chemengineering3020048.

47. Blow, M.L.; Yeomans, J.M. Anisotropic imbibition on surfaces patterned with polygonal posts. *Philos. Trans. R. Soc. A Math. Phys. Eng. Sci.* **2011**, *369*, 2519–2527, doi:10.1098/rsta.2011.0030.

48. Nguyen, V.; Park, W.G. A volume-of-fluid (VOF) interface-sharpening method for two-phase incompressible flows. *Comput. Fluids* **2017**, *152*, doi:10.1016/j.compfluid.2017.04.018.

49. Connington, K.; Lee, T. A review of spurious currents in the Lattice Boltzmann Method for multiphase flows. *J. Mech. Sci. Technol.* **2012**, *26*, 3857–3863, doi:10.1007/s12206-012-1011-5.
50. Moebius, F.; Or, D. Inertial forces affect fluid front displacement dynamics in a pore-throat network model. *Phys. Rev. E* **2014**, *90*, 023019, doi:10.1103/PhysRevE.90.023019.

51. Lallemand, P.; Luo, L.S.; Peng, Y. A lattice Boltzmann front-tracking method for interface dynamics with surface tension in two dimensions. *J. Comput. Phys.* **2007**, *226*, 1367–1484, doi:10.1016/j.jcp.2007.05.021.

52. Fakharil, A.; Bolster, D.; Luo, L.S. A weighted multi-relaxation time Lattice Boltzmann Method for multiphase flows and its Applications to partial coalescence cascades. *J. Comput. Phys.* **2017**, *341*, 22–43, doi:10.1016/j.jcp.2017.03.062.

53. Cho, J.S.; Kim, K.H.; Han, S.; Beag, Y.W.; Koh, S.K. Hydrophilic surface formation on polymers by ion-assisted reaction. *Prog. Org. Coat.* **2003**, *48*, 251–258, doi:10.1016/S0300-9440(03)00098-5.

54. Dietze, G.F.; Ruyer-Quil, C. Films in narrow tubes. *J. Fluid Mech.* **2015**, *762*, 68–109, doi:10.1017/jfm.2014.648.