 In today's lecture what we will do is we will formally derive the condition for the onset of viscous fingering. So we just laid the foundation for Darcy's law etc. in the last class and we will go about the formal procedure. So the procedure is basically the same as what we have been following earlier. We write down the model equations, find the base state, do the linearization and get a relationship between the growth constant and the wave number okay and that basically is what we have been doing to find out under what conditions and instability can occur. So we will follow the same procedure again but then the specific problem that we are going to apply to is the viscous fingering problem okay and the idea is that we have let us say a vertical geometry. So we are having viscous fingering or the Saffman-Taylor problem. The z axis is vertically upwards okay, gravity is downwards and let us say this is fluid 1, this is fluid 2, this has properties mu1, rho1 and this has properties mu2, rho2 and let us just solve this problem okay. Actually you can have different permeabilities also in both the geometries but I think what we will do is we will keep the permeability the same throughout okay. So we will just say that the permeability here is k and the permeability here is k okay. So now what we want to do is we are going to use Darcy's law and Darcy's law basically tells you that the z component of velocity, the z component of velocity is uz is going to be given by – k divided by mu times the gradient of p plus rho gz because this is plus here because g is acting in the negative direction. You have to check with your notes in the last time what we had, it was that the component of g was coming. So now g is in the negative z direction and therefore this is plus rho gz okay. So when I go back to the Navier-Stokes equation I have –db by dz plus rho gz that is the minus sign and so I can group these two together. And this I can write as plus gradient of phi that is what we did last time. So basically what I am doing is I am writing the velocity component uz as a gradient of a potential and the potential is essentially we compare these two equations. There is also a scalar – k by mu times p plus rho gz that is your potential okay. Now the base state whose stability we are interested in finding out. The base state is base state has a flat interface and this guy is moving let us say at a constant velocity. So you are pumping liquid and let us say that the velocity is uniform across the channel. So we will just look at the velocity as being uniform and some capital V okay. So we have a uniform velocity in the z direction okay. So what is going to happen is the base state is one where the interface keeps moving. So the interface keeps moving like I said earlier is going to be an unsteady state problem. So what we can do is we can convert this to a steady state problem by just working in a reference frame which is also moving at velocity v. Supposing you are sitting on this reference frame and let us say at time t equal to 0 this is at z equal to 0. I am going to define a new coordinate z bar which is z – vt. So this is my moving reference frame okay moves vertically upward at a velocity v. Now so in this reference frame what is going to be the location of the interface? The location of the interface is going to be given by z bar equal to 0 okay. So z bar equal to 0 represents the location of the interface okay. So basically the interface is moved up but then when it moves up some distance I am looking at z – vt as the location z bar is 0 always and in this moving reference frame clearly what is the base velocity? The steady state velocity whose stability your interest in finding out what is that? That is also 0 because you are moving the fluid is moving as a plug v. You are also moving along with it. So in that reference frame the base velocity of the steady state velocity is 0 okay. The steady state velocity equals 0 in this reference frame right. So now we do the usual stuff which is do the linearization about this base velocity okay which means the u vector I am going to write it as the base velocity which is 0 plus epsilon times u tilde okay and this of course is a vectorial tilde. Now just like we say that the actual equation is going to be u is going to satisfy Darcy's law. So that means u tilde will also satisfy Darcy's law okay and therefore we have what? u tilde is going to be given by I am going to work in terms of this potential gradient of phi tilde okay. That is there is a perturbed potential phi tilde whose gradient is going to give me my u tilde okay. So u tilde represents the perturbed velocity which is of order epsilon and phi tilde will also be of order epsilon remember that. So phi tilde is of order epsilon and what you also have to keep in mind is that you actually have u1 u2 phi 1 phi 2 phi 1 for the potential in one liquid phi 2 for the potential in the other liquid. So that is something you need to keep in mind. So I have for each liquid each phase I have ui tilde equals gradient of phi i tilde and this is of order epsilon because these are my perturbations okay. Clearly u tilde has to satisfy my continuity equation for divergence of u tilde has to be 0, divergence of u has to be 0 for divergence of u tilde has to be 0 at the order of epsilon. So divergence of u tilde has to be 0 and I can combine that with this for each of the phases and which means that divergence of u tilde is del square of phi i must be 0 for each phase okay. So basically what this means is del square of phi i tilde must be 0 for i equals 1, 2. This is my linearized equation. Now we would like to actually solve this equation, actually partial differential equation x, y, z and time okay. So although we have kept the thing as independent of time in the differential equation, the time dependency is going to come through what? Through the boundary condition. The boundary condition when the interface is going to get deflected, you can have the location of the interface changing with time and you will be using the kinematic boundary condition which has the time derivative. So through that boundary condition, the time dependency comes in okay. So whenever you are talking about the linearized problem, it has to be an unsteady state problem. Only then you will know whether it is growing with time or not. Now we have to go back and solve this analytically right. So we have to make those simplifications which we have done earlier which is we are going to assume that the system extends to infinity in the x and y direction okay. System will extend to infinity in the x and y direction. So which means I can now seek solutions which are periodic in x and y of the form e power i alpha x i alpha y y okay and what I will do is I will go back to my interface and z bar is going to be when I am going to be deforming it, z bar equal to 0 is my base state okay. So z bar equals epsilon times h of x, y, t is the location of the perturbed interface okay. z bar equal to 0 is the base state and I am going to give a perturbation it is an epsilon h. So if you want to write it in terms of this implicit function, you would write it as z bar minus epsilon h of x, y, t equal to 0 because this is the starting point for getting the kinematic boundary condition and all that okay. So now we are going to assume h to be of the form h star e power sigma t e power i alpha x x plus alpha y y. So the disturbance I am giving to the interface is periodic. I can give an arbitrary disturbance, I am resolving it in different Fourier modes okay and I am trying to find out which of these Fourier modes is going to increase or decrease. So like resolving a vector on different components okay. So this is it and this is the growth with time and this is h star which is the amplitude. So if this is going to be the form of the disturbance of the boundary, then clearly my phi i also is going to have the same periodic form okay. As far as x and y is concerned, only the z direction I have to find out what is going to happen okay. So as far as the phi i is concerned, phi i tilde is concerned. So phi i tilde is going to be of the form f i of z bar, t e power i alpha x x plus alpha y y okay. So remember phi i is pde in x, y and z. x and y components are like this. z and time dependency is captured here. I am going to substitute this form in my del square equation. This is what we have been doing earlier. When I substitute this in the del square equation, what do I get? I get d square f i by dz square okay. In fact, this I should possibly write it as dou square. Minus, when I differentiate with respect to x 2 times, I will get minus alpha x square minus alpha y square. So minus alpha x square plus alpha y square times f i equal to 0 okay. Is this clear? All I am doing is substituting this in the del square equation because del square phi i is equal to 0, we get this for each of the i's. And clearly I am going to just combine alpha x square plus alpha y square into some alpha square okay. So alpha square define alpha square as alpha x square plus alpha y square. Let this be the case, then I have d square f i by dz square minus alpha square f i equal to 0, which means f i remember will be of the form, some constant into t because that constant will now be a function of time because that is what where my time dependency is being absorbed okay. So I am going to call this g i of t times e power alpha z. In fact, I should be slightly careful. Let me just do one thing. Let me just use something else. So clearly this is the solution to this equation. I mean this is a linear equation with constant coefficient assume solution of the form e power alpha z. Put it here, this is what you get. So now I need to find out, I need to get these constants. Clearly very far away, I mean I am having also a situation where there is no boundary condition in the z direction and minus infinity and plus infinity. I want things to be bounded and minus infinity and plus infinity right. So at z equal to plus infinity and that is where my first fluid is. First fluid is on the top. F1 is bounded. So that means this guy has to be 0 implies a1 of t equals 0. Here a1 of t is dot 0 and z equal to actually should be z bar no. This should be z bar. Yeah, this should be z bar because I am working in the moving reference frame. A1 of t must be 0 okay and clearly at z equals minus infinity, f2 is bounded which implies b2 of t equals 0 okay. So essentially what this means is it means phi1 of t equals what? F1 is b1 of t e power i minus alpha z e power i alpha xx plus alpha yy e power sigma t okay. I mean phi1 I wrote as f times this. f is the z dependency I am now breaking up into this exponential in this and that okay. And phi2 of t is a2 of t e power alpha z bar. I need to remember to put this bar on top okay. So what I have done is just use the fact that phi is defined as follows and f i is going to be no. There is no e power sigma t. Yeah, so f i is that and the time dependency I am yet to find out. Yeah okay, that is right. So I think everything is fine. So this is valid for z bar greater than 0. This is valid for z bar less than 0. Now what? I am going to find these 2 fellows b1 and b2 right and of course I need to now start using my boundary conditions. My boundary conditions which are going to be necessary are my kinematic boundary condition and my normal stress boundary condition. Like I told you although we have viscosity here we have a potential flow situation okay which means it is something like an invasive flow but then viscosity is present. So I am trying to get the best of both worlds okay. So we have to use the normal stress boundary condition and the kinematic boundary condition. So rather than me go back and derive what this kinematic boundary condition is on first principle you know how to do it. The kinematic boundary condition comes by looking at this and saying that the material derivative is 0 okay and then you can equate terms of order epsilon. So now what we will do is the kinematic boundary condition boundary condition gives what the vertical component of velocity equals dh by dt okay. The vertical component of velocity is now d phi 1 by dz that is my perturbed. In my perturbed reference frame v is my base state which is when things are flat. So if I am looking at quantities which are of order epsilon okay what do I get? At order epsilon I will get the perturbed velocity in the z direction w1, w2 okay will be related to dh by dt. Now this is something you guys have to go and derive. I am not going to sit and do this you have done this for the earlier problems okay. So just find the put the kinematic boundary condition using the same method as before and what you will get is the vertical component of the velocity is d phi 1 tilde by dz equals d phi 2 tilde by dz because this is remember velocity is equal to gradient of potential that is how I had. So the velocity in the z direction will be d phi 1 tilde by dz that is my vertical component that is my these 2 have to be equal at the interface okay and that must be equal to dh by dt that is my kinematic boundary condition okay. And this is equal to dh by dt and this is at order epsilon because h is already at of order epsilon this is also order epsilon. So now I am going to substitute for h in terms of this equation here and I am going to get when I differentiate this with respect to time I am going to get sigma times h star times e power sigma t times e power i alpha xx plus alpha yy that is my d phi 1 tilde by dz this is my d phi 2 tilde by dz. In fact remember guys this is not only phi 1 of t this is also a function of x, y and x okay this is phi 1 completely okay. So what I am going to do is I am going to use I know phi 1 and phi 2 I am going to find out d phi 1 by dz from here I will get minus alpha times that. So this is yeah things are fine I am going to substitute here d phi 1 by dz I am going to find out d phi 1 by dz and find out b1 in a2 okay. So b1 of t so what is d phi 1 by dz bar b1 okay remember that is the function of time and then there is a minus alpha e power minus alpha z bar okay times e power i alpha xx plus alpha yy that is the derivative I have okay must be equal to sigma h star e power sigma t e power i alpha xx plus alpha yy okay. Now clearly this guy is the same as that so that cancels off since we are looking for nonzero there and the question I have now is z bar remember is of order where is the kinematic boundary condition going to be applied it is going to be applied at z bar equal to epsilon h the kinematic boundary condition is applied at z bar equal to epsilon h okay. So now if I want to substitute this at z bar equals epsilon h and what I would have is something containing epsilon here. So basically you do a Taylor series expansion and that is your domain perturbation method that you saw earlier. So you have e power minus alpha epsilon h you evaluated at the base state plus epsilon h right. So what I am saying is e power minus alpha z bar can be written as this thing evaluated at 0 which is 1 okay plus the derivative of this which is minus alpha times e power minus alpha z bar multiplied by z bar okay. So point I am trying to make here is basically I am just explaining the domain perturbation method which you people have seen earlier this term is of order epsilon. So when I substitute this here this is not going to contribute and I am going to essentially evaluate this at z bar equal to 0 that is the story okay. So I am going to evaluate this at z bar equal to 0 although ideally I am supposed to evaluate this at z equals epsilon h okay. So when I evaluate this at z bar equal to 0 I mean this is basically approximate it to this this is of order epsilon and therefore I have minus alpha b1 equals sigma h star e power sigma t okay. So b1 of t is nothing but sigma h star e power sigma t divided by alpha with a minus sign. You can do the same thing for the everything okay you can do the same thing with the other fellow d phi 2 by dz because I just use this equal to that I am going to use this equals this and then I am going to substitute this expression for a2 this minus is not there. So what I am going to get if you believe me is a2 of t equals plus sigma h star e power sigma t alpha that is the time dependency that I have okay. So let me just write this thing neatly once so phi 1 of t what is that phi 1 of t turns out to be b1 of t which is minus sigma h star e power sigma t by alpha okay that is my expression for phi 1 and phi 2 what I need to do is I have basically related my amplitude a1 and a2 in terms of h star what I have to do now is find out conditions for which h star is non-zero. So the only thing remaining for me to use is the normal stress boundary condition. The normal stress boundary condition is going to be basically staying and remember we are going to be working we are looking at the limit of surface tension not existing no surface tension okay we can include the effect of surface tension but right now for simplicity we just say surface tension is not there. If surface tension is not there that means p1 must be equal to p2 the pressures are equal at the interface that is the simple thing and for all practical purposes we are not going to worry about the normal contribution because of the viscosity because we are just saying this is something like a potential flow okay. So what I am doing is I am going to say p1 equals p2 at the z bar equals epsilon h this is the normal stress boundary condition. So what is p1? I have gotten things in terms of phi 1 and we know the relationship between pressure and the gradient right I mean what did I write? Minus k by mu times p plus rho dz right equals phi is this right please check. So if that is the case that this implies p1 equals minus mu 1 by k phi 1 minus rho 1 dz okay I am just I am now I have to make sure I put the 1 and the 2 in the right place because that is going to make a difference right. So if I make a mistake here I am going to really mess up. So p1 is going to be for the first fluid mu 1 by k1 phi 1 minus rho 1 dz and similarly for p2 what I have to do is I have to put p1 equals p2 because that is my pressure which is going to be equal. If p1 if there is surface tension then p1 equals p1 minus p2 equals sigma del dot n that curvature I have to put okay. So now I am not going to worry about surface tension this is with surface tension equals 0 phi 1 is yeah this is rho 2 yeah this is rho 2 you are right. This is the total pressure the actual pressure of the interface okay. So now what I found is phi 1 tilde the perturbation remember these guys what I found out all the written phi 1 here these are the perturbed quantities okay. So yeah I put a perturbation here but for some reason I forgot the perturbation here. So these remember are actually perturbations okay these are all perturbations. So what is the relationship between phi 1 and what is the base state phi 1? The base state phi 1 is uniform velocity V which means what is the base state phi 1? It is Vz because V phi 1 by dz is V that is my vertical government of velocity phi 1 is therefore V phi 1 is my base state plus my perturbation plus epsilon phi 1 tilde okay phi 2 is V plus epsilon phi 2 tilde and what I need here is the actual potential you understand. So what I have to do is put P1 equals P2 put these 2 guys equal I I will tell you what I am going to be doing I am going to be putting z equals epsilon h because this boundary condition is evaluated at z equals epsilon h okay. So I will get an h here h star here phi 1 and phi 2 already have things in terms of h star okay I am going to equate these 2 guys and I am going to use the condition that I want a non-zero h star and that is going to give me a relationship between sigma and my properties and that is my stability condition okay that is the plan. So that is the way we are going to go about doing this but phi 1 is also not right is not it this is there is a V z bar because the derivative of this with respect to z bar is my base velocity which is V pardon me z bar it is not there in the new reference frame what should be a constant phi 1 is just a constant in the z bar frame phi 1 is a constant in the z bar frame yeah 1 second that is not right okay I will tell you why because in the moving reference frame yeah I am not changing my reference frame I am keeping my reference frame as it is. So this is let me do one thing let me go through the algebra first and then I will come back and address this question okay. So one more question which I have to address so what I am going to do yeah I need the V because if there is no V that means as the condition has to have the V in it there is a reason for this let us come back to that so let us just do the algebra right now okay I will explain to you why the thing has to come d phi 1 by d z no what do I need to do I need to substitute this back here p1 equals p2 is mu 1 by k phi 1 plus rho 1 g z bar equals mu 2 by k phi 2 plus rho 2 g z bar okay so I am just removing the minus sign I am just saying these two have to be equal that is my normal stress boundary condition and mu 1 by k what is phi 1 I am saying it is V z bar plus epsilon times phi 1 tilde which is over here phi 1 tilde is minus sigma by h star e bar sigma p by alpha y okay equals mu 2 by k V z bar what is this phi 2 so that is these two have to be equal okay and what I am going to do is this has to be used at z bar equals epsilon h okay now the same stuff here z bar equals epsilon h I am going to substitute for h in terms of my h star exponential i alpha x x i alpha y y okay and this is epsilon equals h star e power i alpha x x plus alpha y y times e power sigma t I am going to substitute this for z bar then what do I get I get a term here which is of order epsilon okay because z bar is of order epsilon I get a term here which is of order epsilon this is of order epsilon this is of order epsilon this guy has e power minus alpha z bar and there is already an epsilon multiplying it so what I need to do is I need to evaluate this at z bar equal to 0 because the next term the Taylor's expression will give me a higher order term so basically I am going to evaluate these two terms at z bar equal to 0 that is my domain perturbation method okay and I get these guys will go off and all these guys will have e power sigma t e power i alpha x x i alpha y y so all this e power sigma t will push off you understand so basically what I am saying is when you substitute z bar equals epsilon h in all terms put e power minus alpha z bar as equal to 1 okay in the middle term because I am using the domain perturbation method and you cancel off all this e power sigma t e power alpha x x alpha y y okay what you are going to be left with this mu 1 by k times v times h star minus sigma by alpha h star plus rho 1 g h star equals mu 2 by k times v times h star v times h star plus sigma h star by alpha plus rho 2 g h star. So basically this h star occurring all of this I want h star to be nonzero because only then my perturbation is nonzero and that gives me a condition which is going to relate sigma and alpha okay we have mu 1 by k into v minus sigma by alpha plus rho 1 g equals mu 2 by k into v plus sigma by alpha plus rho 2 g okay I am going to keep all my sigma's to one side move all my sigma's to one side and I get sigma by alpha remember sigma by alpha is together multiplied by mu 1 plus mu 2 by k I am moving this guy to that side and now mu 2 by k so sigma okay let us leave this as it is okay this is the relationship between sigma and alpha the growth rate and the wave number in order for you to have a disturbance the sigma by alpha is going to be related by this clearly as alpha increases sigma increases linearly what is the condition for stability sigma must be positive right sorry sigma must be negative for stability sigma must be positive for instability that is if you are going to have this kind of a perturbation which is going to grow the viscous fingering sigma must be positive okay so basically what it means is sigma must be positive implies unstable that is fingering will occur so when is okay let us assume that the densities are equal okay clearly both density and viscosity are playing a role but to begin with to find out the effect of the individual things let us assume that the densities are equal if rho 1 equals rho 2 clearly mu 1 must be greater than mu 2 okay for fingering to take place that means if you have oil and if you are pumping water which is your second fluid first fluid is oil which is there and remember 2 is my second fluid which is water so the viscosity of oil is greater than the viscosity of water you are using water as my fluid which is most likely to happen okay then you will have viscous fingering you will have instability so basically what it means is if the more viscous fluid is being driven out by a less viscous fluid we will have fingering of course even if for example I did this problem in a vertical frame that is why this G showed up supposing you do not have this thing in a vertical frame you have this thing in a horizontal situation if you have the flow in the on the displacement in the horizontal direction then that gravity term is not going to show up okay so even if the densities are different the density is not going to make a difference for a horizontal flow but this is the guy which is going to decide whether you are going to have fingering or not so essentially if you have a heavy viscous liquid where you are trying to push it through a push displaced it using a less viscous liquid you will get fingering but if you have a water which is let us say less viscous and you are trying to displace that with oil which is more viscous the interface is going to remain flat okay so basically I just wanted to elucidate the role of viscosity here okay as the one which is actually causing the fingering I want to give a physical explanation all these equations are good but at the end of the day you need to have a physical explanation so let me do that and see this is my flow okay and this is my flat interface okay so I am trying to give a physical explanation you give a small perturbation and what is the story here? Mu low and let us say mu high okay and this is the direction of the flow supposing you give a small perturbation because of which the interface gets deflected the question we are asking is is this deflection going to increase or is it going to decrease and become flat okay so suppose there is an interface deflection will this amplify or reduce that is the question if it amplifies it is unstable if it reduces it is stable let us quickly see what is going to happen supposing this is the deflection remember I am having the same pressure drop at these 2 ends okay the pressure here is uniform pressure here that is uniform now look at this situation here if you look at this small section the effective viscosity in the middle section is going to be lower than the effective viscosity here effective viscosity means you can take something like a weighted average here this let us say 50 50 and let us say this is 60 of this liquid 60 of the lower viscous liquid and 40% of the high viscous liquid clearly the effective viscosity here is lower than the effective viscosity here so you have the same pressure gradient if the effective viscosity is lower this guy will have a tendency to move faster okay this guy will so this guy will keep getting amplified okay the point is mu effective is lower in the thin band okay and hence it moves faster and the thing gets amplified okay but in the reverse case I mean you can do this of course suppose you give a small perturbation of this kind and this is mu high and this is mu low what do you see in the thin band mu effective is high so the velocity is going to be lower and this guy catches up okay is flat so my point I am trying to make here is I mean one just like you have competition between different forces and stuff like that you have competition between different effects so you have to do the mathematics you have to also try to understand the result physically whether it makes sense okay so that is basically the moral of the story here. Now regarding this phi we will have to see I will have to come up with an explanation tomorrow that is why it is happening I just want to finish this up so tomorrow we will discuss why v z bar has to be there but clearly you see if v is 0 then my condition is not going to be dependent on so I knew for sure we had to have the weed but I need to come up with the explanation now.