Slip theory-based numerical investigation of the fluid transport behavior on a surface with a biomimetic microstructure

This paper presented a numerical simulation of drag reduction on a superhydrophobic surface with a groove structure. The computational fluid dynamics (CFD) method was used to analyze the effect of the groove microstructure when droplets impacted a superhydrophobic surface. The simulation results revealed three main characteristics: (1) The distance the droplet spread was larger along the direction parallel to the groove and smaller perpendicular to the groove; (2) Two protruding small spheres were formed at the edge of the droplet along the groove direction; (3) During retraction, the droplets presented a narrow, cross-shaped morphology. The effect of the groove structure in the two-dimensional microchannel on drag reduction near the wall was analyzed based on slip theory by the coupled level set/volume of fluid (CLSVOF) method. The air in the superhydrophobic pit formed a low-velocity vortex, which made the fluid roll on the air surface. The rolling on the surface produced a velocity slip at the gas–liquid interface. In addition, the superhydrophobic surface had an obvious drag reduction effect in the laminar flow state, but the drag reduction effect in the turbulent state was not ideal and even increased the flow resistance at the wall.

In 1997, Barthlott et al. first observed the microstructure of a lotus leaf surface. They believed that the self-cleaning behavior of the lotus leaf was due to the microstructure on the surface and the epicuticular wax layer; this self-cleaning behavior is called the 'lotuseffect' (Barthlott & Neinhuis, 1997;Guo & Liu, 2007;Nakajima et al., 2000). Since this seminal work, surface superhydrophobicity has been investigated extensively. Reiner et al. investigated the wettability and selfcleaning properties of three kinds of superhydrophobic surfaces: silicon wafer specimens with different regular arrays of spikes hydrophobized by chemical treatment, replicates of water-repellent leaves of plants and commercially available metal foils that were additionally hydrophobized by means of a fluorinated agent. They found that water drops with kinetic impact energy could clean the surface perfectly (Fürstner, Barthlott, Neinhuis, CONTACT Chunbao Liu liuchunbao@jlu.edu.cn & Walzel, 2005). Their results showed that the coupled factors on the surface of these organisms exhibited a good boundary layer control mode in contact with droplets or fluids; thus, these surfaces had both drag reduction and superhydrophobic properties (Bixler & Bhushan, 2013). Hence, superhydrophobic surfaces have important applications in flow drag reduction due to their excellent hydrophobicity (Aljallis et al., 2013;Hyungmin & Sun, 2014;Lu, Yao, Hao, & Fu, 2010;Rothstein, 2010). The mechanism governing laminar flow drag reduction on a superhydrophobic surface was mainly the velocity slip on the interface (Choi & Kim, 2006;Ou & Rothstein, 2005). The slip effect is a prominent property of superhydrophobic surfaces. The slip effect reduces the shear force between the viscous fluid and the wall and makes the flow near the wall closer to the ideal flow. Cottin et al. found that a superhydrophobic surface with a microstructure similar to the surface of a lotus leaf could significantly reduce flow resistance. Their experimental and simulation results showed that a superhydrophobic surface with this particular microstructure could exhibit an obvious slip phenomenon (Cottin-Bizonne, Barrat, Bocquet, & Charlaix, 2003). Choi et al. found that the surface of 'nano turf' with superhydrophobic properties had a slip length of tens of microns (Choi, Westin, & Breuer, 2003). Li et al. found that when the microstructure of a superhydrophobic surface reached the micron level, the slip length was close to 50 μm (Li, Di, Li, Qian, & Fang, 2007).
Computation fluid dynamics (CFD) is a new interdisciplinary subject that integrates fluid mechanics and computer science. Based on the computational method, the approximate solution of the fluid control equation can be obtained by using the fast computational capability of a computer. Mou et al. applied CFD techniques to ascertain wind pressure distributions on and around various squared-shaped tall buildings (Mou, He, Zhao, & Chau, 2017). Some researchers have used the CFD method for studies involving modeling hydrogen production and numerically simulating the use of natural gas in a dual-fueled diesel engine (Akbarian et al., 2018;Faizollahzadeh Ardabili et al., 2018). Ramezanizadeh, M. et al. used the CFD method to analyze the performance of thermosiphons in heat exchangers based on experimental data (Ramezanizadeh, Nazari, Ahmadi, & Chau, 2019). In this paper, the CFD method was used to simulate the motion of droplets on the surface of bionic grooves.
Despite the consensus, due to the several factors, such as the complexity of turbulent flow, the uncertainty of the experiments, and the simplification of numerical simulations, the turbulence mechanism of drag reduction on a superhydrophobic surface continues to cause confusion. Studies have produced some contradictions, and the intrinsic mechanism of drag reduction has not been fully elaborated. A superhydrophobic surface with a groove structure was designed in this paper by choosing the slot-edge microstructure on the leaf surface of Setaria viridis as a bionic prototype. A computational model was constructed considering the tiny groove structure, and a relatively fine structured grid was generated. The motion morphology of the superhydrophobic surface fluid was numerically simulated and analyzed by volume of fluid (VOF) and interface tracking. First, the motion morphology of the droplets impinging on the superhydrophobic surface with a microstructure was calculated. The microchannel with a superhydrophobic structure was then established, and the flow drag reduction effect of the biomimetic structure under laminar and turbulent conditions was analyzed to reveal the drag reduction mechanism on the bionic superhydrophobic surface. This research is also a cross-disciplinary study that integrates the knowledge of bionic surfaces, interface science, fluid mechanics, cross-scale computing technology and other disciplines. Therefore, this study involves many disciplines and principles and has important scientific significance.

