Xin Zhao,Chao Wei,Shi-hua Yuan

Science and Technology on Vehicle Transmission Laboratory,Beijing Institute of Technology,Beijing 10081,China

Abstract: A new slip velocity model based on molecular potential theory and macro-force analysis,which is applied in Couette flow with pressure gradient,is built up.The model is validated by later being introduced in hydrodynamic system to predict film distribution,which shows a good agreement with experimental data obtained from multi-beam intensity-based tests with Fe and Cu ball materials under accurate controlled temperature,load and different wall velocities.Results show that the slip length for Fe case is ignorable so it seems like no slip,but for Cu case,the slip length is large to make 20% discrepancy with no slip simulation and also behaves shear-dependently.Moreover,during the experimental cases when both Fe and Cu ball velocities rise from -133 mm/s to 1 330 mm/s,the slip velocity changes its direction with entrainment velocity and thus contributes to first enhance and then diminish the hydrodynamic film,but due to slip length on Cu case varying largely than that on Fe case,the film from Cu case and from Fe case has a clear cross-point between ub=80 mm/s and ub=220 mm/s (ub is the ball speed).The results above support strongly that Cu surface will lead to stronger slip than Fe case because of its smaller solid-liquid interaction,and obviously slip will influence hydrodynamic characteristics prominently.

Key words: Slip,Couette flow with pressure gradient,hydrodynamic characteristic

Introduction

Although more than two decades of intense research on hydrodynamic slippage on various substrates for different flow systems proves advantageous application on biology,chemistry and engineer[1-3],most of them focus on microfluidic system[4-6].Microfluidic is chosen as a typical system for slippage research due to its microscale similar with slip length (quantifying slip level),which makes slip effect notable[7],on the other hand,it is easier to investigate single factor on slip in microchannel due to its simple Laminar flow in Couette cell and simple cross-section geometry shape of channel[8-9].

According to these researches,the slip length can be designed to reach several micrometers on hydrophobic patterned surface with hybrid (liquid-gas and liquid-solid) interfaces.The scale is large enough to make effect on flow characteristics in macroscale flow system.Many theoretical attempts to introduce giant slip for application in macroscale flow system has been tried,and the slip effect on hydrodynamic features has been predicted[10].Aurelian et al.[11]introduces Navier slip model with empirical value of slip length into Navier-Stokes (N-S) equation to simulate the performance of fluid bearings with patterned surface.The empirical slip length presents a tendency of slip effect on hydrodynamic but not an accurate prediction.Davis and Eric[12]believes slip origins from liquid-solid surface force and adds the force into N-S equation to obtain the flow behavior in patterned surface.The method is limited by complicate surface force calculation of structured surface,thus at present only some types of regular surfaces can be applied,such as fair[12]and mesh[13]surfaces.Recently combination of slip length calculation based on molecular dynamics simulation (MDs) and hydrodynamic features from N-S equation is proposed to build up an accurate and general hydrodynamic model with slip effect.MDs is limited to nanometers in space and nanoseconds in time due to computational impracticality.However,we find the micro fluid layer (MFL) method[14]has similar accuracy with MDs in nanoscale,and it can also extend to micrometer scale system,providing possibility to be applied in macroscale hydrodynamic system.

Many macroscopic flows between two surfaces can be considered as the general Couette flow,driven by moving walls and an applied pressure gradient.However,the MFL model developed in our previous work[14]is only valid for the simple shear flow,and it requires a modification of micro and macro forces to be applied to a general Couette flow.In this paper,a modified MFL model in Couette flow with pressure gradient (C-Pressure flow) is proposed and then it is combined with N-S equation to build a uniform hydrodynamic model for the sake of prediction of hydrodynamic characteristics considering slip effect in macroscale flow system.The modified MFL model can be used to analyze both Couette without pressure gradient (Couette flow) and C-Pressure flow,and results show the pressure gradient can actually influence slip.The hydrodynamic film calculated by the uniform hydrodynamic model shows a good agreement with point-by-point film measurement by Multi-Beam Intensity-Based test bench in two different cases whose ball materials are made of Fe and Cu,respectively,confirming the validation of the modified MFL model.With the ball speeds vary from-133 mm/s to 1 330 mm/s,the slip velocity changes its direction with entrainment velocity and thus contributes to first enhance and then diminish the hydrodynamic film.Specifically,because of different solid-liquid forces,the slip length for Fe case is only about 0.2 nm so it seems like no slip,but for Cu case,the slip length is about 2 nm,which makes 20%discrepancy with no slip simulation.And for Cu case,slip length also behaves obviously shear-dependently.

