Nanoflow over a fractal surface

This paper investigates the effects of surface roughness on nanoflows using molecular dynamics simulations. A fractal model is employed to model wall roughness, and simulations are performed for liquid argon confined by two solid walls. It is shown that the surface roughness reduces the velocity in the proximity of the walls with the reduction being accentuated when increasing the roughness depth and wettability of the solid wall. It also makes the flow three-dimensional and anisotropic. In flows over idealized smooth surfaces, the liquid forms parallel, well-spaced layers, with a significant gap between the first layer and the solid wall. Rough walls distort the orderly distribution of fluid layers resulting in an incoherent formation of irregularly shaped fluid structures around and within the wall cavities. C 2016 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). [http://dx.doi.org/10.1063/1.4958975]


INTRODUCTION
The applicability of micro-and nanoflows to a spectrum of disciplines has kindled academic interest on the effect of the material properties and surface roughness on the flow field. Previous Molecular Dynamics (MD) studies on atomically smooth surfaces have shown that increasing the mass, 1 bonding stiffness, 2 and wetting properties 3,4 of the walls reduces the boundary velocity. Other molecular models studied the effect of surface roughness on the fluid slippage by inserting protrusions periodically on an initially smooth plane. [5][6][7][8][9][10] The unanimous conclusion is that surface irregularities reduce the boundary velocity. Subsequent studies have further shown that increasing the depth or frequency of these protrusions further decreases the velocity in the proximity of the walls. 7,9,10 A more complex molecular model 11 used fractal theory to emulate the stochastic nature of surface roughness, inherently found on the faces of materials. 12 The authors varied the complexity of the geometry by adjusting the fractal dimension, a property quantifying the ability of a fractal to fill up space. The paper concluded that increasing the complexity of the surface irregularities lowers the boundary velocity, even if the average depth of the roughness remains constant. The fractal-based model, however, was limited to two dimensions, thus motivating the investigation of three-dimensional roughness effects on the fluid slippage. In brief, investigations of nanoflows over rough surfaces are scarce. Furthermore, to the best of the authors' knowledge, an investigation of three-dimensional nanoflow over a rough surface has not been performed yet; thus, this is the motivation for the present study.
We have employed Molecular Dynamics (MD) to model a three-dimensional nanoflow of liquid argon confined by two solid walls. Fractal theory is employed to design a surface roughness that is anisotropic in the coplanar directions. Surfaces of different roughness depth and wetting properties, i.e., the strength of the solid-liquid interactions, are considered.
We show that the roughness with an average depth of approximately 0.2 nm reduces the boundary velocity by approximately 88%. We attribute this decrease to the more diffusive environment close to the walls, resulting from atoms deflecting on the surface irregularities and colliding with the streamwise moving particles. Although deepening the roughness accentuates the velocity decrease, the reduction is small compared to the initial decline. Increasing the wetting properties of the solid also decreases the interface velocity, regardless of the surface roughness; however, for rougher channels, the effect is more obvious. Furthermore, we have compared the three-dimensional model with its two-dimensional equivalent. This comparison reveals that the slip-length decays more slowly as a function of the depth of roughness in the two-dimensional model. In turn, this shows that the roughness in the spanwise direction is significant and that two-dimensional models cannot accurately capture the physics of flows in rough nanochannels.
The information provided in this paper can be used as the basis for reduced models for the boundary velocity in micro-and nanofluidics. Furthermore, the modeling aspects employed for the design of the walls can be used in the broader field of tribology, where numerical models have shown that the surface geometry plays an important role, 13 as well as on the design and/or optimization of fluidic networks. [14][15][16]

Simulation method
Our model consists of liquid argon confined by two solid walls. The roughness on the inner solid surfaces differs between cases ( Figure 1). The dimensions of the simulation box in the x, y, and z directions are L x = 18σ, L y = 43.5σ, and L z = 18σ, respectively. The walls are normal to the y direction, and the periodic boundary conditions are used in the coplanar directions, i.e., x and z.
Before introducing any roughness, each wall is a collection of 11 (1 1 1) Face Center Cubic (FCC) lattice planes with density ρ = 4 ρσ 3 . We then introduce roughness using the multivariate Weierstrass-Mandelbrot (W-M) fractal function, 17,18 Equation ( that n max → ∞, we can select an appropriate, finite integer value using the relation, where L max is the sample size and L min is the period of the smallest cosine. The equation ensures that the wavelengths of the cosines will fully span the range from L min to L max . The second sum over m in Equation (1) superposes M ridges to create the final surface. The constant γ (γ > 1) controls the frequency density, i.e., smaller value of γ results in a smaller interval between frequencies. Φ m, n is a matrix of random phases for each value of n and m. The parameter D s is the fractal dimension and quantifies the ability of the fractal to fill up space; in the three-dimensional case, its value ranges between 2 < D s < 3. In this study, we use D s = 2.5 which models a highly irregular geometry. The parameter C determines the maximum amplitude of the cosines and is given by where G is the roughness constant and is the parameter we vary to control the depth of the protrusions.
Starting with perfectly smooth walls, we calculated the W-M function with its origin at the center of each surface. Solid atoms on the inner side of the calculated curve, i.e., closest to the liquid, were deleted.
After preparing the walls, we randomly placed liquid atoms between them. Due to the complex nature of the geometry, we used dynamic Voronoi tessellation on the position of the atoms to calculate the volume of the channel. 19 The number of liquid atoms varied between cases to keep a constant density of ρ = 0.84 mσ −3 , where m is the mass of a liquid particle.
The particle interactions were modelled using the 12-6 Lennard-Jones (LJ) potential, where ε is the depth of the potential well and quantifies the strength of the interaction; σ is the van der Waals radius; and r i j is the distance between particles i and j. The LJ parameters for the liquid interactions are ε l = 1ε and σ l = 1σ. The wetting properties of the channel refer to its ability to stay in contact with the liquid and are properties of interest to this study. We adjust the solid's wettability by varying the value of the solid-liquid interaction, ε sl , between 0.2ε ≤ ε sl ≤ 0.6ε. The molecular diameter between the solid and the liquid is fixed at σ sl = 0.75σ. The solid atoms are fixed onto their initial lattice sites through the spring potential, where k is the spring stiffness. For all the cases considered here, k = 500 εσ −2 .
To control the temperature of the system, the velocities of the solid were rescaled every time step. We did not tamper with the liquid atoms as this can result in unphysical behaviour. 20 Following an initial equilibration phase, the temperature of the system was set to T = 0.72 εk −1 b . To develop a Poiseuille flow, we applied a force equal to 0.02 εσ −1 in the x direction on all fluid molecules. 21 The time step of the simulation was τ ≈ 2 fs. For the time evolution of the system, the Newton's equations of motion were integrated using the Verlet algorithm. For the simulations we used the LAMMPS molecular dynamics simulator (http://lammps.sandia.gov/). 22 After an initial equilibration phase, we calculate the slip-length by extrapolating the average velocity profile past the beginning of the channel wall until it vanishes. Due to the irregular nature of the channel walls, we calculate the slip-length as the perpendicular distance between the vanishing point of the velocity profile and the centre-line of the rough wall geometry.

RESULTS AND DISCUSSION
Introducing roughness onto an initially smooth surface significantly reduces the velocity of the liquid particles flowing next to it ( Figure 2). Specifically, our simulations reveal that, regardless of the wetting properties of the walls, the slip-length, the distance beyond the wall where the extrapolated velocity profile vanishes, decreases exponentially with increasing average roughness depth (Figure 3(a)). Although roughness models with periodically inserted protrusions found similar trends, we observe a faster decay, with 88% of the decrease occurring when transitioning from a smooth wall to a roughness of approximately 0.2 nm average depth (G = 0.005σ). Following the initial decrease, the slip-length continues to drop more gradually.
In agreement with previous studies, 3 the boundary velocity also decreases as the strength of the solid-liquid interaction increases (Figure 3(a)). To correlate the effect of the solid's wetting properties on the effect of the surface roughness on the boundary velocity, we normalise the slip-length by its value in the smooth channel, L 0 (Figure 3(b)). As we increase the average depth of roughness, stronger solid-liquid interactions result in more liquid atoms entering the cavities leading to a greater reduction of the slip-length. However, as slippage is an interfacial effect, liquid atoms trapped in the deeper layers of the walls contribute less towards the reduction of the slippage compared to liquid atoms in the top layers. Therefore, as the roughness increases beyond the first few solid layers, the strength of the solid-liquid interactions becomes less important.
To shed light onto the underlying cause of this exponential dependence, we consider the behaviour of the liquid particles close to the channel walls. Next to smooth walls, the liquid forms well-defined layers, the first of which is located approximately 0.15 nm away from the solid (Figures 4(a) and 4(b)). This gap suggests minimal contact between the liquid and wall atoms, allowing the liquid to glide over the solid with little resistance (as illustrated in Figure 2).
On the contrary, in rough channels, the liquid intrudes the solid walls, as demonstrated by the additional, small density peak on the left of Figures 4(a) and 4(b). Particles deflect on surface protrusions and collide with streamwise moving particles causing the velocity decrease. The increased number of collisions also breaks up the density layers found next to the smooth walls (this is illustrated in Figures 4(a) and 4(b) by the lack of definition of the density oscillations in rough channels). Furthermore, as the depth of the roughness increases, the number of deflected particles also increases which further decreases the slip.
The spatial averaging of one-dimensional density profiles filters out potentially important information in anisotropic systems. Using three-dimensional iso-surfaces, we observe that in the presence of surface roughness, the liquid forms a discontinuous, intertwined collection of irregularly shaped structures, oriented at different angles from the wall ( Figure 5). The iso-surfaces also show that the liquid density inside the wall cavities is at least five times greater than what the one-dimensional profiles show: the one-dimensional profiles show a maximum density of ρ = 0.25 ρσ −3 whereas the iso-surfaces indicate densities of up to at least ρ = 1.2 ρσ −3 .
A minor digression: the iso-surfaces of Figure 5 correspond to a single value of the density. Considering the full spectrum of liquid densities, the liquid structure is even more chaotic. However, the plots become visually unappealing and complicated to interpret which is why we did not consider a wider range of values.
The findings of the present study compare well with the results of Ref. 23 showing an inverse correlation between the slip-length and the percentage of solid surface that is not parallel to the flow FIG. 4. (a) One-dimensional density profiles of the fluid across the width of the channel for a solid with low wettability (ε sl = 0.2ε), for different roughness depths. As irregularities are introduced on the solid walls, the density profiles lose in sharpness, suggesting the mixing of the liquid layers. A small peak also appears beyond the beginning of the wall indicating the intrusion of liquid in its cavities and the direct solid-liquid contact. (b) One-dimensional density profiles of the fluid across the width of the channel for a solid with high wettability (ε sl = 0.6ε). The local maxima have higher density values than the equivalent cases for the wall with the lower wetting. However, proportionally, the roughness seems to have the same effect on the profiles.
(oblique surface). We compared our MD results with the model of 23 where L s,eff is the effective slip-length; and the angular brackets denote the average of the slip-lengths calculated locally at every point on the surface, using the average velocity profile ( Figure 6). For lower percentages of oblique surface, the results of our model match with those calculated using the model very well. For larger percentages, the two curves deviate slightly, with our results decreasing slightly slower. We believe this is because the highly irregular geometry of our model makes certain oblique surfaces unavailable to the liquid particles.
Finally, we have investigated whether two-dimensional roughness models adequately describe nanoflows over rough surfaces. We use the two-dimensional W-M function, 24 to prepare models with geometrical characteristics corresponding to our three-dimensional cases. Although most parameters remain the same, the fractal dimension is reduced to its two dimensional counterpart D 2D = D 3D − 1 = 1.5. In turn, the value of G was adjusted to keep the depth scaling factor C the same (see Equation (3)).
Similarly to its three-dimensional equivalent, the slip-length in the two-dimensional model decays exponentially with increasing roughness depth ( Figure 7). However, the decrease is more gradual and resembles the behavior of the slip in geometries with periodically inserted protrusions. This indicates that momentum disturbances in the streamwise direction propagate in the spanwise direction through viscous forces and collisions between derailed particles. We therefore conclude that nanoflows over rough surfaces have a three-dimensional character that cannot be ignored. Although two-dimensional models can produce qualitatively correct conclusions, they significantly overestimate the velocity slip.  7. Slip length as a function of the roughness depth for the two-and three-dimensional models with the same roughness qualities. Although both curves decay exponentially, the slip-length decreases much more gradually in two-dimensional models.

CONCLUSIONS
The effects of realistic, three-dimensional, fractal-based surface roughness on the boundary velocity of a nanoflow were investigated. Our results indicate that the slip-length decreases exponentially with increasing roughness depth. Furthermore, the rate of decrease is greater than the predictions of roughness models with periodically inserted protrusions. As the wettability of the channel increases, the effect of the surface roughness on the slippage also increases. We attribute this behavior to liquid particles bouncing off surface irregularities and colliding on spanwise-moving particles, obstructing their intended trajectory. These collisions destructively affect the solid-like liquid planes found next to smooth surfaces. Instead, next to rough surfaces, the liquid particles form irregularly shaped anisotropic structures. Finally, by comparing our three-dimensional model to its two-dimensional equivalent, we show that nanoflows over rough surfaces have a three-dimensional character that two-dimensional models cannot capture accurately.