Wall slip theory
Navier (1827) first proposed the slip boundary condition and defined the slip velocity formula by introducing 'slip length': where L s is the slip length and ∂u x ∂y y=0 is the surface shear strain rate. As a result of the continuous understanding of slip flow, two typical forms of slip have been summarized: molecular and apparent, as shown in Figure 1. Molecular slip refers to slip that occurs at the liquid-liquid interface called the 'stagnant layer' where no flow exists near the wall. In contrast, apparent slip considers that there is a 'singular layer' near the wall, from which the velocity increases gradually. The maximum velocity of the 'singular layer' boundary is the slip velocity.  Superhydrophobicity is caused by solid surface roughness and extremely low surface free energy (Kannan & Sivakumar, 2008). Low surface free energy reduces the viscous force between the liquid and solid molecules near the wall, while the surface roughness structure leads to discontinuous three-phase contact lines and gas-liquid interface in fluid-solid contact, as shown in Figure 2.
Because of the viscous effect of the fluid, the fluid has a larger shear force near the wall, which leads to larger frictional resistance at the wall. The slip velocity reduces the shear force and the corresponding frictional resistance. Therefore, the existence of the slip velocity contributes to drag reduction (Xu, Chen, Mulchandani, & Yan, 2010). On a superhydrophobic surface, a certain roughness and microvoid structure reduce the contact area between the fluid and solid surface and reduce the viscous force between the attachment molecules and liquid molecules.
In the actual superhydrophobic drag reduction test, it is obvious that there exists a gas film composed of bubbles at the solid-liquid interface. The gas film converts a solid-liquid interface into a gas-liquid interface, which reduces the contact area between the liquid and solid wall, thereby reducing the frictional resistance of liquid near the wall. Therefore, in this paper, the wall was set as a rough surface with bulges and dimples, and the interior of the dimples were defined as being filled with air; the slip effect of air was simulated by air-liquid two-phase flow.