1.Theoretical method

1.1 MFL method concerning pressure function

The traditional MFL method in Couette flow[14]is generated only by the relative movable walls induced shear function.However,in most cases,the pressure should be considered because flow in confined region need external energy.The pressure gradient will change flow profile largely.This kind of mixed flow is called Couette flow with pressure gradient (C-Pressure flow).The schematic of C-Pressure flow is shown in Fig.1(a).This defines that the distance between two walls is defined asH.The inside liquid bulk is divided intonlayers along thickness direction.For simple description,the two walls are named as the 0thand the (n+1)thlayer,and distance between theithand the (i+1)thlayer is defined ashi+1.The width of wall is Δl,it is large enough to ignore effect of edge and small enough to be seen as microscale cell.The total number of molecules,N,in each layer is assumed to be equal and large enough,so for point A inithlayer,the fluid layer can be seen as an infinite plane.The force of molecule in point A inithlayer is analyzed in vertical to walls with molecular potential theory and in parallel direction with shear tension equilibrium method.Figure 1(b) presents force analysis in bothzandxdirection for flow cells from Couette and C-Pressure flow.In Fig.1 (b),Fwis the force from wall,F1is the force from liquid,pAis pressure at Point A.We can see two kinds of flow cells have the same force analysis inzdirection,but suffer from different forces inxdirection due to the C-Pressure flow considering the pressure effect.So we can make sure the force analysis between Couette and C-Pressure flow shares the same equation inzdirection.

Fig.1 (Color online) The model schematic of C-Pressure flow effects

The force equation of flow cell inzdirection is complex due to its origin from all the molecules in liquid and wall.In this way,it can be obtained by summing up the force that all molecules in this flow cell suffer from inzdirection.Assuming that the molecules in flow cell share the same force state (it is reasonable because the flow cell is small enough),therefore,the force analysis of molecule in point A is important for the equation of flow cell inzdirection.After we choose the cut-off radius to be 2.5σ[15],the resulting force expression for point A alongzdirection is expressed as follows

at whichFw,jis the force acted on theithlayer from wall,j=1 and 2 represent the force from one and other wall andFl,krepresents the force from thekthfluid layer.

Specifically,forFw,j,it can be calculated in the force form of single molecule with an infinite wall surface,because wall molecules are regarded as equal distribution[16]and the molecules on the surface play a leading role.While for the force from thethkfluid layer,Fl,k,is more complicated because the force from each molecule at different positions must be considered.But as the layer is very thin,the force from layer acted on point A can be regard as force from plane.Thus the expressions forFw,jand forFl,kinkthlayer[17]are given as follows

whereεwandεfare energy parameter of wall and fluid molecules, respectively. When wall characteristics change,the energy parametervaries together.σwandσfare distance parameter of wall and fluid molecule,ρwis the wall density andρf,kiskthfluid layer density,which can be expressed by

whereMk-1andMkrepresent the molecules mass of the (k-1)thandkthlayer,respectively.Srepresents the area for a layer.wherer0is the distance between wall andithlayer,,andrkis the distance betweenithfluid layer and another layer,for the distance between theithand (i+2)thlayer,ri+2=hi+1+hi+2,and also the restricted conditionsfor the mass conservation as well.

Therefore,based on abovementioned analysis,the fluid layer density can be calculated.And the related viscosity for each fluid layer,η1,η2,… ,ηi,η1,η2,… ,ηi, … ,ηn,can be calculated according to relationship between density and viscosity in Literature[18-19].

whereη0andρ0are bulk viscosity and bulk density,pis the pressure.The viscosity has its limitation value,ηt,andz0is usually equal to 0.6 here.Because the density can be too large based on equations above,so the viscosity has its upper limitation,and therefore we useηtfor density that reaches 1 .1ρ0.

And the force inxdirection can be directly deduced from Fig.1(c).Therefore,for a fluid cell which contains point A,its equilibrium stress tension form is given as follows

where pressure gradient in the direction from inlet to outlet is regarded as positive.

Meanwhile,the force of each layer in shear direction is balanced,as illustrated by (for example,ithlayer)

