Resistance of velocity slip ﬂow in pipe/channel with a sudden contraction

A novel approach based on the local entropy generation rate, also named as the second law analysis (SLA), is proposed to compute and visualize the ﬂow resistance in mass transfer through pipe/channel with a sudden contraction component (SCC) at low Reynolds number ( Re ) featuring velocity slip. The linear Navier velocity slip boundary condition is implemented using the explicit scheme. At small Reynolds number, i.e. Re (cid:20) 10 : 0, the ﬂow resistance coefﬁcient of the SCC, K SCC , is found to be a function of the dimensionless velocity slip length L (cid:3) slip and Re (cid:0) 1 , and gradually increase to a constant value at contraction ratio R area (cid:21) 8, reaching a formula K SCC = ( 0 : 4454 L (cid:3) 3 slip (cid:0) 1 : 894 L (cid:3) 2 slip + 2 : 917 L (cid:3) slip + 8 : 909 ) = Re . Over this range of Re , the equivalent length of the ﬂow resistance is almost independent on Re , while out of this range, the equivalent length increases monotonically with Re . Moreover, the dimensionless drag force work around the SCC is negative and reaches a minimum at a critical L (cid:3) slip . The SLA reveals that the regions affected by the SCC mainly concentrate around the end section of the upstream pipe/channel rather than the initial partition of the downstream section reported in large Re turbulent ﬂow, and this non-dimensional affected upstream length increases with L (cid:3) slip . The ﬂuid physics are further examined using SLA to evaluate the energy loss over the entire domain, decomposed as the viscous dissipation inside the domain and the drag work on the wall boundary.


I. INTRODUCTION
The determination of the flow resistance in low Reynolds number (Re) mass transfer through micro and nanoscale pipes and channels with sudden contraction components (SCC) is a vital issue in the process of design and fabrication of efficient microfluidic devices (i.e., micro-electro-mechanical systems 1 , heat exchangers 2 and cooling of electronic chips 3 ). For fully developed laminar pipe flow without velocity slip, it is known that the Darcy friction factor (flow resistance coefficient is proportional to the friction factor in a straight pipe/channel) is 64/Re [4][5][6][7][8] . This value was also reported in the well-known Moody diagram, stressed by Moody as independent on the wall relative roughness in laminar flow 9 . Kim 10 reported that this value can be predicted by the classic macroscopic theory in microscale flow and concluded that scale effects do not influence the flow resistance coefficient.
However, a number of experimental measurements of the flow resistance in the laminar flow regime do not agree with each other and even show significant discrepancies 11 . After performing a series of experiments, Qu et al. 12 observed that the friction factor of laminar flow in microchannels with various shapes and dimensions is larger than the one obtained from the macroscale theory. On the contrary, the friction factor was found to be less than its macroscale counterpart by other researchers [13][14][15][16][17] . There have been no conclusive explanations of these significant discrepancies. Some works attributed these to the flow separation 18 , the dependence of viscosity on the channel size or other physical effects which become significant in small-scale ducts 19 and so on. We believe that the following factors are also involved: Firstly, the non-slip boundary condition on the wall becomes less valid as a significant velocity slip has been a) Electronic mail: xuerui.mao@nottingham.ac.uk observed [20][21][22][23] when the characteristic length of the channel is reduced. This is because the fluid slips when the molecular attraction between the fluid and the solid surface is reduced as the free surface energy of the solid becomes very low. Watanabe et al. 24 conducted experiments on water flow through square and rectangular ducts and pipes with highly water-repellent walls producing velocity slip boundary conditions and observed the maximum drag reduction ratio could reach up to 20% for square ducts and 14% for pipes. This observation offers potential explanations on the overestimation of the friction factor at extremely small scales by the macroscale theory aforementioned in the literature.
Also, the friction factor of laminar flow in micro-level ducts has been reported to be dependent on the relative roughness of the solid wall 25 , contradicting the Moody diagram 9 . Gloss et al. 25 conducted experiments to test laminar flow in microchannels with length 130 µm, height from 20 µm to 400 µm and relative roughness ranging from 0 to 1. Their outcomes demonstrate that the classic macroscale no-slip theory underestimates the friction factor because an increased dissipation rate near the roughness element causes more pressure loss compared with smooth wall surfaces.
Flow resistance is also dependent on the geometry variations, such as an SCC, as has been investigated in macroscale flow with applications including injection moulding, extrusion, thermoforming and polymer melt processing. A wide range of investigations including theoretical analyses 26 , numerical simulations 27 and experiments 28 were conducted for flow in an SCC with various contraction ratios. It should be pointed out that almost all these studies were on the effects of contraction ratio or wall slip on the flow behaviors (i.e., size and intensity of vortex, vortex shape and streamlines) of non-Newtonian fluids. No studies on the flow resistance coefficient of contraction pipe/channel flow of Newtonian fluids with a wall slip boundary condition have been performed to the best knowledge of the authors.
The flow resistance coefficient is commonly computed from pressure loss by solving the Bernoulli or Navier-Stokes (NS) equations. Celata et al. 29 observed that the viscous heating could be applied to determine the resistance coefficient of microchannel no-slip flow and validated this in their experiments. Naterer et al. 30 reported that the loss of available work could be directly related to the entropy generation by comparing the viscous dissipation term in the mechanical energy equation with the volumetric rate of entropy generation equation firstly proposed by Bejan 31 . The method based on the entropy generation is called the second law analysis (SLA), referring to the second law of thermodynamics. The SLA has been widely used in the optimal design of fluid engineering systems 32 , analyzing the flow resistance for no-slip flow, where the drag work on the surface is zero 33 . Compared with the common NS method, the SLA gives a direct physical interpretation of the flow resistance by revealing the distribution of the energy loss over the entire flow domain, as well as the drag work along the surface.
The main objective of the present study is to quantitatively examine the flow resistance coefficient and reveal the physics in flow through a microscale SCC with a linear Navier velocity slip boundary condition (SBC). In Section II, we illustrate the details of the mathematical formulation of the linear SBC and the analytical results of the resistance coefficient of flow in a plane parallel channel from the view of exact solutions of the NS and SLA approaches, followed by validations by numerical simulations. Corresponding simulation details, flow resistance coefficients and two factors causing flow resistance (the wall drag force work and the fluids viscous dissipation) for flow through an SCC are then described in Section III. Finally, conclusions are drawn in Section IV.