Interface capturing method
This paper coupled the VOF two-phase flow model and the level set function to determine the position of the interface by solving the change in the volume fraction of each phase in the mixed fluid to track the gas-liquid interface.
In the VOF method, a fluid volume function is defined to determine the ratio of the volume of the target fluid to the grid volume to capture the gas-liquid interface via numerical solution. In the computational domain where k-phase fluid exists, the function α k is defined as follows: If α k = 1, there is only k-phase fluid in the element grid; if 0 < α k < 1, there is multiphase fluid coexistence in the element grid; and if α k = 0, there is no k-phase fluid in the element grid. In the element grid, the volume fraction of each phase fluid should satisfy the following mathematical relation: Note that α k can be solved by the continuity equation.
The convective transport equation of the function is expressed as follows: where ν k is the velocity vector (m/s) and t is time (s). The level set method was proposed by Osher and Sethian in 1988. The basic idea of the level set method is to define a higher-order function Ø k to capture the motion of the gas-liquid interface, and the zero isosurface of the function Ø k is a dynamic interface at any time. The function is defined as follows: Note that Ø k > 0 indicates the external region of the kphase interface, Ø k = 0 indicates the interface of the twophase fluid, and Ø k < 0 indicates the internal region of the k-phase interface. The convection transport equation of the function is as follows: The coupled level set/VOF (CLSVOF) method overcomes the inaccuracy of the VOF method in calculating parameters, such as the normal direction and curvature, and ensures the mass conservation of the level set method in the process of convection transport and reinitialization of the function, which makes the capture of the phase interface more accurate and sharper (Son, 2003;Sussman & Puckett, 2000). The basic calculation procedure of the coupled calculation of the two methods is shown in Figure 3.

Surface tension model
The surface tension model is the continuum surface force (CSF) model proposed by Alamos. This model is solved by adding the source term to the momentum equation. In this paper, only gas-liquid two-phase fluids are included, so the continuous surface tension F CSF can be expressed as follows: where σ is the surface tension coefficient (N/m), κ is the surface curvature, n is the surface normal vector, α is the volume fraction (0 ≤ α ≤ 1), ρ is the fluid density (kg/m 3 ), and subscripts 1 and 2 represent the main and second phases, respectively. The internal and external pressure difference between both sides of the gas-liquid interface can be expressed by the product of surface curvature and surface tension coefficient: where p 1 and p 2 are the fluid pressures on both sides of the interface (Pa) and R 1 and R 2 are the radii of principal curvature (m). When droplets impact and spread on the surface, the wall adhesion is added to the surface normal phase as follows:n =n w cos θ +t w sin θ wheren w is the unit normal vector,t w is the wall tangential vector, and θ is the wall contact angle (°).

Basic control equations
Fluid flow is governed by three laws of conservation of physics (conservation of mass, conservation of momentum and conservation of energy), and the governing equations are the mathematical descriptions of these conservation laws. The numerical simulation process of droplet motion includes gas-liquid two-phase flow without considering heat exchange, and only the conservation of mass and momentum are solved numerically.
(1) Mass conservation equation For incompressible fluids, the density ρ is constant. According to the law of conservation of mass, the equation of conservation of mass can be obtained as follows: where v m is the average velocity vector (m/s). (2) Momentum conservation equation For incompressible fluid with constant viscosity, the momentum conservation equation is as follows: where ρ m is the average density (kg/m 3 ), P is the static pressure (Pa), μ m is the average dynamic viscosity (Pa·s), ρ m g is the gravitational volume force (N/m 3 ), and F CSF is the momentum source term caused by surface tension. Substance parameters ρ m and μ m can be obtained by the volume weighted averaging method, and the calculation formulas are as follows:

Validation of numerical-simulation-based liquid droplet behavior impacting a bionic grooved surface
The test results showed that the static contact angle of Setaria viridis could reach 150.8°, which indicated that the surface of the Setaria viridis was superhydrophobic. The surface of the Setaria viridis leaf was view under 100× magnification to obtain the SEM image, as shown in Figure 4(a). Figure 4(a)-III shows that there are parallel groove edges on the surface of the Setaria viridis leaf, and the groove edge direction is consistent with the direction of the main leaf vein. We simplified the surface microstructure of the original plant. Referring to the surface micromorphology and microstructure size of the bionic sample, a computational domain model was established (Bormashenko, Pogreb, Whyman, Bormashenko, & Erlich, 2007;Liu & Lange, 2006). The main region of the computational domain was a regular cuboid. To facilitate the meshing of the bionic microstructure region, the shape of the groove with an arc top structure in the bionic sample was simplified to a small cuboid with a regular distribution (Lee, Lim, Dong, Kwak, & Cho, 2013). Considering the computation time and simulation accuracy, the computational domain was divided into a highquality hexahedral grid. To accurately capture the motion pattern of and internal energy changes in the droplets near the wall, a local grid refinement method was used for the microstructure surface of the lower wall. In the simulation, 300,000, 850,000, 2,170,000 and 3,570,000 cells were used to test the grid independence. Based on the contact time between droplets and walls, it was found that the deviation in contact time between 3.57 million cells and 2.17 million cells was only 1.34%. Considering the calculation time and accuracy, a grid with 3,572,068 cells was selected. The computational domain physical model and grid model are shown in Figure 4(b,c), respectively. To accurately analyze the movement of the droplet on the grooved surface, a small droplet with a diameter of 2 mm was simulated using the adapt region function in Fluent software. The droplet impacted the bionic superhydrophobic surface at an initial velocity v = 1 m/s. Table 1 shows the boundary condition and parameter settings. Figure 5 shows the calculation results for the solution process.
For a transient solution, the time step is a very important parameter, and an improper setting may lead to a series of issues. Selecting an excessively large time step may cause divergence and low time resolution, whereas selecting an excessively small time step may cause the iteration count and the computation cost to increase. Therefore, it is very important to reasonably set the time step.
The time step Δt must be sufficiently small to analyze the time-related characteristics, which can be roughly estimated by the following formula: where Δx is the local grid size (mm) and v is the characteristic flow velocity (m/s). The minimum grid size is 0.01 mm, and the characteristic flow velocity is 1 m/s. According to the estimation formula, the time step Δt is 0.00001 s, which can ensure that all grid elements can be iteratively computed. Figure 5 shows the spreading morphology after the droplet impacts a superhydrophobic surface with a groove microstructure. The droplet spread more easily on the rough surface in the direction parallel to the groove than in the direction perpendicular to the groove; the spreading movement of the droplet was hindered in the latter, so the spread diameter of the droplet decreased. In the retraction stage, the droplet was retracted in a regular cross pattern, which was quite different from the smooth surface. The CFD calculation results were compared with the results in Ref. [Kannan & Sivakumar, 2008]. Two results are shown in Figure 6. Obviously, the whole process of droplet movement experienced four stages: spreading, contracting, rebounding and falling, and the error between the CFD results and the referenced test results was less than 5%.
For comparison, we also set the same contact angle on a smooth superhydrophobic surface; Figure 7 shows the results. There is a visible difference between the results of the grooved and smooth superhydrophobic surface. By comparison, three main characteristics of droplet motion on a grooved superhydrophobic surface can be summarized. The summary and analysis are as follows: (1) The distance the droplet spread along the direction parallel to the groove was larger than the distance perpendicular to the groove. The groove microstructure made the gas store up in the area between the grooves, which made the solid-liquid contact become gas-liquid contact when the droplet spread on the solid surface and reduced the contact area of the droplets with the solid wall. When spreading along the groove direction, the droplet was not subjected to any other resistance in the spreading direction except the surface tension. In the direction perpendicular to the groove, under the effect of surface tension and gravity, both sides of the groove structure hindered the flow of the water droplet, which hindered the spreading of the droplet. Thus, the overall spreading motion of the droplet formed a narrow, directional trajectory. (2) Two protruding small spheres were formed at the edge of the droplet along the groove direction. There were two reasons for this phenomenon. One reason is that the squeezing liquid accumulated at the edges of the droplet due to the downward movement of the central region of the droplet. In addition, the movement resistance was smaller along the direction of the groove; hence, the liquid in the central region preferably spread along the groove. Therefore, the liquid accumulation in the direction of the center groove was significantly greater than in the other directions.
(3) The retraction of droplets presented a narrow, crossshaped morphology. The liquid spreading perpendicular to the groove was hindered by the groove structure, and this liquid reached the maximum spreading state and entered the retraction process first. However, the liquid along the groove direction was still in a spreading state due to the small resistance. The anisotropy of the droplet motion and difference in velocity distribution caused a unique morphology of droplet motion on the surface of the groove structure.