And the boundary conditions are supplementary,V1=u0,V2=un+1.

When ∂p/∂x=0,the function only from Couette flow and Pressure gradient disappears.Equation (5)reduces to

According to the modified MFL model proposed above,the results of density and velocity for each layer under condition of C-Pressure flow and then the distribution of them can be finally deduced from the force analysis for flow between two walls.

1.2 Definition of slip length

As the velocity of each fluid layer is separately distributed and the first fluid layer keeps distance with wall,the slip length,b,can not be obtained from extrapolating the velocity profile at the fluid-solid interface.Therefore,slip length is defined by extrapolating the velocity profile from the position in the fluid (with first fluid layer distance from the first solid layer) to where the velocity would vanish within the solid according to Literature[20].The schematics of slip length for C-Pressure flow are shown in Fig.2,respectively.

Fig.2 (Color online) Schematic of slip length deduced from velocity distribution

For Couette flow,the slope fitted by the middle part of the velocity profile is chosen as the slope at the first point because the slope (shear rate) is theoretically the same at any positions for Newtonian fluid.So we first calculate the slope by linear fitting,and then apply this slope value at the first point to get the imaged no slip plate.The distance between imaged no slip plane and the wall is the slip lengthb(Fig.2(a)).For Pressure gradient induced flow (also for C-Pressure flow),as the slope (shear rate) is different at each position.So we use the slop fitted by the first and second points near wall to determine the imaged no slip plane (Fig.2(b)).

1.3 Slip on Couette and C-Pressure flow

We conducted some simulations for various pressure gradients which are easily reached in flow system to confirm pressure can play an important role in slip phenomenon in flow field.The simulation parameters are given in Table 1.Fe and glycerin are chosen as two walls and fluid.The up wall moves with the velocity of 5 μm/s while the down wall stays static.The pressure gradient is set at -500 Pa/m,-250 Pa/m,0 Pa/m,250 Pa/m,500 Pa/m,respectively.

Figure 3 shows the velocity profile between two walls under different presure gradient (PG) values.Each point corresponds to velocity of each fluid layer.The region between first and last points is defined as bulk phase.Clearly,the velocity profile goes up as the PG rises.For Couette flow (PG=0 Pa/m),the velocity profile is linear in bulk phase.The slip lengths in up (b′) and down (b) wall are both equal to 2.6σ,but for C- Pressure flow,for example,in the case of PG=500 Pa/m,the slip lengths for up and down walls are 4.1σand 1.0σ,respectively.To be specific,the positive PG values tend to enhance the slip amount in the up wall but diminish the amount in the down wall,thus the velocity profile is convex to outlet.It is because positive PG tends to increase the velocity on the up side but decrease that on the down side,which eventually make slip length rise on the up and decrease on the down.The negative PG values perform inversely and its velocity profile is concave from outlet.The velocity profile becomes more non-linear as the effect of PG increases.

Table 1 Simulation parameters

As shown in Fig.4,with the rise of PG,the slip length on the up wall increases from 1.0σto 4.1σ,while that on the down wall decreases from 4.1σto 1.0σ.In the case of the same PG value acting on different directions,for example,PG=500,-500,their slip length are both 4.1σand 1.0σ,but in opposite walls because the pressure gradient contributes to increasing velocity in its direction and decreasing that in opposite direction.The slip amount,defined by (bc-p-bc) /bcto quantify pressure effect on slip,wherebc-pandbcare slip length for C-Pressure flow and Couette flow,can reach 60%between Couette flow and C-Pressure flow with PG=500 Pa/m which means the slip length in case PG=500 Pa/m is 1.6 times as that in case Couette flow (PG=0 Pa/m).

1.4 A uniform hydrodynamic model from Navier-Stokes and MFL methods

The different slip lengths from upper and down surfaces are introduced into Navier-Stokes equation[11].

Fig.3 (Color online) Velocity distribution for different PGs

Fig.4 Slip length in different PGs

