Abstract
Microchannel flow shows a fascinating background on a lot of engineering problems. In order to shed a light on the effect of the surface morphology of microchannels on fluid flow, differently shaped and arranged artificial elements constitute channels with different morphology and numerical simulation based on lattice Boltzmann method is conducted. The impact of micro effect is also stressed by comparing the results considering and not considering it in the same channel model. Analysis of flow details shows the difference of the morphology effect on fluid flow, which differs by the shape and density of the elements’ array. The permeability of channels shows a specific relationship with the density of artificial elements, and differences are found between varied shapes and the existence of micro effects. Further research is carried based on more complex channels with arrays of fractal-character artificial elements. As elements in the channel can be divided into main summits and subsummits, their different roles of the effect on the fluid flow is investigated. The result shows that the permeability will not change if main summits are kept in channels while all subsummits are removed to make a distinct simplification of the morphology. This discovery is furtherly ensured numerically by a test on a channel created with the profile of a rough rock surface. The finding for morphology effect on fluid flow can supply a reference for the prediction of the permeability of complex channels or fractures.
1. Introduction
Flow in microchannels is a foundational problem in both fluid dynamic research and a wide range of engineering applications, including the traditional engineering works on hydraulic machines [1, 2] and the modern manufacturing of microfluidic devices [3–5]. As the prosperity of microelectromechanical system (MEMS), the study on the microchannel flow is motivated by the authentic flow pattern which is crucial for understanding and resolving some challenges in microfluidic techniques [6–8]. Meanwhile, with the ever-increasing demand for fine-grained flow characters in modern industry, such as the exploitation of unconventional resources, some special channels like nature fractures in the shale which may play the dominating role in the matter transport should be treated seriously [9–12].
Although the channel flow problem with different concerns has been studied from different aspects of simplification for more than two centuries [13, 14], remarkable results are still reported constantly with the development of theories and experimental techniques [10, 15–17]. The flow channel has a series of surface properties like roughness, corrugation, temperature, and other physical chemistry differences [9]. The surface morphology, whether it is small as the roughness or large as the bend, is one of the main features of the channel which has a great impact on the flow behavior. Theoretical researches based on fluid dynamics theories come as early as a nonmorphology effect considered Navier-Stokes solution to predict the flow rate of a single fracture channel named local cubic law (LCL) [18]. Considering the complex channel boundary effect, Berman studied laminar flow in channels with porous boundary and gives a quantitative description for flow properties [13]; Green and Naghdi [19] utilized the theory of directed fluid sheets to investigate the viscosity flow in boundary moveable channels. Based on a series of simplifications, these meaningful theoretical works depicted the main configuration for the flow with complex channel boundaries. Simplifications on the morphology of channels usually lead to the revealing of flow details impossible which is in otherwise the foundation of available analytical results. However, details of the flow field, especially impacted by morphology effect of channels, is very essential for modern engineering.
Based on simple models, it is possible to reveal basic physical patterns of the fluid flow, and nowadays, as the numerical simulation becoming popular and easy accessing, complex flows can be studied deeply with consideration of different channel characters [20, 21]. One of the powerful numerical methods for complex flows is the lattice Boltzmann method (LBM), which has a mesoscale-oriented background and fits for capturing flow details. The high computational efficiency is also one of the advantages it holds. Due to these talents, the LBM is widely applied to investigate the flow inner complex channels under different conditions [22–24]. Tan et al. [25] studied the rarefied gas flow with grooved walls using LBM; with interests on the drag reduction ability of the surface structures, they found that grooves with appropriate size and geometry can reduce shear forces effectively which is a good reference for airship designing. Inspired by their research, the importance of surface morphology effect on flow should also be emphasized in a constrained space.
For microchannel flows, the flow rate is one of the important results concerned by researchers, and it is convinced that the wall structures or wall roughness should be considered as an accredited factor for it [11, 26]. A number of models for predicting flow rates in the channels with complex boundary structures have been proposed [27–30]. According to the differences of the structure and roughness on the channel surface, approaches to model morphological properties of channels have been utilized in literatures, including the aperture model [31], the fractal model [32, 33], the joint roughness coefficient [34], and etc. However, the discussion is still ongoing and a comprehensive and detailed research on the flow through a surface-complex channel is still in great need [35–37]. To reveal the mechanism of the surface morphology effect, another way to simplify the geometrical property is to create artificial elements attaching to the channel, and this method is often used for its flexibility and controllability [38]. With kinds of controllable elements, a freely created parametric channel is useful to investigate the flow mechanism locally and globally.
An important fact worth noticed indicates that when it comes to the microchannel, which has a high surface-to-volume ratio, the microstructures on the channel surface will occupy a relatively larger proportion than of the macroscale. Due to this scale-changing effect, naturally or artificially generated microstructures like obstacles or roughness may play a more influential role in the flow behavior [39]. And this is also one of the factors why the microscale flow is distinctive from the macro.
Given all the importance above, to comprehend the mechanism of flow through microchannels with different morphological surface structures, numerical simulations to examine cases with a series of boundary structures are fulfilled in this paper. Lattice Boltzmann method (LBM) is chosen as the main simulation technique to deal with the Kn effect of the microchannel and flows passing complex geometries on the channel surface. Remarkable effects caused by the surface structures are well captured, and further analysis with more numerical experiments is disposed of.
2. Modeling and Method
2.1. Channels with Surface Structures
The evaluation of geometrical characters of channel surface always shows challenges. Whereas the fractal model is convinced as one of the best approaches which can be used to obtain analytical conclusions [40] and the statistical model shows close relationship to real channel surface [41, 42], alternatively we choose the artificial elements approach which is self-defined in details and shows advantages for the specific investigation from quantities and can work for the revealing of laws of local fluid flow. Although the channel model with artificial elements has been used for a long time with different concerns, we show an unlike way of utilizing it by adding or removing each element solely. In this case, more examples are calculated which increase the computation cost several times. However, it is worth because more data gives a detailed reveal for morphology effect.
A series of simple Euclidean geometries are added on one side of the wall to construct a complex channel. Rectangles, semicircles, and triangles are selected for each channel, respectively. To have a better understanding of how these surface structures exert their influence on the flow behavior, the minimum fluid-pass area is fixed realizing by setting every elements with same height (see Figure 1). So the differences on how each kind of geometry influences the fluid flow will be mainly affected by their geometric shapes. And the differently shaped elements are also set at the same longitudinal position of the field.

