 So, welcome to today's lecture on multi-phase flows. What we will do today is look at another problem in fluid mechanics and the idea is that this particular problem will be different from the earlier problem of natural convection in the sense that the earlier problem we only had one phase okay and although towards the end we will discuss the problem of the hypothetical problem of the liquid with the being surrounded by gas we did not allow the interface to deform the interface was remaining flat. Today what we will do is we look at a problem where we will allow the interface to deform okay and what that means is you would have to use things like the kinematic boundary condition which we derived a few lectures back and also the other boundary condition is the normal stress boundary condition. So in that sense this problem is one level more complicated because we are going to consider a truly multi-phase flow problem with 2 liquids and an interface which is actually deforming okay and the idea is that we would make assumptions again to try and get a analytical site into the problem. This problem is called the Rayleigh Taylor problem okay and this is essentially an instability which is going to be driven by a density stratification. So for example, so when you have a density stratification heavy liquid on top of lighter liquid is inherently unstable you all know that if you have 2 liquid layers and these are immiscible liquids okay immiscible. So you can think in terms of water which is usually denser than most of the other liquids lying on top of another organic solvent okay. So this is inherently unstable and what is going to happen is the water will have a tendency to come into the phase which is below and then the oil phase is likely to rise up. So this and if it had been the other way if the lighter liquid is on top it is going to stay as it is okay. So what this means is the instability is going to be driven primarily by the density differences and that is something which we have to make sure we include in the model okay. So basically what I am saying is density differences drive the instability and this has to be retained this has to be included the model. So this is the geometry that we are going to be looking at this is liquid 1 which has density rho 1 and viscosity mu 1 this is liquid 2 which has density rho 2 viscosity mu 2 and we have a situation where this liquid is just lying on top of the other I mean you have been very carefully you added this other liquid and then these 2 liquids are having this flat interface. What we want to do is we want to ask the question if I were to disturb this interface how exactly is the system going to behave okay and x, y and z as my coordinate axis z is the vertical direction and to keep life simple what we will do is we do not want to have to worry about all these boundary conditions. So we are going to assume that the liquids are going to extend to infinity in all directions in the x direction in the y direction and the z direction okay. So the liquids extend to infinity in the x, y and z directions and what I want you to keep in mind as we work out this problem is the closed analogy similarities with what we did for the Rayleigh-Banard problem the natural convection problem because I think that is basically what you are going to be doing whenever you are solving a problem okay. So the analysis procedure is the same only thing is now we have to worry about this thing interface reflection. So the additional feature here is we allow the interface to deform okay and which we did not do in the earlier problem. So what that means is we need to use the kinematic boundary condition okay for instance. So that is the level of complexity. What we will do is we want to consider as usual a base state and then give a perturbation around the base state right. So what is going to be the base state? The base state is going to be one where everything is at rest this liquid is at rest and that liquid is at rest that means all the velocity components are 0 okay. So the base state is that of a static fluid. I mean both the phases are static which means I am going to say u1ss equals v1ss equals w1ss equals 0 subscript what tells me this liquid subscript 2 tells me it corresponds to the velocities of the second liquid okay. And similarly that means that is the state I have one liquid rushing on top of the other completely stationary okay. And what we want to do is we ask the question is this stable or unstable? Of course you already know the answer in the sense that if the guy is denser at the top you expect it to be unstable okay. But then we want to go through the calculation and see how exactly is this instability going to manifest or is there a condition which comes because in addition to the density difference there is also going to be a surface tension of this interface at this interface which we need to look at. So does this have a stabilizing influence? So we are going to look at the influence of the surface tension the density difference and yeah maybe even the viscosity but to begin with we are going to assume that the viscosity is not really going to play an important role because viscosity is a friction and what you are focusing on is an instability which is driven by what is happening at the interface. So what is more important is for you to include the surface tension effect. So to begin with what we will do is we will assume that the 2 liquids are inviscid because what viscosity will do is only going to slow down things. So viscosity is going to possibly change your growth rate of your disturbance and it is going to modify the growth rate it is not going to really change whether it is going to be positive or negative okay but that we will see. So we will the base state how do you find the base state? We found the velocity being stationary but we need to get the pressure gradient okay. So if you want to simplify the equation of motion in the x direction okay what do you have? Rho du by dt equals minus dp by dx plus mu del square u okay. Since the velocity is 0 this reduces to dp by dx being 0 okay. At uss equals 0 we get dpss by dx equals 0 which means pss is a function of y and z only because there is no change in the x direction okay. If you look at the y component similarly in the y direction we get dpss by dy equals 0 okay because the gravitational field is in the z direction and what does this imply? This implies that pss is a function only of z okay. So in the base state that is perfectly understandable because you only have the pressure gradient in the vertical direction in the horizontal direction there is no pressure gradient okay that is what you conclude from these equations of motion and as far as z direction is concerned what do we get? We would get 0 equals minus dp1ss by dz plus z I am showing it going upwards okay. So the gravitational field is in the negative z direction minus rho 1g this is for the first liquid okay and that means the first liquid is extending from minus infinity less than z less than 0 whereas in the second liquid I have which extends from. So this is clearly the hydrostatic pressure gradient which is what everybody understands. So all I am going to do is I am going to integrate this out and I am going to get that pressure is going to vary linearly with z. So from the first equation what do I get? p1ss equals minus rho 1gz plus a constant c1 okay p2ss equals minus rho 2gz plus a constant c2 okay. So that is basically what I get and you know that the pressure as you go up has to decrease and that is what is happening is becoming as z increases is becoming more and more negative. Now we have to determine these constants c1 and c2 okay and for that what we need to do is have some kind of a boundary condition. So we are going to take we take at z equals 0 p1 equals p2 of course at the interface p1 the pressure will be equal p1 equals p2 equals 0. Now what is the motivation for this? Of course you can take it to be any arbitrary constant okay. The idea is that flows are going to be driven by pressure differences okay. So it is not the actual value of the pressure which actually is going to drive the flow. When you have a pipe flow if you have an inlet pressure of 80 atmospheres, outlet pressure is 70 atmospheres you are going to have a drop of 10 atmospheres okay and with incompressible liquid you have a particular pressure gradient that is the flow. Suppose you have 50 atmospheres and 40 atmospheres you again have the same pressure gradient the flow is going to be the same because the pressure gradient is what is important absolute value of the pressure does not matter. So idea is that you can choose one of these pressure points as a reference point and you can calculate what the other value is okay. So that is basically what we are doing here. We just choosing at z equal to 0 the reference value of the pressure to be 0 okay and keeping this as the reference value we are going to find the relative value of the pressure at the other points and then see what is going to happen. So this you can just say is a reference value okay and this is permissible since pressure gradients drive the flow and the actual pressure does not matter okay. So with this simplification at z equal to 0 pressure is 0 what do I get I get basically C1 and C2 will be 0 okay. This yields P1 equals minus rho 1 gz and P2 equals minus rho 2 gz. Now that we found the steady state the steady state is characterized by the velocity and the pressure and we found that. We need to now find out a bit stable. So what do we do to find out the stability of this state? We have to give this perturbation right. So we write the actual variables. So okay what I am going to do is I am going to make one more assumption here. We will assume the 2 liquids in message. So the way I would justify this is I like to make this assumption and see what the analysis yields. If I am not happy I will just go back and relax this assumption include the effect of viscosity and go through with the analysis okay. So that is what you should be doing whenever you do any problem. I mean you start with a simple problem and then see if you are getting any insight. If you are not getting any insight you made too much of a simplification. You add a level of complexity and then you keep on building. It would have been stupid for me to assume that the 2 densities are equal and keep the differences in the viscosity. Because what I am expecting to drive the flow is the density difference. So I am going to retain the density difference but just to make sure that the algebra becomes easier and I am going to neglect the effect of viscosity okay. So this basically is going to simplify mathematically okay and then do my calculation. I get some result. At the end of the day the result is anywhere close to the actual problem then I say fine I mean maybe this assumption was not very bad after all. Maybe the effect of viscosity is not important okay. But if it turns out to be you know different from what I actually seen in experiment then I come back and say well maybe this is the problem. Maybe I assume this thing to be invisible and actually the viscosity effect is important and I have to include it okay. That is the only way to go about doing it okay. So right now my motivation is I do not want a second order equation. So I just want to just keep it simple and so that is the motivation. And let us see if I can get away with it okay or get some insight at least into the problem. So this simplifies the mathematics okay and if the predictions are close to the experiment then we can possibly justify the assumption okay. Of course just because it is close to the experiment does not mean it is right. I mean maybe you just got lucky okay. So I mean you have to be careful. If they are not matching then you know for sure something is wrong. If it is matching you never know if you are right okay that is always a problem okay. So let us not just say that just because this is matching the experiment everything is perfect. Everything is likely to be perfect not necessarily perfect okay. So now what do we do? We find the stability by giving a perturbation. Same step find a steady state. Find do a perturbation. Find the linearized equations and solve. That is what we did last time. We are going to do that is what we have been doing for the last 3-4 classes okay. So what is this thing u1? The actual variable is written in terms of u1 ss plus epsilon u1 tilde. Similarly for everything now remember once I denotes the fluid first fluid or the second fluid similarly for all the variables okay. Similarly for all variables and so is p1 tilde and u2 equals u2 ss plus epsilon u2 tilde and so on okay. And just to tell you that these things are infinitesimal I put that epsilon in front of it. So what do we do now? We substitute all this in my governing equations because the governing equations are going to be valid for u1 v1 u2 v2 okay. The actual variable I am going to go back to my equation of continuity and write whatever u is there as u ss plus epsilon u tilde okay. And we have to do this for all the equations. So the equation of continuity implies d by dx of u plus d by dy of v plus d by dz of w equal to 0. This is for the first fluid and I am going to substitute u1 ss u1 as u1 ss plus epsilon u1 tilde and this means and u1 ss is 0. So when I substitute this back inside here at order epsilon what do I get? d u1 tilde by dx equals 0. Now we go to the equation of momentum in the x direction. The name is Stokes equation in the x direction. What is that? That is rho1 d u1 by dt plus u1 d u1 by dx plus dy dx of maybe we should just write this as p1. Let me just write this as p1 and then I will substitute this thing here. This is in terms of the actual variables. I have not done any perturbation. Now I am going to substitute for all the u1 ss plus epsilon u1. So what do I know? This is going to become rho1. u1 ss of course is 0. So I get u1 tilde by dt multiplied by epsilon because u1 is u1 ss plus epsilon u1 tilde. So u1 ss is 0. So u1 is epsilon u1 tilde. So I have epsilon u1 tilde here and what about this guy? This is going to give me epsilon u1 tilde times the derivative of epsilon u1 tilde which means this is epsilon squared okay. So this is a second order term and this is not going to contribute. This is not going to contribute because this is a second order term. This guy is not going to contribute again for the same reason. All I am saying is this is epsilon squared u1 tilde d u1 tilde by dx plus epsilon squared v1 tilde d u1 tilde by dy plus epsilon squared w1 tilde d u1 tilde by dz okay and all these guys just go off because they are of higher order okay and what does this become? This d by dx of p1 ss plus epsilon p1 tilde. We know that dp1 ss by dx is 0 okay. So that is not the p1 ss is 0. It is that the derivative is 0 and therefore this reduces to at order epsilon I get rho1 b u1 tilde by dt equals minus dp1 tilde by dx. Now you can do the same analysis in the y direction. So let me just write that thing down neatly a bit. Let me write those 2 equations neatly because I want to preserve this so that I do not make any mistakes. This is just the equation of continuity okay. It is already derived and then the other one is dx. If you did the thing in the y direction either you are going to agree with me or if you do not agree with me you have to work out the problem and they said you do not agree with me. Another equation in the y direction you would get rho1 dv1 tilde by dt equals minus dp1 tilde by dy. Everything is going to be the same okay. You will get the convective terms are going to give you a second order term. So that is not going to contribute. These guys are going to contribute at order epsilon and in the z direction you only have one small complication in the sense that you have the gravity term coming okay. But then the gravity term remember I will just explain that and this is similar to the energy balance we had for the Rayleigh-Bernard problem. In the z direction you will have rho1 dw1 tilde by dt plus the order of epsilon square terms if you do not worry about will be minus dp by dz of p1 ss plus epsilon p1 tilde minus rho1 g and at steady state dp1 ss by dz is equal to minus rho1 g okay. So this is going to balance of that guy and so that goes off okay. I did not write the gravity term in the x and y direction because the gravity does not exist in the x and y direction. Here the gravity does exist but this gravity is going to balance dp1 ss by dz because dp1 ss by dz remember was equal to minus rho1 g. So this again simplifies to minus so this again so basically this and this are 0 from my steady state and together okay is not that rho1 g is 0 is that together there is 0 and I get the same equation for w1 tilde also. The point I am trying to make here is you guys have to sit down and make sure you do each and every term properly and then do the calculation okay and I did not want to just say you get the same equation because there is a small settle point here. This is all for the one phase. You will have similar equations for the other phase okay. You have same equations for the other phase. So basically these are your equations which are going to tell you how the perturbations and the pressure are basically related to each other okay. We have similar equations in phase 2. Now in the really very problem in fact if you remember I tried to write the expansion first I had to convert it to the partial differential equations to ordinary differential equations and then reduce the number of dependent variables okay and then we found that we have got stuck because I had assumed the sin alpha x dependency and then the 2 velocity components were actually out of phase sin and cosine and then you had a problem. So I really could not proceed and then somebody said we should use e power i alpha x as a way out. So that is what we will do now. So that what we did was we reduce the number of dependent variables first u v w p and then we converted it to ods. But now what I am going to do is I am going to do the expansion in terms of the independent variables first and then do the elimination of this thing. So what we will do is we are going to seek the solution. So both are actually equivalent. We will seek the solutions as u1 tilde remember what is u1 tilde function of x, y, z and t okay. We are going to seek this as u1 star of what variable and need to impose the boundary conditions in the z direction. What boundary condition? The one of the interface some continuity of velocity whatever it is okay. So I need to have that direction to be determined using the boundary condition okay. So that is the guy I am going to keep here and in the other 2 directions I have assumed them to be spanning to infinity in the x and y directions okay. So I am going to assume look for periodic solutions there but instead of looking for things in the form of sine and cosine I am going to look for things as e power i alpha subscript xx multiplied by e power i alpha subscript yy and the time dependency is on the form e power sigma t okay. So because these equations are linear I am assuming that this is the form of the disturbance. Our objective again is to get this relationship between the growth rate and the wave number because that was the thing we finally got even in the Rayleigh-Bernard problem sigma versus alpha Rayleigh number versus this thing. So what does alpha x represent? It represents some kind of a spatial frequency in the x direction. Alpha y represents the spatial frequency in the y direction. Alpha x is in the x direction, alpha y is in the y direction okay and that is my growth rate here sigma. So this is exactly what we did last time but the point is instead of using sine and cosine I am just using the general form of the Fourier transform not a Fourier sine transform not a Fourier cosine transform. What I am going to do is I am going to substitute this in and I am going to make the same thing for all the variables okay u1, v1, w1, u2, v2, w2. So all variables follow this form okay. So here I know all the 8 variables I am just leaving it as it is. So now what we need to do is substitute and get ordinary differential equations for u1, v1, p1 star okay and the way we do it is substitute it in the equation of continuity. When you substitute this in the equation of continuity I will get the derivative of u1 with respect to x. The derivative of u1 with respect to x is i alpha x multiplied by the entire thing okay. So this is going to give me from the equation of continuity I get i alpha x times u1 star times e power i alpha xx plus alpha yy plus sigma t okay. Then dv by dz is going when I have v I have v1 star okay sorry dv1 by dy. I have v1 star of z of course and when I differentiate with respect to y I get i alpha y. So okay and when I differentiate with respect to z for w1 I would get dw1 star by dz times e power i sigma t equals 0. Now clearly the problem which I had last time of sine and cosine coming in is not there. I have the exponential term everywhere derivative of the exponential gives me the exponential that means this is clearly nonzero and I can knock this thing off and what I have is i alpha x u1 star plus i alpha y v1 star plus dw1 star by dz equals 0. I want to come back and do the same thing for these 3 equations. I want to basically now my job is to get this ordinary differential equation and use the boundary conditions and find out what the stability condition is going to be okay. So what we are going to do now is just substitute the form for u1 here what do I get rho1 when I derivative this differentiate that with respect to time I am going to get a sigma multiplied by u1 star etc etc okay. So I am using this condition now some equation number maybe this is equation number 2 okay this is equation number 1. So this is from equation of continuity which is 1 and from 2 what do I get rho1 differential with respect to time I get sigma multiplied by e power sigma t times plus i alpha xx plus alpha yy times u1 star of z okay equals minus dp1 by dx tilde that is going to be pressure is also going to be of the same form okay when I differentiate that I get minus of i alpha x times p1 star times e power sigma t plus again the fact that such a form is admissible is coming because of the fact that this particular term is common for both and what this gives me is that rho1 sigma u1 star equals minus i alpha x p1 star that is what this is. So what we can do is we can extend the same argument to the third equation here and what we will get from 3 when I differentiate with respect to v1 I will get rho1 v1 star sigma and when differentiate with respect to y I will get minus i alpha y times p1 star okay that is the equation which I get from this thing and what about the from 4 which I have not yet written over there but that is going to be equation number 4 this guy is equation number 4 okay when I differentiate with respect to time I get rho1 w1 star sigma equals now I am differentiating with respect to z so I get minus dp1 star by dz because exponential terms do not have z is the p1 star which is unknown which is a function of z so I have that here. So now what I have done is I have converted it to an ordinary differential equation instead but I still have 4 variables the u, v, w and the pressure with the subscript 1 and the star okay and what I am going to do is I am going to eliminate okay what we can do is we have to eliminate let us say the u and the v component of velocity okay and one way for you to do this elimination of the u and the v component of velocity is we can write from this what is u1 star and from this what is v1 star in terms of pressure okay I can find u1 star v1 star from these equations substituted in my equation of continuity then I will be getting rid of u and v get everything in terms of pressure and w so I have one equation in terms of pressure and w I have another equation in terms of pressure and w and I can go back and eliminate pressure again maybe and get only an equation in w1 okay. So that is the basic steps that I am going to follow and so what do I get u1 star equals from this equation here minus i alpha x p1 star divided by rho 1 sigma okay v1 star is minus i alpha y p1 star divided by rho 1 sigma that is what I found from these 2 equations. Now for u1 star and v1 star I am going to substitute in this equation okay so I have i alpha x multiplied by u1 star which is again an i alpha x with a negative sign minus i alpha x whole squared by rho 1 sigma minus i alpha y the whole squared sigma plus dw1 star by dz okay. So equals 0 I hope that is fine all I have done is doing a little bit of algebra here and this is going to be minus i squared so that is going to be plus 1. So I get alpha x squared I am going to take the rho 1 sigma along with my dw1 by dz okay so I get minus alpha x squared or is it plus it is plus alpha x squared plus alpha y squared times p1 star plus rho 1 sigma dw1 star by dz equals 0. So this is the equation which relates pressure and w1 okay. I have also another equation which relates pressure and w1 right here and what I am going to do is I like to keep my velocity because I like to have my conditions on my velocity, my kinematic boundary conditions. So rather the eliminate velocity I am going to eliminate pressure. I eliminate pressure by differentiating this with respect to pressure I will get dp1 by dz second derivative and I will substitute for dp1 by dz from this equation okay. So that is what we are going to do differentiate that the last equation with respect to z and I get alpha x squared plus alpha y squared times dp1 star by dz plus rho 1 sigma dw1 star by dz equals 0 correct yeah and you know dp1 star by dz is minus of rho 1 w1 star sigma. So this gives me minus of 0. So clearly rho 1 sigma cancels off and what I have is d square w1 star by dz squared minus alpha squared w1 star equal to 0 okay. What is alpha squared? Square is alpha x squared plus alpha y squared. So this again is to tell you say when we did the really bad problem I assumed that things were actually 2 dimensional that we had roles. We did not have any pattern in the other direction in this one of the directions along with the axis of the role was extending. So even if I not made an assumption I would have gotten 2 wave numbers like this and finally I would have gotten a composite wave number alpha okay. So it really does not matter because that is only a mathematical complexity not a physical complexity. See what we want to do is make sure we retain the right physics in the problem. So this is my equation for the one phase if you did the same thing again for the other phase you would get similarly we get d squared w2 star and d is a squared minus alpha squared w2 star equals 0 okay that is for the other phase. One represents the one phase two represents the other phase. We know the solutions to this equation alpha is a constant right this linear equation with constant coefficient second order. We know how to solve this problem what is the solution w1 star is a e power alpha z plus b e power minus alpha z that is the solution to the differential equation w2 star c e power alpha z plus d e power minus alpha z. Now our job is to find the constant a, b, c and d and for this we need boundary conditions okay. So what are the boundary conditions we are going to have remember this first liquid is extending from 0 to infinity to minus infinity clearly what we expect is when there is some kind of an instability at the interface we expect that far away at z goes to minus infinity in the first fluid as z goes to plus infinity in the other fluid the velocity components are going to go to 0 okay. So basically you expect the velocity far away from the interface to be finite and bounded it cannot become infinite. So that is going to help us determine two of these constants for example as z w2 is from 0 to infinity. So as z goes to infinity I want the velocity to be bounded which means this guy should be present because this is minus sign and this will go to 0 this has to be absent. So that tells me c has to be 0 here okay similarly in the other fluid in the lower fluid we will have z is going to minus infinity should be present okay this will be present and this will be absent something like that. So is not that right that is right is not it I thought I said the same thing twice. So that basically helps you determine two of the constants we need to determine the other two constants and that is where the boundary conditions at the interface come in and that is where we are going to use the kinematic boundary condition and the one more condition. So the other condition which we are going to use is the normal stress boundary condition. This normal stress boundary condition is going to be preferred over the tangential stress boundary condition. The reason for that is that we made this thing inviscid see the two conditions which have to be satisfied at the interface both the normal stress boundary condition as well as the tangential stress boundary condition. How is it that we do not need both we need only one the reason is we assume that the fluid is inviscid. So we actually had a second order problem but because I have assumed it to be inviscid my problem from second order has become first order. So I need to let go of one of the boundary conditions. The boundary condition which I am going to let go of is the tangential stress boundary condition because there is no continuity of tangential stress because normal stress is going to be present even in the absence of viscosity. So I retain the normal stress boundary condition and let go of the tangential stress boundary condition. What we will do is we will use those conditions get these constants and get the distribution curve tomorrow.