II. RESISTANCE OF VELOCITY SLIP FLOW IN A 2D PLANE CHANNEL
A. The Navier velocity slip wall boundary condition The incompressible steady Newtonian flow through a plane parallel channel as shown in Fig.1a is governed by the NS equations 34 : where U is the velocity vector, U ⊗ U = U • U T , p is the pressure, ρ is the density and µ is the dynamic viscosity. Considering that the mass flow rate is small enough to activate the partial boundary velocity slip, the normal component of U on the wall is zero (see Fig.1b ), and the tangential component on the wall, denoted as U wt , is a function of the stress vector tangent to the wall τ wt : where τ wt represents the l 2 norm of τ wt , the negative signs indicate U wt and τ wt are in opposite directions, and S l and S nl are the linear and nonlinear slip coefficients, respectively 35,36 .
The present work will focus on the linear SBC, while the nonlinear counterpart commonly used for non-Newtonian fluids where the viscosity is dependent on shear rate 37 can be implemented similarly. For Newtonian flow, the fluid motion near the wall can be assumed as local Couette flow, therefore, the tangential wall stress vector τ wt is where µ is the dynamic viscosity, C is the central point of the mesh cell next to the wall, Z is the central point of a boundary face, δ CZ is the vector between this two points and δ CZ,⊥ is the projection of δ CZ in the direction of n. U w and U C are the velocity vector at Z and C, respectively and U wt and U Ct are their tangential components . I is the third-order unit matrix and n is the normal vector to the boundary face. Substituting Eq. (4) into Eq. (3a), the linear SBC can be obtained as where L slip (L slip = S l µ) is the velocity slip length 38,39 .
The tangential velocity of the fluid in the vicinity of the wall representing the Navier slip boundary conditions can be modelled via the following explicit scheme: where the superscripts i and i −1 denote values of corresponding quantities at the current or former time step, respectively. In the explicit scheme, relaxation is commonly implemented to improve the stability and convergence property when using large spatial or temporal intervals.

