 In the previous lecture, we calculated friction factor for fully developed laminar flow in somewhat regular ducts like a circular tube or an annulus or a rectangular duct or an annulus sector. We used simple algebraic methods as well as the Fourier method. In today's lecture, I am going to consider ducts of even more complex cross section, such as a triangular duct and ultimately a duct of any arbitrary cross section. So, let us begin with the duct of triangular cross section and apply one method called the Cantorovitch variational method. The figure shows here a triangular duct of base 2 b and height 2 a with an included angle apex angle 2 phi governing equation would be d 2 u dx square plus d 2 u dy square equal to 1 over mu dp dz, which I can non dimensionalize in this manner u star equal to u divided by 4 a square 1 over mu dp dz and x star equal to x divided by 2 a and y star equal to y over 2 a. x is measured from the apex and so is the y. So, the origin is the apex of the triangle. Of course, this is a Poisson's equation with the right hand side equal to constant, but notice that the boundary conditions are given at x star equal to 0, which is the apex itself and x star equal to 1, which is the base and u star is also equal to 0, where y is equal to plus minus f x and that f is a function of x itself. So, we have a boundary now, which does not have a constant y, but a y which is a function of x and in this case a linear one with m equal to tan phi this angle is 2 phi. So, phi is the half angle for such an equation solution can be obtained by variational method due to Kantorovic. Thus, let u star which is a function of x and y be given as f square minus y square multiplied by some function of x star. The objective is to find out what is f x star by satisfying the boundary conditions and the equation. We will restrict our attention to 2 phi less than 90, so that m is always less than 1. Remember m is equal to tan phi. The variational principle is stated like this. The variational i is 0 to 1 minus f to plus f, the equation itself multiplied by variation in u star dx star dy star is equal to 0. Since, dx by dx star is equal to m equal to tan phi is equal to constant and in this case d 2 f dx star square will be 0 and then you will see that the equation itself u star is equal to f star f square minus y square f. Therefore, du star by dx star will be 2 f d f by dx into f plus f square minus y square into d f by dx star and remember f is equal to m times x star. Therefore, this simply becomes 2 f m into m square x star into f plus f square minus y square d f by dx star and d 2 u dx star square will become 2 m square f plus 2 m square x star d f by dx star plus f square minus y square d 2 f star by dx star square plus d f by dx star into 2 f d f by dx star and that we saw is equal to 2 m, this is equal to 2 m square x star and therefore, you will see that this equation d 2 u dx star square would be written as this is f star square minus y square f double prime plus this is 2 m square x star which is essentially 2 m f and so is this 2 m f. So, this becomes 4 m f d f by dx star which is f dash plus 2 m square f likewise d 2 u star by d y star square will simply be equal to minus 2 times f. So, if I substitute for this and this in this equation minus 1 delta u star and if I further then carry out the integration first with respect to y star then I will get 4 by 3 f delta variation of 0 to 1 4 by 5 f square f double prime plus 4 m f f prime plus 2 m square minus 1 whole square into f minus 1 into f dx star is equal to 0. This shows that the term in the square bracket itself must be 0 or 4 by 5 f square f double prime plus this term minus 1 equal to 0. If I now define f star is equal to f minus 0.5 divided by m star minus 1 then since f is equal to m x star you will see that I get x star square plus f double prime plus 5 x star f prime into this plus line. So, I get a second order equation in f star second order equation in f star. This equation can be cast in this form 1 over x star cube remember this equation can be cast in this form 1 over x star cube d by dx star x star phi over d phi by dx star into this phi by 2 function of m and f star equal to 0. The solution is f star is equal to f minus 0.5 which we have said before is equal to a x star raise to r 1 plus b x star square raise to r 2 or f itself is just 0.5 m square minus 1 raise to minus 1 a x star raise to r 1 plus a star raise to r 2 and u star is then given by this expression. Remember u star is related to f in our equation in this manner and therefore we have obtained f and therefore we have obtained y. So, where r 1 will turn out to be this and r 2 will turn out to be this constants a and b are now to be determined from the boundary condition u star equal to 0 at x equal to 0 and x equal to 1 that is the apex and the base of the triangle. Condition at the base x star equal to 1 gives a plus b equal to that I put x star equal to 1 then I get this a plus b equal to 0. For m less than less than 1 because we are restricting to 2 phi less than 90. So, 10 phi or m would always be less than 1 r 2 will be less than 0 and therefore at x equal to 0 this would tend to infinity which is of course unacceptable and therefore b itself must be 0 as a result u star turns out to be a function with x star raise to r 1 minus 1. Integration of this gives us the u bar star and that evaluates to this quantity r 1 plus r 1 further the hydraulic diameter for this particular triangle can be shown to be equal to 2 m into m plus square root of m square plus 1 raise to minus 1 and hence friction factor must be multiplied by Reynolds number is a function of m that is the tangent of the half angle and a function of r 1 as you remember is again a function of m. So, in this particular problem after considerable mathematical manipulations we have shown that the friction factor would be function of the included angle. Here are some solutions. So, I begin with 85, 75, 60, 50, 40, 30, 20, 10 and 5 the value of m that is the tan phi is given the evaluated value of r 1 is also given the hydraulic diameter values are also given and here are the friction factor versus Reynolds numbers that are evaluated for each geometry. A special case is 2 phi equal to 60 remember 2 phi equal to 60 would straight away give us an equilateral triangle 2 phi equal to 60 gives us an equilateral triangle for which the friction factor is 13.33 it is an often quoted value what is less often quoted are the values which I have shown for other angles methods of this type are not general for different ducts such as elliptical or triangular with rounded corners for example, which is very common lay encountered in practical heat exchangers different strategies that is the trial function whose variation is to be taken has to be invoke and then considerable algebra is required before one gets a solution and for each duct one has to treat it as a special case and go through the solution developing the solution. Our interest now turns to more complex ducts in which case such methods do not apply and therefore, we seek now a general method which can be applied to ducts of arbitrary cross section absolutely arbitrary cross section of which a triangle or an ellipse or any other even a circular duct would be a special case we begin with that method. So, consider for example, a duct of arbitrary cross section and the coordinates of the boundaries z b y b are known the flow is in the x direction that is into the plane of the screen and only requirement is that the domain must be singly connected. Singly connected domain implies that once I put a pencil at any point on the domain and I go anticlockwise or clockwise does not matter I must return to the point of origin without lifting the pen I should be able to return to the point where I started from without lifting the pen. So, that is the only constraint on this method in other words if I had for example, is a circular rod inside a solid rod inside that would not apply for this method. It has to be a singly connected domain whose boundary coordinates are known z b y b the governing equation would be d 2 u d z square plus d 2 u d y square equal to 1 over mu d p d x equal to a constant. But if I define u divided by minus 1 over mu d p d z equal to u star minus z square plus y square by 4 then substitution of this in this equation would readily show that a Laplace equation results in other words whose right hand side is 0. And since u is 0 on the boundaries u star would be equal to or u star boundary would be equal to z b square plus y b square by 4 u star b of course, is a fictitious function which takes finite values on the boundaries. How do we solve this Laplace equation on a domain which is completely arbitrary? The solution for such a Laplace equation is given as u star z y is equal to sum of coefficient c i multiplied by some functions g i z y. So, c i's are the coefficients to be determined and g i's are prescribed by exploiting a very special property of the Laplace equation and that property is the following. For any positive integer n the real and imaginary parts of a complex variable z plus i y raise to n are each exact solutions of the Laplace's equation. So, if I assigns successively values n equal to 0, 1, 2, 8 etcetera I would get the first 17 solutions which usually suffice for most complex geometries. You are of course, welcome to take even 10 or 12 and generate 21 or 25 solutions. So, the functions g i would read as follows and I will show you how they are evaluated z plus i y raise to n is each a solution. So, if I take n equal to 0 I would simply get g 1 equal to simply 1 and that is for n equal to 0. As you can see here, there is no imaginary part nor a real part and therefore, g 1 would simply equal to 1 for any value of n because n is equal to 0. If I take n equal to 1 then the real part will give me z and the imaginary part will give me y. As you can see here on the screen g 2 is equal to z and g 3 is equal to y. So, these two solutions are for n equal to 1. If I take n equal to 2 then you can see I will get g 4 equal to z square minus y square and g 5 will be equal to simply 2 z y because this is the imaginary part and this is the real part because z plus i y square is simply z square minus y square plus i times 2 z y. So, we take both these as solutions and this is for n equal to 2. Likewise, you can go on taking so for example, g 6 and g 7 here are solutions for n equal to 3, g 8 and g 9 are solutions for n equal to 4, g 10 and g 11 are solutions for n equal to 5, g 12 and g 13 are solutions for n equal to 6, g 14 and g 15 are solutions for n equal to 7, g 16 and g 17 are solutions for n equal to 8 and as I said earlier you could take n equal to 9, n equal to 10 as much as you like and in each case you will get two extra terms in the equations. So we do have g i z y's are known and the only thing unknown is c i. So, suppose I choose in a particular duct let us say 16 boundary points then I only take first 16 terms in that expression and my task then is to determine c i equal to 1, 2, 3 up to 16. So, 16 coefficients are to be determined from 16 boundary conditions u star z b y b which is nothing but z b square y b square by 4 which is known because we know the coordinates of the boundary points equal to c i g i z b y b the functions are evaluated at the boundaries themselves. So, essentially you get 16 equations and 16 unknown c i which are to be evaluated one way in which this can be done is by the method of l u decomposition. This is a metric solution method followed by forward elimination and backwards substitution procedure which you must have studied in your numerical analysis course. Other method for this type of equation sets is what is called Gram-Schmidt orthonormalization. I will be following l u decomposition procedure all the results that I will show are with l u decomposition method. I will now apply this method to a variety of ducts. So, let us take first of all elliptical duct as shown the origin 0 0 its major axis is A and its minor axis is B and in this particular case I have taken A equal to 2 and B equal to 1 and as you can see I have chosen 17 boundary points they need not be equally spaced although I have shown I have taken them to be equally spaced but they need not be. So, 1 2 3 4 5 6 7 8 9 10 11 13 14 16. So, I have essentially taken 16 points. So, I know the coordinates of 16 points because the equation of an ellipse is simply z square by A square plus y square by B square equal to 1 that is the equation of an ellipse. So, I know exactly if I choose the values of z then I know the values of y B for each of these points. Here are the evaluated constants C i. So, you can see I have taken 16 boundary points 2 0 1.5 z B equal to 1.5 1.5 0 minus 1.5 and in each case y B has been evaluated the right hand side which is z B square plus y B square by 4 is also mentioned here and after evaluation by L U decomposition method what do I find? I find that only constants C 1 and C 4 are finite all other constants are identically 0. In fact, no matter what value of A and B you will take you will always find this for a special case of an ellipse that only C 1 and C 4 are finite the rest of the coefficients turn out to be 0 which makes for a very elegant solution. For example, the final solution then is simply remember U is equal to U star which is C 1 plus C 2 C 1 plus C 4 into the function G 4. Do you see that G 4 was z square minus y square and therefore, the solution is 0.4 plus 0.15 z square minus z square plus y square by 4 U equal to U star plus z square minus plus y square by 4 and we found U star from the analysis. So, U itself becomes 0.4 minus 0.1 z square plus 0.4 y square in this case integration is done from 0 to A 0 to y B U dz dy divided by 0 to y B dz dy and y B from the equation of an ellipse is simply 1 minus z B square by a square raise to half. So, you carry out these two integrations you will see that will turn out to be 0.2 and as a result in this particular case for A B equal to 1 and A equal to 2 the solution is U over U bar plus 2 minus 0.5 z square minus 2 y square. If I were to generalize the solution is like this as I said although C 1 and C 4 are functions of A and B they are the only ones that remain finite all others are 0. So, the general solution is this U bar would be then given by this to calculate friction factor into Reynolds number product which is d h square by 2 U bar over 1 over mu d p d x which we have shown in our previous lecture that and d h square for an elliptical duct is given as 4 times area divided by vetted perimeter area of course is simply A B into 5, but the perimeter is given by this polynomial expansion and lambda here is A minus B divided by A plus B. If you wish to see this you can look up any mathematical table in which elliptic integrals evaluate perimeter equal to that. So, here are the results for variety of values of B and A B equal to 1 in all cases, but A is chosen to be 1 1.25 1.67 2 2.55 and 10 which is a very oblong ellipse and you will see in each case C 1 and C 4 have been evaluated also given are the values of d h by B. The special case is when A is equal to 1 and B is equal to 1 that is a case of a circular tube and as you would expect in this case C 4 is 0 only C 1 is finite and the f r e is equal to 16 as we expect, but for all other cases f r e goes on increasing the more oblong the ellipse the higher is the friction factor. These are very useful solutions particularly because many heat exchangers use elliptical tubes for heat transfer. Let us look at another problem that is a triangular duct with rounded corners rather than sharp corners. Manufacturing wise rounded corners are often more conducive to manufacture and we wish to account for the effect that rounding would have. The rounding radius in this case is B and the unrounded side is A. So, it is really an equilateral unrounded triangle in which rounding is given at the three corners. So, B is the rounding radius A is the unrounded side the origin is here and y and this is x and I choose let us see how many points I have chosen in this case these are the points chosen. As you can see I have concentrated points here where the curvature is taking place and only one point on the whole entire side. Now of course, here I am not giving you the value of Z B and y B they are very obvious, but notice this that in this case all coefficients 1 2 3 4 5 for 17 points which I have chosen all coefficients are finite none of them is 0. Although some are small like this one for example is minus 0 371 whereas, this is minus 58 which is very large. So, they are very uneven magnitudes, but nonetheless some of the coefficients are very small and some coefficients are very large. So, all the 17 terms in our expansion would have to be retained to express the solution. To evaluate friction factor and Reynolds number then you must evaluate the u bar from this solution which has to be done on a computer because hand integration becomes very difficult. So, you have to do that on a computer and here are the results. So, in each case I have taken A equal to 1 and the rounding radius first of all is 0 which of course, corresponds to an equilateral triangle with sharp corners 0.05, 0.1 and 0.167 and you can see I have evaluated cross sectional area divided by A square perimeter divided by A and therefore, hydraulic diameter divided by A. Now, you can see that for a sharp corner it is 13.33 which we had also seen in the Kantorovich method, but as the rounding radius increases you see that the friction factor versus Reynolds number seems to increase, but based on hydraulic diameter where the hydraulic diameter itself goes on increasing with the rounding radius. Let us take yet another problem, this is a circular segmented cross section. So, the duct itself has a flat side and a part of a circle part of a circle and the apex angle which is that apex angle is theta. So, the duct boundary is a flat side and a curved side. Now, you can imagine such a geometry would be extremely difficult to handle by any other method because it is so irregular. Again, I have chosen 16 points 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 and 16 and solutions are obtained for a prescribed value of the apex angle theta and radius B. Here are the results, theta equal to 90, 60, 45, 30 and 10. In each case the hydraulic diameter to B the radius is given. You can see now in this case C 3, C 4, C 7, C 11 and C 15 are the only finite values. All others are almost 0. When you actually print out from your computer after l u decomposition, you may find these values to be of the 10 to minus 8 or 10 to minus 9 which are not worth considering at all because they are so very small. So, only 1, 2, 3, 4 and 5 coefficients are finite in this particular case. Again, using these 5 coefficients you evaluate the average velocity u bar which enables you to calculate friction factor multiplied by Reynolds number product. Theta equal to 90 is a very special case because it represents a semicircular duct of a semicircular cross section. As you can see from this figure, if theta were 90 then I would have a perfect semicircular duct. So, then the value is 15.76 but then the value goes on decreasing but not so much. Nonetheless, these accurate values are very, very useful when we compare with experimental data. We want to be sure that calculations and predictions compare very well. So, what we have got now is a very, very general method. You can in fact program the entire method for a variety of ducts. All you have to do is to give boundary point coordinates for different types of ducts and that is what I have done. I have got a single computer program which carries out which has all the 17 functions stored in it and it has a composition subroutine as well as it has the subroutine to calculate u bar which requires numerical integration. So, once that is done all that is required is for each different type of a duct. All you need to do is to give coordinates of the duct boundaries. The only constraint of course is that the duct that you choose must be singly connected. It must have singly connected cross section. So, as you can imagine the method does not change with the duct shape. It simply requires bound whereas, unlike Kantorovich method or the Fourier series methods and so on and so forth where in each case you have to develop different functions to satisfy the boundary conditions and so on and so forth. Whereas, here you do not have to worry about that. You simply give the boundary points, coordinates and that is it. The method is general and can be applied. Moreover, the method can also be extended to heat transfer and that we shall see in couple of lectures from now. So, lecture number 18 you will be able to see that. The method has been developed in the paper by Sparrow and Hage-Shake in Trans-Hasmay journal of heat transfer 1966. Several non-circular ducts of variety of cross sections that I mentioned in my earlier lecture moon shaped duct, sign duct and so on and so forth results for them are given in this companion by R. K. Shah and London. Laminar force convection in ducts advances in heat transfer volume 15 1978 and might I say that the results which I presented here have been corroborated in this companion.