 In the last lecture, we converted the partial differential equation of a two-dimensional boundary layer to an ordinary differential equation, which we call the similarity equation. The purpose of this lecture is to show how to solve such a similarity equation by shooting method. In this lecture, I will also present some solutions to the velocity boundary layer equation. Let us see our interest is to solve this third order ordinary differential equation with the boundary condition f 0 equal to b f, which we defined as the blowing parameter into 2 m plus 1. f dash 0 is equal to 0, which is the no slip condition and f dash infinity equal to 1 is the boundary condition. So, our solution is for the velocity profile f dash eta u over u infinity will be function firstly of the pressure gradient parameter, because u infinity is a function of m and b f which arises from the manner in which v w the wall velocity varies. Secondly, our velocity distribution in y direction would be given by that. So, once we get this once we solve this equation, we get f f dash variations from which we get the velocity distribution v as well as u as a function of y. Just by way of reminder, I may say eta, the similarity variable eta is nothing but eta is y times under root u infinity divided by nu x, where u infinity is c x raise to m. Therefore, you will see this becomes nothing but c by nu x raise to m minus 1 by 2. That is what the variable eta is and this is the stretching parameter. Similarity solutions are possible only when u infinity varies in this fashion. We say m is a pressure gradient parameter, why? Because, remember d p infinity by d x is simply minus rho u infinity d u infinity by d x and that is simply is equal to minus rho c x raise to m, c into m into x raise to m minus 1 or that is equal to minus rho c square m into x raise to 2 m minus 1. So, you will see the pressure gradient is totally determined by whether m is positive or negative and therefore, we call this m as the pressure gradient parameter. This is called the similarity independent variable and f dash eta is equal to u or u infinity. So, now, we want to solve this third order equation, third order ordinary differential equation. We essentially split it up into three first order differential equations as shown. Before I do that, let me just tell you the f dash eta solution gives the coefficient of friction as a function of the Reynolds number. So, the shear stress at the wall would be mu times d u by d y is equal to 0. So, tau wall will be mu times d u by d y at y equal to 0, which I can also write mu times u infinity d u by u infinity by d y at y equal to 0, which I can also write as mu u infinity d f by d y d f dash by d y at y equal to 0, which I shall write as mu u infinity d f dash by d eta into d eta by d y at eta equal to 0. As you know, since eta is equal to y into u infinity nu x, you will see this different becomes mu times u infinity under root u infinity by nu x, which is d eta by d y into f double prime 0 and we define c f x as the shear stress divided by the kinetic energy of the free stream. Therefore, this will become 2 times f double prime 0 into 2 Reynolds x to the power of minus 0.5 and that is what I shown here, the local skin friction coefficient would be defined as that. So, once we know, we have solved the problem, we can determine f double prime 0 and sometimes it is also of interest to find out average skin friction coefficient over length l. Well, all you do is integrate tau all x from 0 to l and divide by the length l and you will get that. Therefore, we must determine f double prime 0 for which we have no boundary condition at the moment. Further parameter of interest will be listed in a later slide. So, as I said, we split up the third order equation into three first order equations. The first is the d f by d eta f dash and I can solve it from eta equal to 0 to eta equal to infinity because I know the boundary condition f 0 b f 2 over m plus 1. This is known because I have already specified the value of b f and I have already specified the value of m. d f dash by d eta will be equal to f double prime. This is the second ordinary differential equation with f dash 0 equal to 0 because this is a known condition. Similarly, f double prime by d eta will be f triple prime and f triple prime will be simply from our ordinary differential equation. It will be this quantity with a negative sign and that is what I have written here. But to solve this equation, I do not know the initial condition at f double prime at eta equal to 0 and that is an unknown condition. So, how do I do that? I will show that later. But since these three equations are ordinary differential equations, I can always solve them by Ringe-Kutta method with which you are already familiar. From eta equal to 0 to eta equal to eta max, I say eta max in lieu of infinity because it is no point going on solving for a length longer than is necessary. Usually, values of eta max of 3 to 10 will suffice depending on the value of b f and m. So, eta max is a fictitious quantity. It represents infinity condition because in a numerical calculation, you must give some value to infinity and that is given as eta max. What would be the solution algorithm? Select the value of m and b f which is the given. You want to solve this problem for a given pressure gradient and given blowing parameter. Select eta max and step change d eta. You are solving from eta equal to 0 to eta equal to eta max. The conditions known here are f 0 f dash 0, but f double prime 0 is not known. So, this eta equal to 0 to eta max is divided into very small number of steps. We call them d eta. Each step is called d eta and we intend to obtain solution at different values of eta along up to n. So, solve three equations simultaneously by Runge-Kutta method by first guessing a value of f double prime 0. Now, if I choose a value which is not correct, then what do I expect? I expect f dash eta max to somehow go to 1 at correct value of f dash is 1, but what I will find? If my guess is not quite correct, I will find I end up like so. I will end up like this. Obviously, I must refine my guess for f double prime 0 and in the next guess I may find that I come out in that position. So, this is guess 1, this is guess 2. Neither is f equal f dash eta max is not equal to 1. So, I check if f dash eta max is equal to 1. Since it is not, I must refine these two first two guesses and I go on doing that by this little formula. It is simply a linear interpolation of the errors on the both sides so that the third guess will give me an error in f double prime 0, which is more and more accurate. If I say f double prime 0 as a phi, then phi at k plus first iteration it would be equal to phi k plus 1 minus psi k phi k minus phi k minus 1 psi k over psi k minus 1 and where k is the iteration number and psi is the value of f dash eta max predicted at eta max that is the value psi. So, I sense which way to move for f dash double prime 0 by observing the value of f dash eta max eta equal to eta max. If of course, I now find that f dash eta max is equal to 1, then I say assume value of f double prime 0 is a correct one. At convergence then print values of f eta f dash eta and f double prime eta and note down the value of f double prime 0, which is required for calculation of the scheme friction coefficient. Here I show a typical calculation for a set of assuming m is equal to 0 and b f is equal to 0. So, the equation that I am really solving is simply f triple prime since m is equal to 0, m f f double prime equal to 0. This is the equation I am solving for m equal to 0 and the boundary condition is that f dash 0 is equal to 0, f dash infinity is equal to 1 and f 0 which is equal to b f into 2 by m plus 1. Since this is 0 and b f is equal to 0, f 0 is also equal to 0. So, I have taken in this case after some experience that eta max should be about 7. So, I have to eta equal to 7 because I am solving for m equal to 0 and b f equal to 0 which is essentially a flat plate boundary layer without suction or blowing. I have divided this into 300 steps. I am not suggesting that 300 is really required. You can get away with smaller number, but I took a value of 300. So, I have d eta is 1 divided by 300 of eta max and since I do not know the solution, I made an initial guess of f double prime equal to 0.02. So, this was my first guess and I found eta max f dash eta max equal to 0.123. Obviously, this is not correct. I then refined the guess to 0.07 and I found that it is 0.342. With these two solutions, I let our iteration refinement formula to take over and then I find that the third guess was 0.22 and it became 0.761, fourth guess 0.306, 948, fifth 0.329, 0.997 and sixth gave me 0.33071 and f dash max is indeed 1 which means that this value of f double prime 0 is a correct one. Now, when you are solving this for the first time, obviously you do not know f double prime 0. You do not know what eta max to take. So, the best thing is to increase the value of eta max, decrease the value of eta max and in each case determine the value of f double prime 0 corresponding to f dash max equal to 1. So, that our solution should be independent of the number of points taken in the domain and the value of eta max chosen. Of course, you do not want eta max to be too large. At the same time, it should not be so small that you do not resolve the entire viscosity affected region properly. Because I started with a very poor guess of 0.02, I needed six iterations to discover the correct value of f double prime. In this case, therefore, Cfx which is 2 times f double prime 0 rex to the half will be simply 2 times 0.33. So, that is equal to 0.6614. Cf bar will be 1.28. It is also possible that the equation that I wrote here, there is a series solution possible for this equation and that series solution gives Cfx equal to 0.664 and Cf bar equal to 1. So, we are very close. Our shooting method has given very close answer to the series solution. Now, you can see how the profiles look like for the flat plate boundary layer without suction and blowing. You will see the value of f dash goes from 0 to 1 very asymptotically and that is u of 0.664. And where f dash becomes constant, f becomes linear and the value of f double prime, you can see 0.33 and when you go towards the n, you will see that the second derivative of f or the du dy is also going to 0, which is what we expect. So, that also is a fair correct. Actually, u over u infinity becomes very close to 1 at about eta equal to 5. So, we took eta max equal to 7. Well, that is good enough. Even if I took eta max equal to 6, I would get exactly the same results. Normally, one does not want this to be too large. The difference between the edge of the boundary layer and the eta max chosen should not be too large because unnecessarily you are doing computations which need not be done. Now, in order to interpret the solutions in boundary layer theory, we define certain thicknesses. One of the problems with boundary layer theory is that remember u over u infinity goes very asymptotically to 1 and one cannot really say what is exactly the thickness of the boundary layer. It is a notional quantity. Boundary layer thickness is a notional quantity and we associate that with the notion that we shall say boundary layer thickness delta is the value of y, where u over u infinity f dash eta is equal to 0.99. That is sort of a convention that is followed. But a more exact representation would be displacement thickness delta 1 and we define that as 0 to infinity 1 over rho u divided by rho infinity. Now, you can see what does that represent? It simply represents the idea that when a boundary layer is formed because of this viscosity affected region, if the boundary layer did not exist, everywhere the fluid will be flowing at the velocity u infinity. But because of the boundary layer, less fluid is flowing through this area and therefore, we define delta 1 as 0 to infinity rho infinity u infinity minus rho u divided by rho infinity u infinity dy. You can see therefore, this quantity represents the amount by which the mass deficit has occurred because of the low velocity region in the vicinity of the wall. In our case of course, since we are saying uniform property flow rho is equal to rho infinity and therefore, we will simply have 0 to infinity 1 minus u over u infinity dy. If I say this is equal to 0 to delta plus delta to infinity 1 minus u over u infinity dy, then you will notice that this is really beyond delta, the boundary layer thickness, u is in fact equal to u infinity and therefore, that is equal to 0. So, I get this is really over u infinity dy. This thickness is called the displacement thickness. So, even if I choose eta max which is much greater than delta in dimensionless form, it does not matter because all the integrations after that are 0 and therefore, displacement thickness is a far more precise quantity than the quantity delta which is a sort of a fictitious value at f dash eta equal to 0.90 and therefore, these thicknesses are much more reliable indicators of the thickness of a boundary layer. Likewise, momentum thickness is defined as mass deficit multiplied by u. It simply tells you the amount of momentum deficit that has been caused by the presence of the boundary layer. Again, if rho is equal to rho infinity, then this simply becomes u over u infinity 1 minus u over u infinity and since u over u infinity is f dash and dy is related to eta, I can convert these definitions to this form delta star is equal to delta over x re x half delta 1 star delta 1 x. This is fairly straight forward to show. For example, I can define delta 1 as 0 to delta 1 minus u over u infinity dy and remember eta is y times under root u infinity by nu x. Therefore, d eta is simply dy times under root u infinity by nu x. So, if I change these values to eta, then you will see this becomes 0 1 minus f dash into d eta into under root nu x by u infinity and this will become eta at delta. Eta at delta is what I shall call delta star will be simply, sorry eta at delta will be simply delta under root u infinity by nu x, which is nothing but eta max in our case in a way or less than eta max. So, you will see now I get delta 1 under root u infinity by nu x equal to 0 to infinity or whatever value you want to call it 1 minus f dash d eta. If I manipulate this delta 1 by x into under root u infinity x by nu and that is equal to delta 1 by x r e x to the half equal to 0 to infinity 1 minus f dash d eta. This quantity I define as delta 1 star and analogously delta star delta by x u infinity by x by nu. This is the value of eta at f dash eta equal to 0.99. Delta 2 star likewise is delta 2 by x under root u infinity by nu x is equal to 0 to infinity f dash 1 minus f dash d eta and in this case in actual numerical integration after the converged solution is obtained, I simply replace this infinity here by eta max and this also I replace by eta max. No harm done. So, I get values of delta 1 star delta 2 star delta star and I recover the value of c f x as 2 f double prime 0 r e x double prime 0. So, these quantities are evaluated after the solution has been obtained. I have obtained several solutions for u infinity equal to c x raise to m which is our velocity profile. So, all positive values of m represent accelerating flow, all negative values of m represent the decelerating flow or the flow with an adverse pressure gradient. This is the flow with a favorable pressure gradient. Remember, beta is equal to 2 m over m plus 1. This is the wage angle and I also mentioned the wage angle. m equal to 0 represent the flat plate solution and I already showed you 0.33 is the value of f double prime delta star where f dash infinity is 0.99 turns out to be 4.9 delta 1 star is 1.7 to 7 delta 2 star is 0.663. As the flow accelerates, then compared to m equal to 0, you will see the boundary layer thickness is reducing as it should because the flow is accelerating now and therefore, the viscosity affected region is thinner, but in a decelerating flow, the viscosity affected region becomes thicker. Thinner the boundary layer thickness, sharper is the velocity gradient and that is what you see here. There is an increase in the f double prime 0 value which would increase the skin friction coefficient. On the negative side, very interesting thing happened that at value of m equal to minus 0.091 which corresponds to beta equal to minus 0.2, f double prime is 0. This is the separating boundary layer and its thickness will be 7.42 compared to 4.9 for a flat plate boundary layer. If you decrease the value rather of m still further, you will simply get a recirculation region near the wall and of course, that is not admissible in boundary layer theory. Incidentally, look at this solution for m equal to 0 with a skin friction coefficient which is 0.33, f double prime 0 is 0. This boundary layer has been measured for its skin friction and velocity profile. The velocity profiles were made by Russian scientist Niquoratze in 1942 and by Lipman and Dhawan in 1951 and since then several others have verified these measurements. On a boundary layer development on a flat plate, all of them predict excellently these values of 0.33 and 4.9 which is the boundary layer thickness. I have already showed you that for m equal to 0, there is an exact solution which matches very well and so does it match with the experimental data for m equal to 0. Experimental verifications have also been done for accelerating boundary layers and decelerating boundary and separation is indeed predicted when m equals minus 0.091 which is the separation pressure gradient parameter. The commons are for m equal to 0, delta star is about 0.5 and f double prime 0 is approximately 0.33. For an accelerating flow, m is greater than 0 and the boundary layer thickness now reduces but the skin friction coefficient increases. You can see the velocity profiles obtained for m greater than 0 compared to the value of m greater than 0. You will see for m equal to 0, the boundary layer thickness is indeed just about 5 or say 4.9, 4.9, 4.92 whatever value but the boundary layer thickness is the velocity gradients are indeed very sharp as the acceleration takes place. On the negative m side, m equal to decelerating flow, the boundary layer thickness goes on increasing and at m equal to minus 0.91, you see the velocity gradient at the wall. This is eta, this is u over u infinity is in fact 0. You can see the 0 gradient very clearly. If you had to reduce m further, separation will definitely occur. This is the threshold value for separation to occur. What is our conclusion then? Adverse pressure gradient causes flow thickening compared to a flat plate boundary layer whereas the favorable pressure gradient causes flow thinning. Now, we shall look at effect of suction and blowing. Now, remember our B F is defined as V w divided by u infinity into r e x to the half and that should be a constant. That is what we have said. So, V w over u infinity under root u infinity x by nu should be constant or V w into under root x by nu infinity should be a constant. Therefore, V w should be proportional to under root u infinity by x or that is equal to proportional to x raise to m minus 1 divided by 2 and that is what I have shown here. That is what I have shown here. V w should be proportional to x raise to m minus 1 by 2. Now, I am only going to consider two cases of m. The solution can be obtained for any value pressure gradient, but you will see now that if m is equal to 0 that is on a flat plate, if I want to obtain similarity solution for which m is equal to 0 u infinity is a constant, then V w must vary for m equal to 0. V w must vary as under root x. For m equal to 1 V w must be constant because m is equal to 1. Similarity solutions are possible for these two cases only when in case of m equal to 0 V w varies as 1 over under root x and for m equal to 1 V w is constant. Both these solutions have great relevance in gas turbine cooling technology because near the leading edge of the blade, the cooling air is injected through the nozzle through the leading edge at almost constant rate. Therefore, m equal to 1 case is very well taken care of by this V w as very constant whereas, along the suction side there is a region where there is almost a constant pressure gradient m equal to 0 and that is where you inject fluid at decreasing rate as 1 over under root x. Therefore, these solutions have affinity to the situations obtained in gas turbine blade. Remember this is the leading edge. So, this corresponds to m equal to 1. So, here are the solutions for flat plate m equal to 0 and I have varied B f on the negative side up to minus 2 and on the positive side. B f less than 0 represents suction. B f greater than 0 represents blowing into the flat plate boundary layer. Now, you can see what suction does. Suction reduces the thickness of the boundary layer and therefore, the skin friction increases. Thinning of the boundary layers always increases skin friction. Along with delta star even the displacement thickness and momentum thickness is also reducing and as a result the skin friction coefficient is increasing. On the blowing side obviously, the boundary layer thickness compared to no blowing would go on increasing as you can see. With a decrease in skin friction and in fact, at point B f equal to 0.612 shear stress equal to 0 occurs which means separation has occurred due to blowing. We shall see now the case of m equal to 1. Very similar suction increases the skin friction whereas, blowing reduces skin friction. But notice that in case of a stagnation point flow no matter how hard you blow even up to B f equal to 1 there is no indication of separation at all. It will require a very large value of B f in fact, unreasonably impractically large value of B f to really cause separation in a stagnation flow. That is understandable that if I want to cause there is an oncoming flow and if I want to cause by separation of the boundary layer then the flow I must blow almost at the same rate as the oncoming flow is. So, even up to B f equal to 1 separation does not occur. Here are the velocity profile for the flat plate and stagnation point flow. Let us look at the stagnation point flow first. You will see B f equal to 0 curve is there and the velocity has become almost equal to 1 at m eta equal to about 5. Because of the suction negative values of m minus 0.5 and minus 1 there is a thinning of the boundary layer with increasing gradient of velocity near the wall. So, that increases the skin friction. On the blowing side however, there is a thickening of the boundary layer but even with B f equal to 1 there is no indication that separation would occur. On the other hand, if you look at the solutions for flat plate boundary layer you will see that now this is the solution for B f equal to 0 which is going up to 5 and these are the solutions for suction obviously thinning the boundary layer. But look at what happens for 0.5 the boundary layer thickness is almost become double of what it was for B f equal to 0 and at 0.612 you get a zero gradient velocity profile at the wall indicating a separating profile. We have seen how velocity boundary layer solutions are influenced by the pressure gradient parameter m as well as the suction and blowing parameter B f. How thinning and thickening of the velocity boundary layer takes place under the influence of these? We have also seen the merits of shooting method. It involves selection of eta max and so once a program has been written for the general purpose we can go on varying the values of m and beta f as per wish and can generate a large number of solutions for further use.