B. Resistance coefficient derived from the exact solution of the NS equations
The velocity profile for fully developed laminar flow through a plane parallel channel as shown in Fig.1a with linear SBC can be obtained from Eq. (1) and Eq. (2): where H is the height of the channel, and this velocity profile is in agreement with references 40 . Then the average velocity can be calculated from In general, the pressure loss over a streamwise length L for all types (laminar or turbulent, circular or noncircular, smooth or rough surfaces, horizontal or inclined pipes) of fully developed internal flows can be represented as 41 where K is named the flow resistance coefficient, which measures the efficiency of the mass transfer.
Considering that the pressure gradient is constant, or ∆p/L = dp/dx and combining Eq. (8) and Eq. (9), K for fully developed laminar flow in a plane parallel channel with linear SBC can be obtained as where Re is the Reynolds number based on the average velocity and the channel height H.
Similarly the flow resistance coefficient in a fully developed laminar pipe flow with SBC can be obtained: The first factor on the right of Eq. (10) and Eq. (11) represents the friction factor. Clearly in Eq. (11) when the no-slip condition is imposed (L slip = 0), there is K = 64/Re, a well known formula for the macro-scale laminar pipe flow. Moreover, the wall friction coefficient ( f ) can be defined as the ratio of the shear stress and fluid velocity at the wall, and subsequently there is f = µ/L slip from Eq. (4) and Eq. (7). This indicates the wall friction coefficient is inversely proportional to the velocity slip length at constant fluid viscosity, in agreement with the experimental work 42 .

C. Resistance coefficient calculated from SLA
As an alternative to the definition, it has been proposed that K can be obtained from SLA for non-slip flow 33 : where T is the temperature,Ṡ irr, D is the generation rate of entropy due to thermodynamic irreversibility, and TṠ irr, D corresponds to the rate of the exergy loss according to the Gouy-Stodola theorem 43 .
When the SLA method is used to calculate K for the slip flow from the view of the balance of energy, the skin friction drag does negative work on the fluid (which is zero at the no-slip condition). In other words, the head loss of fluid is caused by two factors: viscous dissipation and skin friction drag. Therefore, the definition of the resistance coefficient should be adjusted to take into account this term as where W drag represents the contribution of the wall skin friction on fluid, and it can be calculated from Eq. (4) and Eq. (7) as Bejan 44 proposed that the volumetric rate of entropy generation consists of the volumetric rate of the heat transfer irreversibility and the fluid friction irreversibility. For steady, isothermal and incompressible channel flow without or with linear SBC considered in the present work, the local entropy generation rate is only from the fluid friction irreversibility and can be computed froṁ The overall entropy production rate can be obtained by integrating the local entropy generation rate: Substituting Eq. (14) and Eq. (16) and into Eq. (13), the flow resistance coefficient is the same as Eq. (10) derived from the NS equations. Similarly the resistance coefficient in a fully developed laminar pipe flow with SBC can be obtained and it also matches that from the NS method (see Eq. (11)).

D. Validation
Numerical simulations are performed to validate the derivations in Section II C. The plane parallel channel shown in Fig.1a with L = 30 mm, W = 0.2 mm and H = 2 mm is adopted. A circular pipe with length and diameter L = 30 mm and D = 2 mm, respectively, is also considered. The density of water is ρ = 997.05 kg/m 3 , and its dynamic viscosity is µ = 8.90 × 10 −4 Pa · s. The average velocity for no-slip or slip flow is U m = 0.50929 m/s and the slip length is L slip = 10 −4 m. ∆t = 10 −4 s was chosen as the time step to ensure that the maximum Courant number is less than 1. The coupling of pressure and velocity is addressed by the PISO algorithm.
The streamwise velocity and coordinate system are nondimensionalised by average velocity U m and channel height H, respectively. Fig. 2 shows the dimensionless streamwise velocity profiles along the height direction of the plane channel (velocity distributions in a circular pipe are not shown here for conciseness). The simulation results are in excellent agreement with the analytical outcomes. Furthermore, the explicit scheme leads to well agreed velocity distributions. The resistance coefficients for fully developed slip flow in the plane parallel channel and the circular pipe calculated from the NS equations and the SLA approach are in good agreement with each other, as shown in Table I. This validates the analytical derivations presented in Section II C and demonstrate the accuracy of the aforementioned explicit scheme of the linear SBC.