(a)

(b)

(c)
Considering the fluid flow through a void-rich medium is generally estimated by Darcy’s law (see equation (1)), which is practical in engineering application and valid for low Reynolds number flows (approximately 1~10). Darcy’s permeability in the law is chosen to evaluate the flow through the complex channel (see equation (2)). Darcy’s law indicates that the flow rate is proportional to the pressure difference and inversely proportional to the fluid viscosity, and the permeability shows intrinsic properties of the medium or the channel. where and are Darcy’s permeability and real permeability, respectively; is the dynamic viscosity of the fluid; is the kinematic viscosity of the fluid; is the mean density; defines the mass flux at the exit; is the length of the channel allied with flow direction, and is the aperture of fracture; defines the pressure difference between the entrance and the exit, is the area of the exit.
2.2. Lattice Boltzmann Method
Lattice Boltzmann method (LBM) is the development of the lattice gas automata (LGA) to provide an effective tool for computational fluid dynamics [43]. LBM is based on the special discretization of Boltzmann equation, which the fundamental equation in nonequilibrium statistical physics, and significantly different from the traditional computational fluid dynamics methods. For the mesoscale, nature of LBM provides it with another solution for fluid simulations and gives it edges on revealing flowing details [44, 45]. Meanwhile, the computational solution in LBM is quite simple for coding and extending.
The widespread D2Q9 model of LBM is adopted for the simulation. In this model, the single relaxation time Boltzmann-BGK equation replaces the original Boltzmann equation which has a complex collision term [46, 47] (see equation (3)). A 2nd-order truncation of the Maxwell velocity distribution function is used to solve the equilibrium distribution function [48] (see equation (4)). And the continuous velocity field is discretized into 9 directions on 2D (see equation (5) and Figure 2). where is the particle distribution function, and is the equilibrium distribution function; the subscript is the index for discrete velocity; , , and are the position, time, and time step, respectively; is the discretized velocity vector of the fluid element. where is the local flow density; is the velocity of fluid; is the sound speed in lattice unit and equals to (, δx, and δt are unit lattice distance and unit lattice time, respectively), and is the weighting coefficient for different : 1/3 (for ), 1/9 (for ), and 1/36 (for ).