Analysis of drag reduction on a superhydrophobic surface based on slip theory
Based on the wall slip theory, the drag reduction principle of a two-dimensional microchannel slip surface was analyzed by theoretical analysis and numerical simulation.
Our study was the first to consider fluid transport behavior on superhydrophobic surfaces based on slip theory. By establishing a bionic surface model of a superhydrophobic surface slip, the VOF two-phase flow model was used to simulate the flow state and drag reduction effect of the fluid flowing on a superhydrophobic surface. Figure 6. Comparison of the results between the CFD simulation and the test in Ref. [Kannan & Sivakumar, 2008].
The superhydrophobic surface of the microchannel model consisted of a number of microgrooves. To ensure the full development of the flow, a smooth wall was set with a certain distance between the superhydrophobic surface and inlet. Moreover, to avoid affecting the flow state of the superhydrophobic surface, a certain flow distance was set at the lower boundary and the outlet boundary. The composite boundary conditions of the solid-liquid contact and gas-liquid contact simplified the film form of the superhydrophobic interface to a certain extent and reduced the difficulty of the numerical simulation while ensuring the accuracy of the simulation.
A quadrilateral orthogonal grid was used for twodimensional microchannels with regular microprotrusions. To improve the calculation accuracy, a local grid refinement method was used near the gas-liquid interface. Prior to the numerical simulation, grid independence was verified. With the import and export pressure drop as the testing standard, four grid schemes of 5 × 10 4 , 1 × 10 5 , 1.5 × 10 5 and 2 × 10 5 were used for verification. The results showed that the calculated deviation in the grids with 2 × 10 5 and 1.5 × 10 5 cells was only 0.74%. Considering the calculation accuracy and time, the grid number was 201,735. Figure 8 shows the two-dimensional microchannel computational domain model and grid model. Table 2 shows the boundary condition and parameter settings.
The outlet pressure was monitored during the calculation. When the monitoring curve was stable and the residual error was lower than the convergence  standard, the result of the solution was considered to be convergent. Figure 9(a) shows the gas-liquid two-phase distribution cloud map defined in the initial state, in which the red part is the water phase, the blue part is the gas phase, and the two-phase interface is a straight line. Figure 9(b) shows the gas-liquid two-phase distribution cloud map of the surface of the microstructure after flow stabilization. The boundary line of the gas-liquid twophase changes from the initial straight line to a concave parabola, which is similar to the result in the superhydrophobic experiment (Ou, Perot, & Rothstein, 2004). After reaching steady state, the groove of the superhydrophobic surface can effectively repose the air, forming a gas-liquid contact surface with less shear force, thereby generating a velocity slip at the contact surface, effectively reducing the frictional resistance and viscous resistance during the fluid flow and achieving the effect of drag reduction.
The gas-liquid boundary line in the pit at position A in Figure 9(b) was extracted. A polynomial fit of the coordinates at the interface produced the result in Figure 9(c). Figure 9(c) shows that the gas-liquid interface had a distinct curvature and could be approximated as a parabola. This curvature was the result of the combined effect of the presence of air in the groove and the surface tension of the liquid. According to the extracted data, a mathematical expression of the boundary line distribution was obtained. The fitting function is shown in formula (15). According to the polynomial, the maximum depth of the air-liquid interface depression in the pit was 0.995 μm.
where R 2 is the fitting accuracy. The red curve in Figure 10(a) is the gas-liquid boundary. The upper part of the boundary is air, and the lower part is water. The boundary between the air and the water is obviously an arc. Due to the presence of air in the pit, the flow of the superhydrophobic surface can be divided into two categories: the nonslip flow on the solid-liquid interface and the slip flow on the gas-liquid interface at the pit. Figure 10(b) shows that the air residing in the superhydrophobic pit forms a low-speed vortex, which makes the fluid roll on the air surface when moving on the superhydrophobic surface. Figure 10(c) shows the velocity distribution in the longitudinal direction (i.e. the direction of the arrow in Figure 10(a) near the gas-liquid interface at the middle position in a single pit), and Figure 10(d) shows an enlarged view of portion A in Figure 10(c). The fluid velocity at the air interface is not zero, which indicates the existence of slip velocity. Farther away from the gas-liquid interface, the fluid velocity increased parabolically, which is consistent with the trend of theoretical velocity distribution. Under the action of a superhydrophobic surface structure, both the velocity slip and flow structure changed. The vortex that formed inside the groove opposite to the direction of the external flow field slowed the change in the flow state of the superhydrophobic surface, making the fluid flow near the wall more stable and the drag reduction effect more obvious.  The dimensionless pressure drop ratio is introduced to measure the drag reduction effect of the superhydrophobic surface in microchannels.
The dimensionless pressure drop ratio can be calculated as follows: where p N is the theoretical pressure drop of the nonslip surface (Pa) and p s is the actual pressure drop of the slip surface (Pa). Figure 11(a) shows the relationship between the inlet velocity of the microchannel and the dimensionless pressure drop ratio of the microchannel in the laminar flow. The dimensionless pressure drop ratio is always greater than zero, which indicates that the existence of Figure 11. Relationship between the inlet velocity and dimensionless pressure drop ratio: (a) Relationship between the inlet velocity and dimensionless pressure drop ratio of the microchannel in the laminar flow regime (108.5 < Re < 651); (b) Relationship between the inlet velocity and dimensionless pressure drop ratio of the microchannel in a turbulent flow regime (10850 < Re < 65100).
a superhydrophobic surface can effectively reduce the energy loss of the fluid at the wall in the laminar flow. However, with increasing inlet velocity, the dimensionless pressure drop ratio decreases, and the drag reduction effect of the superhydrophobic surface also decreases. Figure 11(b) shows the relationship between the inlet velocity of the microchannel and the dimensionless pressure drop ratio in the turbulent flow state. As the inlet velocity increases, the dimensionless pressure drop ratio fluctuates at approximately zero, which indicates that the drag reduction effect of the superhydrophobic surface in the turbulent state is not obvious. This phenomenon may be related to the selected microstructure. It is well known that bioinspired shark-skin riblets can provide a maximum drag reduction of nearly 10% (Bechert, 1997). However, they could achieve an effect with any size. In their study, the optimal riblet geometries were defined according to the applications (Hough, 1980;Walsh, 1982). Our previous studies also indicated this fact (Liu, Sheng, Yang, & Yuan, 2018). Hence, although the drag reduction effect of the microstructure selected in this study was not obvious in a turbulent flow, this does not mean that a superhydrophobic surface has no drag reduction effect in a turbulent flow. This study focused on the flow numerical simulation method.
Even as the inlet velocity increases, the dimensionless pressure drop ratio also has a negative value, which increases the energy loss on the solid surface.
In summary, the superhydrophobic surface had a relatively obvious drag reduction effect in the laminar flow state, and the drag reduction effect in the turbulent state was not ideal and even increased the flow resistance at the wall. Therefore, future research should focus on a superhydrophobic surface that can achieve drag reduction in a turbulent state to develop the application of superhydrophobic drag reduction.