A. Problem descriptions
After investigating the efficiency of mass transfer in a straight channel, the flow resistance coefficient through an SCC as schematically illustrated in Fig. 3 will be examined using both the NS and the SLA approaches. The total flow resistance consists of the loss in the upstream large channel, downstream small channel and the SCC.
By definition, the flow resistance coefficient K SCC induced by the SCC can be written as where ∆p up,• and ∆p down,• denote the pressure loss of fully developed laminar flow in the upstream and downstream channel without SCC, respectively, and they can be determined from Eq. (9) and Eq. (10). U 2 is the average velocity at the outlet and will be used as the characteristic velocity, and ∆p is the pressure difference between inlet and outlet. K SCC can also be obtained from the SLA method through re-writing Eq.
where TṠ irr, D,up,• and TṠ irr, D,down,• represent the entropy production rate of flow that is not affected by the SCC in the upstream and downstream channel, respectively. These terms are calculated from Eq. (16) when considering no-slip or linear SBC. W D,up,• and W D,down,• are the work contributed by wall skin friction, and can be obtained from Eq. (14). A is the area of the cross section of the downstream channel. Also, the overall flow resistance coefficient K total of this SCC could be determined directly from Eq. (13).
The height of the downstream channel H 2 is set to 2 × 10 −5 m and used as the characteristic length, resulting in L * slip = L slip /H 2 . H 1 /H 2 represents the contraction ratio and will be denoted as R area . The upstream and downstream lengths of the channel are L up /H 2 = 25 and L down /H 2 = 50, respectively, which have been found to be long enough to accommodate the flow affected by the SCC. Moreover, the wall shear stress and the drag force work will be nondimensionalised by ρU 2 2 and ρU 3 2 H 2 2 , respectively. The time step for each case is adjusted to ensure that the maximum Courant number is less than 1. The flow resistance coefficient K SCC of the SCC computed from Eq. (18) and Eq. (19) is shown in Fig. 4. Unlike the resistance coefficient in the macro turbulent flow which is only dependent on the cross section ratio and insensitive to Re, K SCC in the present micro channel drops rapidly with increasing Re. Fig. 4b shows the effect of the contraction ratio R area on the flow resistance coefficient K SCC . At prescribed Re and L * slip , one can find that K SCC increases slightly with the contraction ratio R area . Through data fitting, the resistance coefficient over low Re (e.g., Re ≤ 10) can be modelled as where coefficients (a, b, The dimensionless equivalent length L * equi increases with respect to L * slip , but is insensitive to Re when Re is small enough as shown in Fig. 5, where R area is fixed at 4:1 and the numerically computed K SCC is adopted. This can be also demonstrated analytically by substituting Eq. (20) into Eq. (21), which results in a Re-independent L * equi . However at higher Re (e.g. Re ≥ 10), the analytical fitting equation underestimates the resistance coefficient K SCC and consequently L * equi , as confirmed by the numerical results in Fig. 5. Variations of L * equi at other R area are similar and are not shown here. Above we reported effects of Re, L * slip and R area on the flow resistance coefficient K SCC and dimensionless equivalent length L * equi . In the following sections, the fluid physics will be explored and efforts will be focused on the two sources of the flow resistance, i.e., drag force work and viscous dissipation loss.

C. Wall shear stress vector and drag force work
The dimensionless drag force work of the SCC, W * D,SCC , increases with Re number for velocity slip flow as shown in Fig. 6. Furthermore, there is a critical dimensionless velocity slip length (L * slip,cri ≈ 0.4 for all the four R area studied) where W * D,SCC reaches its minimum. One should note that W * D,SCC is almost always less than 0 (the reason will be illustrated later), indicating that the drag force work tends to reduce the flow resistance coefficient of the SCC according to the definition in Eq. (19). amplitude of the upstream and downstream wall shear stress τ * x for slip flow are far less than |τ * x | for no-slip flow at a prescribed Re, and this is in agreement with the experimental observation 42 .
The dimensionless drag force work per unit area w * d for the SCC with R area = 4 : 1 is plotted in Fig. 8 to explain the negative value of W * D,SCC mentioned above. In terms of the dimensionless tangential wall shear stress work of the SCC in Fig. 6, W * D,SCC can be rewritten as The region of the downstream channel affected by the SCC is very small according to Fig. 8 and w * d,down is almost constant rendering the third term of the right side of Eq.  Fig. 3).