Generally, the D2Q9 model fits for the isothermal systems where the compressibility and temperature changes of fluid can be neglected. Physically, the model can be considered as a simulation system linked with an infinite “thermal bath” that heat exchanges constantly so that the system keeps isothermal [49, 50]. Such a model is valid when the temperature difference caused by fluid flow or temperature influences on the viscosity and heat conductivity are insignificant [51, 52]. For the current problem of the flow in a channel, the D2Q9 is appropriate. The hydrodynamic quantities are related to the discrete distribution functions as below where is the pressure of the fluid.
In the following simulation, we set the mesh size for the channel model in lattice unit, the Reynolds number (Re) keeps approximately of 0.48 and Knudsen number (Kn) of 0.0111 (defined by ) for our simulation. Under this condition, surface structures cannot be ignored comparing to the channel aperture [53, 54].
According to flow regimes divided by Kn, this simulation has stepped into the transformation stage from nonslip continuous regime to continuous regime with boundary slipping [55, 56]. Existing literature reported the slip effect of the regular channels fully [57, 58]. But the complex channels with Kn effect has not been studied thoroughly. Therefore, the fluid flow behavior is simulated both considering and ignoring the Kn effect. This investigation is beneficial to supply a comparison for the flow in complex channels with and without the micro effect. A pressure-difference boundary condition is set to the entrance and the exit of the channel to activate the flow [59]. The bouncing back scheme is set as the boundary conditions for nonslip walls [60]. To consider the Kn effect, which includes the viscosity changes at the main flow and slippery at the boundary, a combination of the correction on the relaxation time and slip boundary condition is utilized.
According to Nie et al. [61], the corrected relaxation time replaces the original one,
Then, the viscosity can be expressed as
In order to capture the proper slippery velocity, Succi [62] implemented a slip boundary condition called DSB which has been widely used to account for the slippery phenomenon. This boundary condition can be expressed as the combination of the standard bounce-back boundary condition and the specular reflection boundary condition. For upper walls, it can be shown as where is the reflection coefficient, and it only depends on the fluid property and wall characters of the channel.
The simulation is conducted with parameters in dimensionless lattice units, and referential physical units according to dimensional analysis are shown in Table 1.
2.3. Model Validation
To make a validation for the model and method introduced above, numerical results obtained by a self-made code are compared with analytical solutions of the flow in a smooth channel. The theoretical analysis is provided by Karniadakis et al. [8], and this unified flow model for the microchannel flow is capable to estimate a wide range of flow regimes from the continuum regime to the transitional. As the model shows great agreement with the DSMC result and other credible simulation results, it is reliable to be used as a touchstone for our model and code. The nondimensional velocity profile of Karniadakis’s analysis [8] can be expressed as equation (12). where shows the profile of nondimensional velocity at the vertical section of flow direction and varies with and Kn; indicates the channel height. When is used in the simulation, the model has a second-order accuracy for a wide range of Kn.
As is shown in Figure 3, smooth-surface microchannels with Kn of 0.0111 and 0.0222 are taken as the benchmark. The LBM result shows good agreement with the analytical solution. The difference on the velocity profile can be observed when using the nonslip boundary condition and the slip boundary condition. For the case that Kn is at the magnitude of 10-2, the flow is typically at the slip flow regime, but the nonslip velocity profile deviates slightly from the slippery as Figure 3(a) shows. This result is inconsistent with some other reported researches [56, 57]. Therefore, it is reasonable that the slip effect is neglected for simplification in the engineering application. But Figure 3(b) implies that the velocity profile with and without the Kn effect show obvious differences, which only doubles the Kn number to 0.0222. Hence, the microchannel flow at the slip flow regime or not reached the transition flow regime should be investigated thoroughly considering both with and without Kn effect to give a credible reference.

(a)

