 In the previous lecture, we were discussing about the fundamentals of fully developed flow and developing flow in channels and let us recapitulate what we studied a little bit before proceeding further. So we considered fully developed flow in a parallel plate channel and considering that we started with the continuity equation and the momentum equations to describe the velocity profile, the fully developed velocity profile which is parabolic in nature. Then we calculated the dimensionless wall shear stress and dimensionless pressure drop and because the wall shear stress and pressure drop are related, their dimensionless forms are also related. The dimensionless wall shear stress is the Fanning's friction factor Cf and the dimensionless pressure drop is the Darcy's friction factor f and we derived that f is equal to 4Cf. This was what we discussed in the previous lecture and we will try to extend the discussions of the previous lecture for flow through circular channels, right. In microfluidics, rectangular channels are more common than circular channels but circular channels are also important in the context of micro needles, I mean which may have circular section. So it is important to understand the fundamental mathematical paradigm by which we describe the circular flow through circular tubes, the fully developed flow through circular tubes which is given by the Hagen-Poisouli law. So to understand that, let us appeal to this slide which is there. So in this slide you can see that we have drawn a cylindrical tube, I mean it is better not to call it a pipe because we are considering small dimension. So we may call it a tube, although like geometry wise there is no great difference, tube of radius r and we have Vr, Vtheta and Vz as the 3 components of the velocity. So here we are using the cylindrical polar coordinate system. The flow is steady so that del Vz del t equal to 0. The flow is fully developed so del Vz del z is equal to 0, it is just like del u del x equal to 0 for a rectangular channel, here the axial direction is the z direction, so the del Vz del z equal to 0 and because of axial symmetry there is no variation with respect to theta. So if you consider all these, then you will see that if you write the momentum balance, the z momentum equation, then the left hand side becomes 0, just like the left hand side in the x momentum equation becomes 0 for flow through a rectangular channel. Similar thing the left hand side which is the acceleration term becomes 0, then right hand side there is a balance between the pressure gradient term and the viscous term. Now the r momentum equation, so far we have discussed about the z momentum equation. In the r momentum equation you have basically all the terms because Vr is equal to 0, all terms are 0 that means del p del r which is the only term remaining that term is also 0. So because del p del r is 0 that means p is a function of z only, p is not a function of r. So that boils down to this equation, I am just going through it, I mean in a very straight forward manner because the same thing just in a different coordinate system we have discussed in details. So I just do not want to spend unnecessary time in describing the same physics for another case where only change is the description of the coordinate system but physically it is the same. So you can see that the z momentum equation the pressure is a function of z only and Vz is a function of r only. The only difference between the rectangular channel or the parallel plate channel and the circular pipe is the fact that the viscous term you express instead of d2 Vz dr2 it is 1 by r d dr r d Vz dr that is just a representation in the cylindrical coordinate system. So this is because Vz is a function of r, this is a function of r only, p is a function of z, so this is a function of z only. So function of r is equal to a function of z that may be possible only if it is a constant. So we take it up from here and go to the board to address this problem. So 1 by r d dr of r d Vz dr is equal to 1 by mu dp dz. What we do? So then each is equal to say c. So we can write d of r d Vz dr is equal to cr dr. We can now integrate it. So if we integrate it now r d Vz dr is equal to cr square by 2 plus c1. So you can see that d Vz dr is equal to cr by 2 plus c1 by r. Now if you consider the channel or the tube the center line is r equal to 0. So at r equal to 0 you must have d Vz dr finite. But this 1 by r dependence is killing that finiteness until and unless c1 is equal to 0. So d Vz dr must be finite at r tends to 0 which implies c1 equal to 0. So d Vz dr so Vz is equal to cr square by 4 plus c2. So we will apply the other boundary condition that at small r equal to capital R what is Vz? Vz must be equal to 0. So that means c2 is equal to minus cr square by 4. So that means you have Vz is equal to c by 4 into small r square minus capital R. So we can now express Vz in terms of V average. So what is V average? Integral Vz into 2 pi r dr 0 to r divided by pi r square right integral Vz df by the area. So c by 4 into 2 pi into r cube minus r square r dr 0 to r by pi r square. So this becomes c pi by 2 by pi r square into r4 by 4 minus r4 by 2. So this is minus cr square by 8 alright. So that means you can substitute that expression for V average in the expression for Vz and once you substitute that you get Vz by V average is equal to 2 into 1 minus small r square by capital R square. So this is the velocity profile. Again the parabolic velocity profile is evident from this analysis. What we will do, we will find out either the dimensionless wall shear stress or the dimensionless pressure drop. Let us find out the dimensionless pressure drop. So you can write V average is equal to minus cr square by 8 and c is equal to what? 1 by mu dp dz that is 1 by mu delta p by l. Now so you can write in place of delta p we are basically working out the same steps that we worked out in the previous class just the coordinate description is a little bit different the geometry is a little bit different. So this is minus hf rho g where hf is the head loss due to friction. So V average is equal to hf rho g by mu l into r square by 8. So hf is equal to 8 mu vl by rho g r square and you can write it in terms of the flow rate by writing V average is equal to q by pi r square. So 8 mu ql by rho g pi r to the power 4. Many times in engineering convention we express the result in terms of the diameter of the pipe rather than the radius. So if you substitute r equal to d by 2 then we get hf is equal to 128 mu ql by rho g pi d to the power 4. You can find out the friction factor f from here the Darcy's friction factor just in the same manner in which we calculated in the previous class for the rectangular channel and you will get f equal to 64 by Reynolds number. So what we have done we have summarized the description for fully developed flow for parallel plate channel and for circular geometry. So let us summarize the important findings as we get back to the slides. So far we have discussed about the fully developed flow and fully developed laminar flow. So this diagram I have purposefully put in this context. I mean you can recall that this is nothing but the Moody's diagram what you see here and we will try to understand that which regime is important for us in microfluidic applications. In microfluidic applications we are so in the x axis there is Reynolds number. In the y axis on the left side there is friction factor either f or cf whatever you want to plot they are related and the right hand side it is either the roughness by the diameter or the diameter by roughness in this particular case it is diameter by roughness that is plotted. So either of the diameter by roughness or the roughness by diameter. Then basically we are operating in microfluidics in the low Reynolds number regime and in the low Reynolds number regime you have laminar flow. So this curve cf so this figure is for a circular geometry. So cf f was 64 by re so cf is f by 4 so 16 by re that is this. So in the log plot it is a straight line okay. So this is the regime in which we are operating and you can see that here the friction factor appears to be independent of the surface roughness. As you go to the high Reynolds number regime you can see that for different values of the relative roughness you are getting different values of the friction factor. So the friction factor is depending on the roughness but for fully developed laminar flow the roughness is apparently having no role to play because all the perturbations induced to the flow because of roughness are attenuated. Perturbations do not amplify in a laminar flow. So that makes it independent of the roughness but still the analysis does not take into account the length scales of the roughness. So the length scales of the roughness themselves can be important as compared to the system length scale when you are talking about microfluidic systems. So it is worth investigating that what is the role of roughness for laminar flow through I mean micro channels of whatever geometry that is the first remark that we have put. The second remark is that you see hf is proportional to muql by rho g pi dh to the power 4 hydraulic diameter to the power 4 in this case the hydraulic diameter is the diameter itself. So that means that the head loss for a given flow rate is inversely proportional to the fourth power of the diameter. So if you want to have a particular flow rate then if you reduce the diameter from 1 meter to 1 micron see 10 to the power 6 to the power 4 that is 10 to the power 24 times the head loss will increase. So you will require tremendous pumping power to drive the flow. So you cannot really achieve that high a flow rate in a pressure driven microfluidic system which you can achieve in a macroscopic system you have to compromise with the q because to achieve the same q you require to have a huge driving pumping power therefore you compromise on q but sometimes it is a blessing in disguise not in all of applications you require a large q. For example you are administering a medicine or a drug to a patient in a very minute and precision quantity. So there you do not require a high flow rate there in fact a precise delivery of the drug is what is important. So high flow rate may appear to be an important and useful thing but it depends on the context where other factors may be more important than high flow rate. So these are the 2 important observations that we have made from our analysis so far. Now whatever we have discussed so far is something which we normally covered in the undergraduate first level course of fluid mechanics but we want to go beyond that in this particular course. So now what we will see is that like we may still consider steady flow but instead of a parallel plate channel it can be a channel where the width is not infinitely large as compared to the depth. So that means the width and depth are comparable so how do you get a solution for such a problem. So we will try to address that I would partially work out this problem and then leave the remaining part of the problem as a homework to you. So we will work it out partially. So we will first make a schematic to understand that what we are talking about. So let us say that you have a channel of cross section like this. Let us say this is y axis, this is z axis and perpendicular to the plane of the board is the x axis and the pressure gradient is applied along the x direction. So there is a dp dx. Let us say that this is a by 2, this is a by 2, this is b by 2, this is b by 2 and we are interested to find out the velocity, the wall shear stress, the dimensionless wall shear stress or the dimensionless pressure drop, the standard parameters which we are interested. Now let us write the governing equation first before going for a solution and to write the governing equation I mean we can do a similar analysis like the parallel plate channel. Still we are assuming that it is fully developed. Let us write the x momentum equation. Let us write the x momentum equation because it is fully developed. The left hand side is 0, steady and fully developed. The right hand side, so let us write the viscous term. Now there will be 2 terms, one for variation along y, another for variation along z. For the parallel plate channel case this term was not there. But for the channel where the width and I mean these dimensions height and width are comparable these 2 terms may be of the same order. So they have to be each has to be written – dp dx. So clearly this does not have such a simple solution as that of the parallel plate channel case. The parallel plate channel case for a fully developed flow reduced to a ODE but here it still is a PDE. Now when this problem first came to engineers, see engineers traditionally have always tried to come up with a working solution and the working solution may be good or bad but it is some solution and sometimes it is important to get a solution rather than getting a very complex solution or no solution. So with that spirit, we will see that how to get an approximate solution to this problem. How does the approximate solution compare with the exact solution? We will work out the exact solution for this problem also and then I will give you a homework that you plot the approximate solution, you plot the exact solution and compare the error between the 2 solutions for different ratios of A by B. Possibly for the parallel plate channel case, you will find that that error is not very high but as the aspect ratio will become close to unity that is A by A is to B will be close to 1 is to 1 that error may be more and more significant. We will see why? This is not, this is x momentum equation, left hand side is, left hand side is the acceleration term, it is fully developed flow, del u del t plus u del u del x plus v del u del y plus w del u del z. What is del u del z? What is w del u del z? What is w? The flow is taking place in the yz plane, right. So w, so this is basically, so you have the flow is in the yz plane means what? The variation of the velocity profile is because of variation in z and y. So if you consider that you will see, I mean you have only one velocity profile, you have only one velocity component that is u and the other velocity components will be trivially 0 by applying boundary conditions by using the continuity equation, we have not worked out that details. So v and w are 0, okay. So in conjunction with that if you apply the fully developed flow consideration the left hand side will be 0, okay. So basically the flow variation is in the yz plane, the axial component of velocity varies in the transverse plane, that is what is the physical scenario, that is the physical concept that we can borrow from the, see why we start with the parallel plate channel? We could have started with this only, the parallel plate channel case gives us a lot of physical insight and we are just borrowing that physical insight directly in this example. Now you will see that if you write y momentum and z momentum equation you will have del p del y equal to 0 and del p del z equal to 0 because those velocity components themselves are 0. So that means p is not a function of y and z, so p is a function of x only, that is what I am saying that it is important, see if you want me to spend more time on deriving this equation by starting from the scratch I can do that. But at this level of maturity I expect that once I have covered the parallel plate channel you should be able to follow how this equation comes, it is just an additional dimensionality that has been added but physically the concept is very similar, okay. So again to repeat where from it follows that del p del x equal to dp dx if you refer to the y momentum and z momentum equations those velocity components are 0 that will give rise to del p del y equal to 0 and del p del z equal to 0 that means p is not a function of y and z. So p is a function of x only that means the del p del x will become dp dx. Now let us come back to the approximate solution, see let us say that I am a person who is not so conversant in advanced mathematics, still I want to solve this problem and want to design an engineering system out of that, so how can I get rid of that, see one of the clues that I can get is that possibly the velocity profile had it been a parallel plate channel would have varied quadratically with y that is what is the solution that we have seen for a parallel plate channel, had it been a parallel plate channel with variation along z that would have also quadratically varied with z. Now there is a combined variation with respect to y and z, so you really do not expect that from a fundamental mathematical perspective it will vary simultaneously quadratically with y and z but as an approximate solution we can consider that that is not from a very rigorous theoretical basis but from an approximate basis we can take a solution of the form assume u is equal to c into y square minus s square by 4 into z square minus b square by 4, so what is the motivation of taking this form, the motivation of taking this form is the quadratic variation of velocity with respect to y in the parallel plate case, had it been a parallel plate channel with variation along z it would have been this. Now this form we are taking quadratic but why this form, this form we are taking to satisfy the no slip boundary condition at y equal to plus minus a by 2 and z is equal to plus minus b by 2, see always keep one thing in mind many times you do CFD solution many times you do approximate analytical solution, many times you do detailed rigorous analytical solution whatever you do the first and foremost thing that we always check as a check to the solution that the solution satisfies your boundary conditions, even wherever we do a CFD problem we see that whatever contours or vectors whatever we are getting those are satisfying the boundary conditions that is the first check that we always make, so here the motivation of taking this form is 2 fold one is to keep a form that we get for the parallel plate channel, so because this form we get for the parallel plate channel possibly if this channel was such that it is a parallel plate channel that is either of a by b or b by a very large or very small then it would be quite accurate, the second motivation is that this will automatically satisfy the boundary conditions at the 4 boundaries, so once we do that this is heuristic remember this is not correct but this is something which is approximate and we will see that under certain limiting conditions it is not that bad, so the velocity profile may be specified completely if you know what is the value of c, so the question is how do we find out the value of c, no, no, no this is not imaginary wall this is the channel is like this, so these are walls these are these are all walls this is not I mean these are not imaginary walls this is a rectangular channel, problem is these are certain sort of mental blocks that we always have like from the second year undergraduate level we have understood that channel means parallel plate channel, so when we are putting a boundary I would say like if you tell to a person who has not learnt fluid mechanics this is something which is more obvious than a parallel plate channel, why should I suddenly assume a parallel plate channel I should better assume that it is bounded by the walls rather than infinite in one direction, so just do not keep any mental block, assume I mean you have all the boundaries not infinitely large in one direction, so this is what is u now to get c what we will do is we will use the integral form of this equation, so what we will do is we will substitute mu okay, so first let us work out the derivatives del u del y is equal to 2cy into z square-b square by 4 this is del 2u del y2, why we are calculating this we will require that to be substituted in the momentum equation, so you will see a very interesting thing let us complete this del 2u del z2 what will be that 2c into y square-s square by 4, if you add these 2 right, so that will become 2c into z square-y square-s square by 4-b square by 4 right hand side will be equal to dp dx and there is a mu adjusting factor, now there is no guarantee that these 2 will be equal right, so we can clearly see that this equation although or this expression although satisfies the boundary condition does not satisfy the governing equation right, so it is not a correct solution, it may be an approximate solution but it is not a correct solution but you may accept a level of approximation by accepting that it will not satisfy the governing equation in its differential form but maybe it will satisfy the governing equation when it is integrated over the cross section, so that is what we will do now, so what we will do is we will just write it a little bit like this 0 is equal to this as 1 by mu dp dx just a little bit readjustment, so this will be 2c into z square sorry y square-z square-s square by 4-b square by 4 is equal to 1 by mu dp dx. Let us integrate this over the cross section of the channel, so dy dz y is equal to what? y is equal to-a by 2 to a by 2 z is equal to-b by 2 to b by 2 okay, this also we integrate dy dz, we can safely keep dp dx out of the integral because p or dp dx varies with x only, it does not vary with y and z, I am not going to enter into the trivial algebra and waste time in this, you can yourself complete this to get the value of c, this is just a simple polynomial integration and at this level it is very straightforward, so once you calculate this it will be c into some number a right, the number a depends on the integral of this y square-z square all these things is equal to 1 by mu dp dx into ab, so c is equal to 1 by mu ab by adp dx, so that completes the description of the approximate velocity profile. Once you get the approximate velocity profile, you can find out the wall shear stress, how do you find out the wall shear stress, so wall shear stress depends on question now is which wall right, so you can calculate the wall shear stress here or I mean here and here it will be equivalent or here and here, when we are talking about parallel blade channel, when we say wall we know that right just one either of the parallel walls, now there are 4 walls, so either of this or of this and then you can non-dimensionalize that to get a friction coefficient, I am not working out that because that is just trivial algebra, now that solution if you work out, you will see that it will satisfy the boundary condition without satisfying the equation, so there is a level of approximation associated with it, question is what is the level of approximation, how accurate it is, is it good for some cases, is it okay for some cases or is it poor for all cases, to assess that we have to go for the exact analytical solution for this problem which we will attempt now, so this is the equation that we are solving, when you first look into a partial differential equation which is the linear partial differential equation, the non-linear terms have been dropped because of the fully developed flow condition, what is the technique which we are always tempted to use, that is method of separation of variables, what prohibits us to apply the method of separation of variables here directly, that is the inhomogeneity of the equation, the right hand side is not 0, right, so had we have the case when the right hand side is 0, it would have been easy and straight forward for us to apply the method of separation of variables, to convert this equation to that form, what we do is that we separate this equation into a collection of 2 problems, what are the 2 problems, problem 1 is a parallel plate channel problem, what is the parallel plate channel problem and maybe let us consider because of symmetry only one fourth of the domain, u1 is a solution which is a parallel plate channel solution, what are the boundary conditions at y is equal to 0, du1 dy equal to 0, at y is equal to a by 2 u1 equal to 0, okay, this is problem 1, problem 1 is the parallel plate solution, see this is one important philosophical way of looking into education, we always try to build up from what we know to solve something what we do not know, so here we see that we have already worked out a problem of parallel plate channel, we are trying to build up a solution of this more complicated problem by starting with that as a basis, so that is problem 1, but there is a problem 2 which is this, problem 1 already takes care of the pressure gradient, so problem 2 does not have to bother about the pressure gradient, so this makes the governing equation at least homogeneous and this is a very standard partial differential equation in mathematics or mathematical physics, all of you know that this is called as the Laplace equation, okay, so Laplace equation in 2 dimensions, what are the boundary conditions at y is equal to 0 del u2 del y equal to 0, right, let me use this side because otherwise the board will not be visible, at y equal to 0 del u2 del y is equal to 0, then at y is equal to a by 2 u2 is equal to 0, right, y at y equal to a by 2 u2 equal to 0, see at y equal to a by 2 u has to be 0 and we want to use u as a superposition of u1 u2, this superposition we can use because of the linearity of the governing differential equation, because the governing differential equation is linear, when u equal to u1 is a solution and u equal to u2 is a solution, u equal to u1 plus u2 is also a solution, so that we can use, so u1 plus u2 is the solution which is 0 because we have already taken u1 equal to 0, u2 also has to be taken as 0, so that u1 plus u2 equal to 0, at z is equal to 0, what is the boundary condition, del u2 del z equal to 0, this is just because of symmetry, centre line symmetry, at z is equal to b by 2, what is u2, what is this, u2 is not 0, the reason is that at z is equal to b by 2, here we have no control over u1, so u1 is not 0 and u1 is a function of y, u1 as a function of y is solved from this problem, so we can say that u2 is equal to minus u1, right, u2 is equal to minus u1 which is a function of y, how are these problems coupled, this u1 as a function of y you will get from problem 1, okay and u is equal to u1 plus u2 is your solution, problem 1 we have already worked out, so I will not waste any time by working it out again, we will see how to work out the problem 2 and then the superposition you can get the solution, we see we can start with the solution in the form of method of separation of the variables because the governing equation is homogeneous, not only that 3 out of the 4 boundary conditions they are also homogeneous, does not matter one of the remaining one is inhomogeneous, still it is okay, so separation of variables u2 we can take as fy into gz, so del 2 u2 del y2 is equal to f double dash del 2 u2 del z2 is equal to f double dash del 2 u2 del z2 is equal to f double dash f g double dash okay, so this is the governing equation that needs to be satisfied, so if you want to satisfy this governing equation f double dash g plus f g double dash is equal to 0, so you have f double dash by f is equal to minus g double dash by g, this is a function of y only, this is a function of z only that means each is equal to a constant, question is and a very important question whether it is a positive constant or a negative constant okay, so we have to figure that out whether it is a positive constant or a negative constant, now here we can show and this I leave on you as an exercise that this is indeed a negative constant, you can arrive at this conclusion in 2 ways, one way is that you use a method of contradiction, you say that no it is a positive constant plus lambda square, you will see that that form will not be able to satisfy all the boundary conditions, otherwise from an intuitive view point if you see if you take this as minus lambda square then g will be coming in the form of sin and cos right or not, f will come in the form of sin and cos and g will be coming in the form of exponential and if you take it the other way then here you can see that f is coming in the form of exponential and g in the form of sin and cos, so in one case it is f in the form of exponential and g in the form of sin and cos, in one case g in the form of exponential and f in the form of sin and cos, now you have to figure out that out of these 2 which one is physically acceptable by considering the boundary conditions that u2 is 0 at both these boundaries but u1 is not 0 at both these boundaries or u2 is not 0 at sorry u2 is 0 at these 2 boundaries but not 0 at these 2 boundaries I am sorry I just may try to mention the other boundaries, so u2 is 0 at these 2 boundaries but it u2 is not 0 at both of these boundaries, so that gives the source of thought which will dictate that which one will be exponential and which will be sin and cos but even if you are not very confident with that and you start with the plus lambda square it will make no difference except that you will come up with a solution that will not satisfy the boundary condition. We will take up this issue in the next class for the time being thank you very much.