D. Disturbed region and local entropy generation rate
To examine the range of regions affected by the SCC, a nondimensional surface integral of the local entropy production rate is defined aṡ where A is an arbitrary cross section area.Ṡ irr,D,• represents the local entropy production rate of flow without SCC and its surface integral can be formulated based on Eq. (16). The non-dimensional length of upstream and downstream channel domain affected by the SCC can be defined as  Fig. 9, L * up,dist and L * down,dist are almost independent on Re and increase with L * slip . In addition, L * up,dist rises conspicuously and L * down,dist remains almost constant with increasing R area .
Contours of TṠ irr,D are shown in Fig. 10 from which we can observe the viscous dissipation loss of fluid at every point in the flow field and the effects of Re, L * slip and R area . Comparing Fig. 10a with Fig. 10b, the viscous dissipation loss further concentrates on the corner as marked by the green circles in these two subfigures at larger Re. One interesting observation is that an "ear" type dissipation loss region appears at the beginning section of the downstream channel for slip flow, as marked by the green circle in Fig. 10c.
Moreover, the local entropy production rate of fluid in the SCC consists of six sections, mainly including upstream and downstream channels and 90 • corners, etc (see Fig. 10d). Conspicuously, TṠ irr,D for flow of the upstream and downstream channel affected by the SCC increases from the center line of the SCC and peaks at the wall (see the arrow 1 and 6 in Fig. 10d), suggesting that the viscous dissipation loss minimizes at the center line and enlarges along the height direction of upstream and downstream channels. TṠ irr,D for the 90 • corner 4 increases along the direction departing from the corner, and the variation trend is opposite to the other 90 • corner 5 where the local entropy generation rate or flow resistance is the maximum in the whole flow field. The growth gradient for TṠ irr,D at these two corners rises dramatically at larger Re. An interesting point is that there is a "source" (surrounded by the green circle in Fig. 10d) of TṠ irr,D , from which the flow resistance increases along with opposite directions ( 2 and 3 in Fig. 10d).

IV. CONCLUSIONS
This work addresses the flow resistance coefficient in flow passing through a pipe or channel with or without a sudden contraction component at low Reynolds number conditions with a slip wall velocity boundary condition which is implemented using the explicit scheme.
The resistance coefficients of the fully developed channel and pipe flow at low Re with linear Navier velocity slip boundary condition are solved with the NS equations and the second law analysis (SLA) approach which is extended to determine and visualize the flow resistance loss of velocity slip flow after considering the work contributed by wall drag force.
Then, numerical simulations are also performed for these two geometry models and the obtained resistance coefficients are in good agreement with analytical outcomes. Effects of Re and dimensionless velocity slip length (L * slip ) on the flow resistance coefficient of laminar flow through the SCC with various contraction ratios (R area ) are revealed based on the aforementioned NS equations and SLA method after conducting a series of numerical simulations.
At Re ≤ 10, the flow resistance coefficient for the SCC, denoted as K SCC , increases monotonously with L * slip and 1/Re and nonlinearly with R area , and at R area ≥ 8, the resistance coefficient can be modelled as K SCC = (0.4454L * 3 slip − 1.894L * 2 slip + 2.917L * slip + 8.909)/Re. Physically, the flow resistance of fluids results from two factors: wall drag force work (equal to 0 for no-slip flow) and viscous dissipation in the flow field related to the local entropy generation rate. These two terms are computed and analyzed individually to illustrate the fluid physics underpinning the flow resistance.
The amplitude of the negative drag work (W * D,SCC ) decreases monotonously at larger Re. Also, there is a critical velocity slip length at which the amplitude of W * D,SCC reaches maximum. In terms of fluid internal viscous dissipation, it minimizes at the center line and peaks at the wall of the upstream and the downstream channel. Moreover, the dissipation loss concentrates around the end section of the upstream channel of the SCC.
Both the dimensionless lengths of the disturbed region of the upstream and downstream channels are almost independent on Re and enlarge at larger L * slip . With increasing R area , the former rises remarkably while the latter remains nearly a constant value.