(b)
3. Results and Discussions
Adopting the model and parameters above, present numerical results will be discussed in this section.
When fluid transports through a channel, flux is observed immediately and shows different values according to the morphological characters and physical properties of the channels. The flux of channels with different structures both under the slip and nonslip condition is shown in Figure 4. Firstly, the flux of the channel decreases with the increment of the proportion of area occupied by structures. The relative difference of flux between each specific case and the smooth channel is shown by the tagged numbers near the symbols. It is denoted that the difference caused by the slip and the nonslip effect is 4.9% of the smooth channel condition, while this quantity can reach about 20% with small elements or about 70% with large elements added in the channel. Besides, the rectangle elements show the strongest flux-reducing ability and the triangle elements show the weakest according to the comparison in Figures 4(a)–4(c). It is also shown that the difference of flux between slip and nonslip condition () declines with the , which means that the richness of elements makes the slip effect less remarkable. Although the speed of declining disperses with different altitudinal and shaped elements, this relationship can be observed globally. Therefore, it is convincing that boundary structure plays an important role on the flux-reducing effect which can reach up to several tenths percent on relative difference while the slip effect only has a lower than 10% influence on the flux. In this case, the morphology effect should be stressed and further investigated by the detailed information of the flow field.

(a)

(b)

(c)
3.1. The Velocity Profile and Permeability
The velocity profile at different sections was captured to make a detailed study of flow pattern. A group of element-rich channels is chosen for the analysis of flux-reducing effects. Flow profiles at the outlet position () of channels with a common relative height () of 0.25 are shown in Figure 5. According to Figure 5(a), it shows that the fluid flowing through a channel with several elements () can have a relatively low velocity when it flows out. As we all know, the mean velocity of the fluid flow through a relatively narrow channel is lower than a wider one under the same pressure-gradient condition. Therefore, the velocity of the element-rich case can be compared with a narrow one, as the section velocity profile in the channel with triangle elements is closer to a narrow channel case with in Figure 5(a). This result signifies the morphology effect of surface structures again. In Figure 5(b), three kinds of elements show diverse impact on the velocity profile. A slight deviation between the slip and nonslip condition is observed. Whether it is under slip condition or not, velocity profiles of three kinds of element-rich channels rank as the triangle one the highest and the rectangle one the lowest, which is consistent with the flux of three cases.

(a)

(b)
As trying to find a certain expression between flux and , real permeability of the channel is used to represent the fluid transfer property. The flux-reducing effect is then described by a dimensionless parameter . The permeability of a smooth channel stays constant as channels here have the same height (). Therefore, the flux reducing is actually reflected by the reciprocal of . Through this minor transfer of parameters, Figure 6 shows the relationship between and . As a linear fit in Figure 6 shows, there is a positive linear relationship between and . This relationship reveals that every single element has the same flow rate reducing effect in the channel. It also shows that two factors will exert influence on the permeability relationship. The first one should be concluded as morphology, as differently shaped elements have the different strength to change the flow rate, in accordance with the flux data (Figure 4) and section velocity pattern (Figure 5). The second factor is the slip effect. Under slip boundary condition, the slope of a relationship line will become more steep, which means that the surface structure’s flux reducing effect will be reinforced if the flow is microscale and steps into slip regime. Quantified flux-reducing effect is described by . With a larger value of as 0.14, the strongest flux-reducing effect is caused by rectangle elements under slip condition of , while the weakest effect is caused by triangle elements under the nonslip condition of .

In Figure 6, a linear fit for relationships between and is performed. Through analysis, the relationship of all conditions can be expressed in a common form as where is the slope of the fitting line and stays constant with the value 1. It can be shown that the relationship, (r, rectangle; c, semicircle; t, triangle) and , exists. This result indicates the importance of considering the slip effect. In the real situation, severe deviation would occur if microstructure and slip effect are not included in the permeability analysis [63].
3.2. Influence of Morphology
The result above just shows the overall properties of each channel, while they have all stressed the influence of the morphology effect. However, the mechanism of how the structure affects the velocity and channel flux is not clear. As the flux-reducing effect will be enhanced nonlinearly with the decreasing of channel height [64], the morphology effect on the channel flow will be more significant. Therefore, further analysis is needed to shed a light on the morphology effect.
Pressure plays a role of driving the fluid’s movement. In a standard Poiseuille flow condition, the pressure quantity declines evenly along the flow direction. However, if elements are added along the channel, the pressure will change accordingly (as Figure 7 shows). Comparing with , which is the pressure distribution of the smooth channel, pressure () is larger when the stream encounters the first two elements, while it is smaller when the stream passes each element. Especially, it is shown that the pressure increase at the upwind side of the element, and a slump occurs at the leeward side. Therefore, the pressure difference between the two sides leads to a force acting on the element. As the structure does not deform with the flow in this study, the element will push the fluid reversely, which leads to locally slowing down. Due to this obstruction, the pressure cannot deliver effectively, so, the pressure drops, the velocity decreases, and the flux reduces, which presents the decline of permeability [65]. Three types of elements show the same properties on the pressure profile, but the minor difference can be observed between each type. The pressure increases most at the upwind side of the first rectangle element both under the condition of the slip and nonslip and decreases most when the flow encounters rectangle elements, which indicates rectangle-shaped elements have the strongest ability to make the pressure to shake. The other two kinds of elements cannot compare with the rectangle of pressure perturbation. Hence, the flux differs in three kinds of channels.