whereb1andb2are slip length of down and upper surfaces,respectively.Equation (9) is solved by finite difference method.So the physical amount ofp,ρ,h,uandvare dispersed.And the slip length dependent on the pressure and distance between two surfaces also varies in different positions.Therefore,slip length can be calculated separately according to local shear stress with (ui+1-ui) /hiand pressure gradient with (pi+1-pi) /dlin each position,where dlis the step length along surface.The flow chart in Fig.5 describes the calculation process of the uniform hydrodynamic model.To begin each circle,the values of load,wall velocities inxdirection (U,V) and fluid parameters are first input in Hydrodynamic section to obtain the pressure and film distribution between two walls.Then the pressure,film and the movable velocity of two surfaces (U,V) are input into MFL section which has installed with parameters for certain fluid-solid system.As the pressure and film are quite different at each position,the velocity profile at each position is also quite different.Then the slip length is deduced from velocity profile at each position and is transferred to Hydrodynamic section as new velocity boundary condition to renew pressure and film values.Then,new velocity boundary,pressure and film are sent to MFL section.The circle will end until pressure value from two adjacent circles keeps almost unchanged.After the process ends,the pressure,film,velocity distribution and other hydrodynamic characteristics will come out.

Fig.5 Flow chart

2.Experiment section

2.1 Test bench

The multi-beam intensity-based test bench can measure point-by-point film thickness at any positions between ball and disc.Each film thickness calculates a slip length at the corresponding position.The experimental condition can easily adjust variation for single variable.

The flow system in multi-beam intensity-based test bench is made up of five systems: ball-disk motion system,motion control system,loading system,data acquisition system and friction torque acquisition system.In the ball-disk motion system,the metal ball used was 0.0254 m in diameter and was controlled by one electric motor.The glass disc was coated with a layer of chromium (Cr),and another electric motor was connected to the torque sensor.The ball and the glass disc locate vertically and can be driven independently (motion control system).The external load (as shown in Fig.6) can be accurately added vertically on the ball based on an improved loading system,which uses a tensiometer to control the load.The gap between the ball and disc is filled with oil,and a small contact zone occurs under an external load.The image of the contact region can be magnified by a microscope and observed through a CCD (data acquisition system).The friction torque can be measured.Otherwise,the whole system is exposed in air to stay steady temperature.

Fig.6 (Color online) Multi-beam intensity-based test bench

In our experiment,both Fe and Cu ball are tested under different rotation velocities.The point-to-point film between ball and disc surfaces in the contact region and the friction torque are measured.

2.2 Experimental procedures

The test cases are conducted as following description: the metal ball (Fe and Cu materials are both used) increases velocity gradually from -133 mm/s to 1 330 mm/s (ub),while the disc is fixed at a velocity of 133 mm/s (ud).The load is set steadily at 15 N in these two cases.The interference figure and friction torque can be measured simultaneously.The related film thicknesses are calculated according to the principle in the literature[23].Each test will be repeated five times to obtain the average value and error.In the test,the glycerin is employed because it is easy to use in simulation; its properties and other material parameters are shown in Table 2.The temperature is controlled at 20°C in all experiments.

Table 2 Parameters of materials

3.Results and analysis

3.1 Slip phenomenon in the interferograms figures

Figure 7 give the interferograms image of the contact region in Fe and Cu ball rotation velocities and steady disc velocity (133 mm/s),respectively.Compared Figs.7(a),7(d),their theoretical entrainment velocity,u= (ub+ud)/2,is equal to zero,thus the interferograms in Fig.7(a) are circle,but the hollow in Fig.7(d) faces to right,which means that the real entrainment velocity also faces to right.A reasonable explanation is most probably slip occurrence on the ball surface (Cu).But whether slip also happens on the disc surface,it can not be judged here.Even though it does happen in ball surface,the amount of slip is quite small.The area of hollow in Fig.7(c) is larger than that in Fig.7(f),which reflects a larger entrainment velocity.

Fig.7 Interferograms with different entrainment velocities (Fe and Cu ball materials)

3.2 Validation of the uniform hydrodynamic model

The multi-beam intensity-based test can obtain the film thickness at any positions between ball and disc.In the meantime,the uniform hydrodynamic model can calculate the film thickness distributed within calculated region,which contains definitely contact region.Therefore,there are theoretically countless film thickness value obtained in each experiment and simulation.Although pressure difference in different positions can influence film thickness in different positions,the film thickness at the center of contact region (center of circle) always shares the same pressure due to the external load keeping unchanged in all experiments and simulations.Therefore,the Fig.8 gives the measured and simulated results of film thickness at the center of contact region with the variation of roller velocity for different ball materials(Fe,Cu) and one group no-slip simulation data.

Fig.8 (Color online) Film thickness at different ball velocities