Conclusion and outlook
This study used the CLSVOF method to track the gas-liquid interface. The movement of the droplet on the superhydrophobic surface with a groove microstructure and the drag reduction principle of the two-dimensional microchannel slip surface were analyzed.
Referring to the surface microstructure and size of the bionic prototype, the computational domain of a numerical simulation of droplets impinging on a superhydrophobic surface was established. Simulation analysis revealed that the droplets were more likely to spread in the direction of the groove than perpendicular to the groove after impacting the groove microstructures, and two protruding small spheres appeared on the edge of the droplets along the groove. In the reversion stage, the droplet motion morphology presented a narrow, crossshaped morphology due to the anisotropy of the droplet motion and the difference in velocity distribution. A twodimensional superhydrophobic microchannel model was established, and the flow state and drag reduction effect of the fluid flow through the superhydrophobic surface were simulated and analyzed by using a VOF two-phase flow model. The parabolic gas-liquid interface was formed when fluid flowed through the rough superhydrophobic surface, and the slip velocity was formed to reduce the friction resistance. In addition, the superhydrophobic surface had an obvious drag reduction effect in the laminar flow, while the drag reduction effect in the turbulent state was not ideal and even increased the flow resistance at the wall.
The drag reduction effect of the superhydrophobic surface under turbulent conditions was not satisfactory. Therefore, the research and development of superhydrophobic surfaces that achieve the effect of drag reduction in the turbulent state are important directions for the application of superhydrophobic drag reduction. Other biomimetic superhydrophobic structures are also a research direction.