With elements on the surface, not only the local flow path will be changed, fluid will also recirculate due to streamline separation [66]. This effect was observed by Brush and Thomson [67] in sinusoidal rough-walled fractures, while as there are obvious intervals between each element, it is not clear that how the separation will be affected by the morphology of elements. The contour and streamline present the velocity field of each channel with different morphology (see Figures 8–10). Streamlines nearing roughness elements show different characters of flow paths. Trigonal-shaped vortexes are found at the front and back side of rectangle elements, which indicates a strong separation of flow, while small vortexes locate at both the corners of semicircle elements. It also shows that obvious vortexes are absent near the triangle elements, which can be inferred as no separation occurs at this area. The rectangle element, whose upwind surface is completely perpendicular to the flow direction, will lead to the formation of the high-pressure zone at the upwind side, and it consequently causes strong vortex locally. For the semicircle element, its geometrical shape is relatively smooth for fluid to pass, so the vortex only occupies a small area at its corner. In the case of a triangle element, due to the smaller upwind angle of triangle geometry, the fluid only forms a very little distinct local vortex. The vortex signifies an independent movement of a portion of fluid breaking away from the mainstreaming and forming a closed circuit.

(a)

(b)

(c)

(a)

(b)

(c)

(a)

(b)

(c)
As the contour of velocity shows, regions with very low velocity are obvious both at the intervals of two elements and at the corners near the elements. To make a further analysis on this property, the velocity profile is extracted from the middle place of the intervals (as Figures 8(c), 9(c), and 10(c) shows). Firstly, the value of velocity at y-direction () signifies the tortuosity of the flow path as the channel aperture expands at the interval between elements. under the condition of slip boundary is larger than nonslip case, and the patterns are similar at intervals of three kinds of elements. has relatively small value and can be neglected comparing with the velocity at x-direction (); thus, will be mainly discussed. The profile of shows different trends in three different shaped cases. There is a long flat stage at the left side of the velocity profile of the rectangle element case, while the stage is not so flat in the semicircle element case. When the shape of elements is a triangle, the velocity increases straightly with at the left side of the profile. The velocity at the place where the height is constant with the elements’ peak is very low. Its value can be compared with the boundary layer at the opposite side. Thus, elements have a strong influence on the fluid flow at intervals. Figures 8–10 also show differences in the velocity value at the position, which indicates the velocity near the rectangle elements is the lowest and the velocity near the semicircle elements comes just after it. All the patterns mentioned shows the velocity at the interval is extremely lower than mainstream and differs with the shape of elements.
3.3. Further Evaluation of the Morphology Effect
Vortex which is independent of the mainstream and has a low speed is generated at the corners of surface structures of the channel. This phenomenon leads to the slowing down of the fluid at the intervals of elements. As always can be seen, the morphology of a fracture or a complex channel profile can be decomposed as main summits and subsummits, and structure-disassemble works can be done based on the wavelet decomposition [68]. In this case, to evaluate the effect of decomposed small structures on the fluid flow is meaningful. As Figure 11 shows, a series of elements are added to the surface of a channel by a strategy borrows from fractal theory, which makes the morphology complex and mimics some real situation. The construction of the profile follows the step that half-sized elements will be created and placed at a half-height distance position away from the original one. Then the gap between two large elements is filled by small elements. This figure mimics some real profiles of the channel surface. As the above analysis shows, the flow of fluid at the interval of two elements is constrained; thus, what changes will be on the flow pattern and the permeability if the structure is simplified by keeping only two main summits on the channel surface. If paddings at the interval of two large elements are removed, the morphology of the channel surface will change largely. Take relative roughness as a tool to evaluate the morphological changes, the relative roughness will decrease by about 32% if small elements are eliminated.