The experimental data of Fe ball is close to both MFL-hydrodynamic simulation with Fe material and no slip hydrodynamic simulation,while the experimental data of Cu ball is merely close to MFL-hydrodynamic simulation with Cu material.At this step,it is hard to judge whether slip happens on Fe material because both MFL-hydrodynamic simulation with Fe material and no slip hydrodynamic simulation seem to meet Fe experimental data very well.Therefore,Fig.9 gives the detailed simulation about slip length of two types of ball-disc systems.Slip length of glass side (Cr) is also listed.It is clear that slip length in Cu surface (about 2 nm) is much larger (almost 10 times) than that in Fe and Glass (Cr)surfaces (less than 0.2 nm) between two systems.Thus,the slip can be ignored reasonably if the surface is made of Fe and Glass (Cr).This can be used to explain the reason that the traditional no slip hydrodynamic model can draw a close curve to the Fe ball-glass (Cr) disc system due to a too small slip amount (almost ignorable).Then we can find the slip lengths for four cases have similar tendency,it first decreases then increases slowly.The lowest point always lies inubbetween 80 mm/s and 220 mm/s.We believe it is because thatubchanges its direction,which changes from the opposite direction of disc to the same direction.This makes the shear rate experiences a tendency of first decrease and then increase.As the velocity of disc is steady at 133 mm/s,when the velocity of ball gradually varies from -133 m/s to 133 mm/s,the shear rate gradually decreases to zero,leading to a decreased slip length.After that,the shear rate keeps increasing again,making the slip length also increase.

Fig.9 (Color online) Slip length at different ball velocities(simulation)

The slip length can directly influence the velocity of fluid at the ball interface.The entrainment velocity,defined by (ub,s+ud,s)/2,whereub,sandud,sare fluid velocity at the ball and disc interface,is also influenced by slip length.Larger entrainment velocity can arise from larger hydrodynamic effect,which generally means larger film thickness.When the velocity of ball varies from -133 m/s to 133 mm/s,the slip velocity is the same with direction of entrainment,thus contributing to larger hydrodynamic effect.After that,the slip velocity changes its direction and diminishes the hydrodynamic effect.Therefore,we can see an obvious cross point between slip case (Cu)and no slip case at aroundub=133mm/sin the magnifying region Fig.10.

Fig.10 (Color online) Film thickness at different ball velocities

3.3 Shear dependent

We further give the figure of simulated slip length on Cu material varying with shear rate (Fig.11).In our case,the slip length is clearly shear-dependent,showing increase with shear rate,as the same shown in Literature[24].Interestingly,although the slip length is shear dependent,but it varies quite slightly from 1.8 nm to 2.2 nm with shear rate largely rising from 2×105s-1to 2×106s-1.This can be an explanation why in many experimental results represent slip is shearindependent.

Fig.11 (Color online) Slip length on Cu material at different shear rates (simulation)

4.Conclusion

In conclusion,a new slip velocity model in Couette flow with pressure gradient is built up.The model is introduced into hydrodynamic method to obtain the film which later is compared with film from multi-beam intensity-based test under two cases with two ball materials (Fe,Cu) under accurately controlled temperature,load and different wall velocities.Both experimental and simulated results strongly support that the slip length on Cu case is obviously larger than that on Fe case.Specially,in Experiment,when both Fe and Cu ball velocities rise from -133 mm/s to 1 330 mm/s,the slip velocity changes its direction with entrainment velocity and thus contributes to first enhance and then diminish the hydrodynamic film,but due to slip length on Cu case varying largely than that on Fe case,the film from Cu case has a clear cross-point with that from Fe case betweenub=80 mm/s andub=220 mm/s (Fig.10).In simulation,when we look detail of the slip length for each case under different ball rotation velocities,the slip length for Fe case is only about 0.2 nm,while for Cu case,the slip length is about 2 nm.The evidences above hold the truth that the smaller solid-liquid interaction contributes to larger slip length.And as for such small slip length,in most case,the no slip boundary condition can still hold true.Furthermore,we investigate the shear effect on slip length in Cu case using our MFL model.Although slip length seems to be shear-dependent,slip length varying from 1.8 nm to 2.2 nm when shear rate varies from 2×105s-1to 2×106s-1,we suspect that the variation of slip length is so small to be detected in experiment.This can be an explanation why in many experiments represent slip is shear-independent.

Acknowledgement

This work was supported by China Scholarship Council.