To make the difference of results between original channels and simplified ones clearly, a large mesh size (about 300,000 grids) is used for the simulation. Figure 12 shows the flow pattern and permeability results of two kinds of channels, and situations including micro effect and not are also offered in this figure. It is obvious that if the micro effect is included, permeability increases a small number, which is in accordance with the expectation. In original channels, which has labels ending with the number 1 in Figure 12, small vortexes can be found at the corner of each valley between elements, while larger vortexes occupy the position of small summits in simplified channels. From the contour of velocity, areas of different speed level can be judged. It is noted that regardless of the structure of channels, the velocity contours are exactly similar before and after the simplification. Furthermore, the permeability shows an impressive result that the simplified channels just have a little bit larger values than the original ones, which indicates that even the morphology of two channels is largely different, they can have the same permeability. The result also stresses the dominant role of main summits on fluid flow through a channel.

The above conclusion is achieved from ideal models. To furtherly verify it and try to put it into application, a real fracture profile obtained from a marble sample [69] is used for a test, as Figure 13 shows. This profile can be divided into two main summits and one subsummit. After simplification, the subsummit will be removed, and the permeability of two channels will be compared. Similarly, the micro effect is investigated. Approximately 400,000 grids are used corresponding to the resolution of the profile, and the length of the simulation domain is scaled according to the size of the profile.

Figure 14 shows the results of the channel before and after the simplification. It is clear that the permeability value under the slip condition is larger than the nonslip condition. After simplification, the velocity field is similar to the original one, despite the vanishing of apparent vortexes. Quantitatively, the permeability changes slightly due to the simplification of about 1%, which is definitely negligible.

4. Summary
By introducing micro effects into LBM simulation, an investigation on the fluid flow through microchannels is conducted. The result is also compared with the case that micro effects are neglected. Based on adding differently shaped artificial elements on the surface, the mechanism of flow affected by the surface structure has been studied.
With the change of density of elements, the permeability exhibits different value which has a specific relationship with the elements’ density. The relationship can be expressed linearly if the reciprocal of the permeability is chosen as the dependent variable. This result exists both under slip and nonslip flow conditions.
Differently shaped artificial elements exert unlike influence which is expressed by local flow patterns. The rectangle element leads to the lowest permeability of the channel and has the most developed vortex structure locally. The semicircle element and the triangle element follows at the flux-reducing effect and the size of local vortexes. It is convinced that these elements on the surface have a strong velocity-slowing effect on the near fluid and changes local pressure greatly. What should also be noted is that under the slip flow condition, which corresponds to the real flow condition in microchannels, elements show stronger performance on all the above effects.
Further researches on the morphology effects are conducted based on complex channels constructed by a fractal strategy and its simplified channels correspondingly. The result shows that even small elements locating at the interval of main summits can change the morphology of the channel largely; however, they only play a small part on the effect of fluid flow. The main summit dominates the flow both in macro- and microchannels. In this case, the morphology effect on the fluid flow of channel structures can be ascribed to the existence of large summits, and complex channels can also be simplified when making a permeability evaluation. An attempt on the channel created by a real rock fracture profile provides a preliminary reference.
Through dimensional analysis, a reference of physical quantities in the model this paper used is listed in Table 1, which is a representative case for microchannels or fractures. The corresponding lattice value for the simulation based on lattice Boltzmann method is also shown in this table.
Data Availability
Parameters used for the modeling and simulation in this work are included in the main text, and the data obtained from the simulation work that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
We thank Dr. Martin Haustein at TU Bergakademie Freiberg, Prof. Jili Feng, and Dr. Wenbo Gong at China University of Mining and Technology (Beijing) for fruitful discussions. We acknowledge the support from the National Natural Science Foundation of China (under Grant Nos. 51727807, 11772064, and 41472259), the National Basic Research Program of China (under Grant No. 2013CBA01504), the National Key Research & Development Program of China (under Grant No. 2016YFC080250504), and CAEP Foundation (under Grant No. CX2019033).