 What we will do in today's lecture is continue our discussion of the breakup of the jet due to surface tension and these are the equations which we had basically derived in the last class and these are the linearized equations. So we wrote down the equation of continuity and the momentum equations. We had made an assumption of the jet being inviscid, okay, no viscosity because viscosity is not the one which is causing the breakup. What is causing the breakup is the surface tension and we need to retain that effect and that effect is coming from the normal stress boundary condition. So gamma here represents the surface tension, is the surface tension. So maybe I have actually written this as a nu. So let me write this as gamma. So that is the gamma, okay. So this represents surface tension and what I have done is all these equations are with dimensions. So in addition to all of this, we need to use the kinematic boundary condition, okay. The kinematic boundary condition we have to derive for this particular problem. How do you derive the kinematic boundary condition? You write the interface f as r-a times 1 plus epsilon f of z equal to 0, okay. And this f of z, f is of course going to be a function of time as well, f is going to be a function of time as well and what we want to do is write df by dt equals 0. That is your kinematic boundary condition, okay. df by dt equals 0 implies the partial derivative of f with respect to time plus v dot del f equals 0, okay. f is a scalar, remember f is a scalar, okay. So what is df dt? It is just minus a epsilon d small f dt because these are independent variables now, okay. And I have minus a epsilon df dt plus v is going to be, this is going to be u r df by dr plus uz df by dz equal to 0. That is your v dot del f term because there are 2 velocity components, u r and uz, assuming theta symmetry in this problem, okay. So there is no theta component. What is df dr? df dr is just 1. So this gives me minus a epsilon df dt plus 1 plus uz, sorry, df dr is 1 so this gives me u r, okay. Plus uz times df dz is the partial derivative of this with respect to z which is, I am going to have to multiply this by minus a epsilon df dz equal to 0 or in other words minus a epsilon df dt plus u r plus uz, sorry, minus uz a epsilon df dz equal 0, okay. Now this is the kinematic boundary condition with the actual variables u r and u z, the actual velocities. So this is, I have not made any assumption here. I am just saying that the surface is of this kind, okay. And what we have to do now is invoke, see I am interested in the perturbation variables. I have to write this in terms of the perturbation variables. The base state velocities are 0. So this I am going to write as epsilon u r tilde and this I am going to write as epsilon uz tilde, right. So u r remember is u r steady state plus epsilon u r tilde and uz is uz steady state plus epsilon uz tilde. So these guys are 0 and when I substitute this here what do I get? Minus a epsilon df dt plus u r is nothing but plus epsilon u r tilde and this is going to be minus a uz tilde epsilon squared df dz equal to 0. This being a higher order term that is 0. I mean that is not 0, this is a higher order term. So I am not going to consider this. This is of order epsilon squared. So that vanishes and what this gives me is u r tilde equals df dt multiplied by a, okay. So that is my kinematic boundary condition. So if this is the term, this is the equation at order epsilon, what I like to do is I like to add this to my set of equations here which is basically trying to tell you that u r tilde is equal to a times df dt. Remember that we have written this deflection of the interface f is dimensionless because r has dimensions of length and a is here. So f is dimensionless, okay. So just to briefly go through, we had derived the normal stress boundary condition at order epsilon in the last class and these are the linearized equations of momentum continuity. Boundary condition we have just derived here. What I did is just use the fact that kinematic boundary condition comes from the material derivative of f, proceed it further, okay. And I have gotten this thing at order epsilon, okay. We now need to solve this but before I solve this as it is, what I am going to do is I am going to make this dimensionless, okay. So let us make it dimensionless and then see what is to be done. So to make this dimensionless, I need characteristic scales for pressure, velocity, length and time, okay. So what is the characteristic length scale? We are talking about a circular jet which is infinitely long which has a radius of a, okay. So the characteristic length scale is going to be a, lc I am going to choose as a, the radius of the jet, radius of unperturbed jet. The other thing which I am going to do is I am going to choose for pressure, the characteristic pressure. The pressure difference remember at base state is given by gamma divided by r or gamma divided by a, okay. So if I have a cylindrical jet of radius a, I can choose my characteristic pressure as gamma divided by a because this is, my jet is known, this is constant, that is also constant. My characteristic pressure becomes something which is fixed, okay. So once the characteristic pressure is known, I can calculate what my characteristic velocity scale is because pressure goes as rho u squared, right. So my u is going to be, u characteristic is going to be given by square root of pressure characteristic divided by rho which is square root of gamma by a rho. Why is this? Because remember what we are talking about is a stationary jet. The actual problem, the base problem does not have any velocity. I am looking at a thread which is actually stationary. So there is no velocity which is characterizing the flow. So I do not have a characteristic velocity from that but what I am saying is whatever is the velocity is going to be induced by the surface tension which is breaking up. So that is the reason the characteristic velocity has the surface tension occurring in the definition, okay. So now I have my characteristic length scale and velocity scale. So my time scale is easy to be found out. My characteristic time scale is going to be, Tc is going to be given by, lc divided by uc, okay. This has units of time and so I would have a divided by this, okay. This is right. Yeah, so we proceed. This power 3 by 2 looks funny, but so with this I want to make the equations dimensionless, okay. And what we will do is we will just go through with the process here. Let us take the first equation. The first equation is rho times du r tilde by dt equals minus dp tilde by dr. This is with dimensions, okay. So if I want to make a dimensionless, I have to take out my characteristic velocity out of this. When I take out the characteristic velocity, I am going to take out square root of gamma by a rho. And I have to find some other symbol. So I am going to call it du r star by dt where ur star is defined as u r divided by u r characteristic, yeah, u r characteristic, okay. u r characteristic is the same as uz characteristic and both of them are, so I am taking out this factor. What about t? Okay, let us do t later. Minus, when I take out my pressure characteristic out, I am going to get gamma by a dp star. And when I take out, make this, I am just converted this to dimensionless. I am going to make the independent variables now dimensionless and I will get rho times square root of gamma by a rho times take out a t characteristic which is here, which is a to the power 3 by 2 times square root of rho by gamma, okay, du r star by dt star. So what I am doing now is making the independent variable dimensionless. I mean you guys can do it faster and minus gamma by a square dp star by du r star. The point I want to make here is that, okay, wonderful, everything is fine. So I have just went from dimensionless, so the stars represent the dimensionless variable, okay, without dimensions. So this is without dimensions. So you can see what happens now is this gives me, this rho cancels with this square root square root, this gives me a square, this cancels with that and this gives me gamma that cancels and I basically got rid of my coefficients which were hanging around, okay. So that is my dimensionless equation. You can do the same thing for all the variables and we can get our dimensionless equations now, okay. I am not going to do it for the other variables but I just want you to know clearly this is unique because that is the only lens scale. This is coming because the pressure difference is given by the surface tension force and I am using that and once this pressure is defined, I can use this to get uc and tc. What we will have when we make a dimensionless is dUr star by dt equals –dp star by dr, okay. dUz star by dt equals –dp star by dz. Come to the next one, I should possibly put a star here, something is wrong, 0 and when it comes to the kinematic boundary condition, you will get something similar, ur star equals star, okay. And this guy here P1, I am going to scale with gamma by A, so when that comes out, P1 star is going to be given by –f double prime. So, you can do that and you can see. So, these equations are my dimensionless equations, remember f is already dimensionless. Yes, I mean, we can take a different but what is the characteristic lens scale in the z direction. The wavelength is what we are going to find out. The wavelength of the disturbance which is most critical which is going to decide the breakup of the drop is what we are going to find out. We do not know what the wavelength is. If I had a characteristic dimensions in the z direction, I could possibly use that. If I use that, then that particular lens scale would come in the differential equation, otherwise it will come in the boundary condition. So, basically if there are 2 lens scales, depending on how you define your characteristic variables, they will come either in the boundary condition or in the differential equation. So, that parameter will appear but either where does it appear, that is the only thing which is going to be different, okay. So, we need to solve these equations and what I am going to do is I am going to drop the stars from now on, okay. So, let us drop the stars but the variables are dimensionless. Remember that, okay. It is just for me to make life easy. Otherwise, I keep forgetting the stars somewhere. Wonderful. So, I am just going to follow what is done in Gary Leel so that it is easy for you to refer. So, this equation can be solved in many ways but remember this f is a function of z and t, okay. And our objective is to find the relationship between the growth rate and the wave length or the wave number, okay. So, f is a function of z and t and I am going to seek f as e power sigma t sin kz. I have got to use exponential but I am just going to use sinusoidal. You can use cosine. It does not matter. What am I doing? I am looking at a periodic disturbance in the z direction, okay. So, this is infinitely long. I am just saying I have a disturbance in the z direction which is infinite which is periodic. This is going exponentially in time because my equations are linear and this is representing the amplitude in some sense. The actual disturbances are going to be arbitrary, okay. And these arbitrary disturbances I am going to be able to decompose them into different Fourier modes. And that is the justification for seeking this periodic function in the z direction. So, what I will do is I am going to find out for different k values which is the one which is going to grow. If it turns out that for all values of k, sigma is negative that means it is not going to grow. It is stable, okay. If it turns out that for some k it is sigma is positive that means it is unstable. So, any arbitrary disturbance is going to be decomposed into a bunch of Fourier modes. We find out which of this is going to grow, which of this is going to be unstable. And that is what we always are going to be doing in this course, okay. And whenever you do linear stability analysis of infinite systems which are extending to infinity in some direction, this approach is used. If you have a finite system then you will use a finite Fourier transform in terms of sin n pi x by l sin 2 pi x by l. But since it is infinite we use the actual Fourier transform, okay. So, the other thing is I am going to jump directly to my normal stress boundary condition which we derived. So, I can substitute this thing for f here and I can find out what is p. So, p1 turns out to be –f –ce power sigma t sin kz – of f double prime which is – this is remember the derivative is with respect to z. So, I am going to get a k square and that is going to be the plus sign and that is going to be given by ce power sigma t sin kz times k square – 1. That is my pressure, okay. What I am saying is if this is the form of the interface, the form of the pressure is going to be given by this. But remember this is at, I have got this from my normal stress boundary condition. So, this is at r equal to 1, okay. This is since this is from the normal stress boundary condition. Now, I know what the value of the pressure is at the boundary. So, what I need to do is I need to get the differential equation for pressure, okay. Remember what are these 2 equations here? I can combine these 2 equations from the Navier-Stokes equations and I can write this as du by dt in a vectorial form as – gradient of p. What have I done here? I am just saying that I can look at these 2 equations, write this in a vectorial form. The r component is going to be du r by dt is – dp by dr. The z component is going to be du z by dt equals – dp by dz, okay. And remember this is nothing but the divergence of u equals 0. So, I am going to take the divergence of this equation. And when I do that, I will get divergence of this scalar, I mean this is a time derivative operator, the divergence is spatial. I can move the divergence operator inside here and I will get d by dt of divergence of u equals – divergence of del p. The left hand side is going to be 0, the right hand side is going to be del square p. So, basically I am going to get del square p equals 0, okay. So, taking divergence we get del square p equal to 0, p1 for whatever reason I put a subscript 1 here, okay. So, that is my differential equation. I need to solve this differential equation. What are the independent variables z and r? So, I mean this is basically going to be 1 by r d by dr of r dp1 by dr plus d square p1 by dz square equal to 0. And since it is a second order equation in z, infinitely long and it is homogeneous, you can seek the solution in the form of some variable separation. You will get some trigonometric function in the z direction and you will get a Bessel function in the r direction, okay. So, basically what I am saying is this solution p1 is going to be of the form a of t cosine kz plus b of t sine kz times i0 of kr, i0 of, yeah. So, I mean I am not going to be doing the math here but you guys can go and check if this is indeed right. Now, the 2 things, this is a second order in z, second order in r and you need to have, you will get 2 solutions, right. It will be i0 of kr and there will be another solution which is k0 of kr. You would actually get a B, another say c and d, 2 constants associated in the r direction. So, these are your independent solutions in the z direction sine and cosine and since it is a variable separation form it is going to be i0 and k0. Now the fact that k0 is unbounded, basically it is going to knock off this contribution. This is knocked off because k0 is unbounded and r equal to 0 and what I am left with is only this. Only i0 is going to contribute. C multiplied by b is some other arbitrary constant. C multiplied by a is some other arbitrary constant. The thing which I want you to remember is the differential equation here has only r and z. So, why am I putting a as a function of time and b as a function of time? The reason is the boundary condition has the time dependency because on the boundary the variable is changing with time. The pressure inside is also going to change with time periodically or in whatever way it is, okay. It may not be periodic, it may be exponential. So, the reason why if this pressure had been independent of time then a and b would have been just constants but because the pressure here is actually changing with time, these guys are not constants but these are functions of time, okay. So what I am going to do is I am going to write this P1 as a multiplied by c is some other constant with let us say e of t times cosine kz plus f of t sin kz times i0 of kr, okay. Now what I want to do is I want to compare this guy has to collapse to this value or this function at r equal to 1 and r equal to 1 this must match this and at r equal to 1 I only have the sin dependency I do not have the cosine dependency. So, what this means is e has to be 0, okay. So, from the boundary condition at r equal to 1 P1 equals f of t sin kz times i0 of kr, okay. So, this periodicity in the z direction is the same as what is being imposed by the boundary condition. The variation in the radial direction has come by the governing differential equation here and the amplitude f of t is something which we need to find out, okay that is still an unknown quantity but what I have done is I basically got in this, okay wonderful. So now I need to get a relationship between f and c and how can I do that? This is the pressure is given I need to use one of these conditions here I do not know velocity yet I know pressure I am just going to put at r equal to 1 this pressure must be equal to this, correct and that will tell me how f of t is related because then I can equate the 2. So, I am going to use this boundary condition at r equal to 1 my pressure is c e power sigma t sin kz times k square-1 equals f of t sin kz i0 of kr am I missing something? i0 of k you are right i0 of k sin kz goes off. So, f of t is nothing but c times e power sigma t times k square-1 divided by i0 of k, okay. So, f is known now I have been able to relate so this is what we are trying to do is these arbitrary constants which come I am going to try and eliminate them and try to get in terms of one constant and then we are going to finally say that that constant cannot be 0 and then get a relationship between sigma and k square, okay. Otherwise you can keep all those constants write a determinant and say that the determinant should be 0 in order to get a nonzero solution. So, either approach is fine but we are going to do it this way. If f is known then what is my pressure? My pressure is going to be given by I am substituting this f back there c e power sigma t times k square-1 sin kz times i0 of kr divided by i0 of k. I am going to use my radial component of velocity balance, okay to find out u r star. So, the plan is this substitute for differentiated with this way to r find out dp by dr I will get d ur by dt and then from this I get u r. Once I know u r I can substitute the value of u r here I already know what f is I have assumed it to be of the form that find out what this is because the kinematic boundary condition is yet to be used. So, idea is I am going to get u r in terms of c I am going to get f in terms of c and everything will be fine, okay. So, let us find out what is dp by dr is going to be c e power sigma t, differentiating i0 gives you k times i1, okay and this gives me k times k square-1 times sin kz times i1 of kr divided by i0 of k that is dp by dr, okay. So, c comes here this is what I get and this remember is d ur by dt is a negative of this quantity here. So, to find out u r I am just going to have to integrate this with respect to time, okay I am going to integrate this with respect to time and I get u r equals when I integrate with respect to time I get the sigma in the denominator, right. So, c times e power sigma t divided by sigma k times k square-1 times sin kz times i1 of kr divided by i0 of k. So, all I am doing is integrating this with respect to time and I get e power sigma t divided by sigma, okay, beta-s, beta-s sign because there is a minus here d ur by dt is minus dp by dr, wonderful. They are pretty much done and this remember at r equal to 1 must be at r equal to 1 u r star equals minus df by dt, okay. This is going to be equal to r equals minus df by dt, minus df by dt or plus df by dt, it is plus, right. Why not write minus here, is it plus, okay. So, what is plus df by dt? Remember f we have already assumed to be of the form c multiplied by e power sigma t. So, df by dt will be sigma e power sigma t times sin kz. So, ur will be equal to this I have to evaluate this ur at r equal to 1, okay and I just check, it is plus or minus, it is plus, okay. So, ur is equal to c sigma e power sigma t sin kz at r equal to 1 and I just go to this equation here and I get c minus c e power sigma t by sigma times k times k square minus 1 times sin kz, okay. So, right now what I will do is I will just erase this plus d and I will justify to you as to why we are neglecting this d, okay. This must be equal to c sigma e power sigma t sin kz from the kinematic boundary condition. So, e power sigma t sin kz cancels and c has to be nonzero, remember, okay. So, what do I get? c has to be nonzero implies sigma square equals 1 minus k square multiplied by k times I1 of k divided by I0 of k, a nonzero k for nonzero c. If you want a nonzero disturbance, basically we are looking for a set of conditions when your linear equation has a nonzero solution, linear homogenous equations have a nonzero solution and if this condition is satisfied, you have a nonzero solution, okay. So, what this tells you is what is the growth rate for different case? So, basically this answers the question for each wave number k, what is going to be the corresponding growth rate, okay. So, this tells you and if you want to plot this function on the right as a function of k, you would be able to get the dependency of sigma square on this. So, let us do that. So, we plot sigma square at k equal to 0, k equal to 1, this thing is 0, that is k equal to 1. In between, k equal to 0 and k equal to 1, this k is positive, okay. If k is greater than 1, sigma square is negative. If k lies between 0 and 1, sigma square is positive. So, for 0 less than k, less than 1, sigma square is positive, which means what we are really interested in is sigma, okay and that means sigma will have be plus or minus square root of a positive quantity. So, we will have a positive value and a negative value, okay. So, this means sigma is positive for this range and again for k greater than 1, sigma square is negative and actual thing is it is neutrally stable because the real part is 0. What this means is if you are going to see and then you also see that the maximum growth rate is going to occur at some point in between. Some wave number in between is the one which is going to grow fastest, okay. So, when you give an arbitrary disturbance, it is going to be made up of different wave numbers or different wavelengths. The wave number which is going to grow fastest is the one which is going to dominate and that is the going to be giving you the indication for what is the break up length for example. Because k remember is wave number which is reciprocal of wavelength. So, this particular thing you can calculate and if I remember right, this is about 0.6, okay. So, what this means is the maximum growth rate occurs for sigma for k equals 0.6 approximately, this is right, 0.6 Jason, you think 0.6 is right, 0.6979, okay, 0.697. So, there is almost 0.7 then. Now, what does this mean? The wavelength then you could observe the lambda critical when surface tension is going to actually pinch and going to break this thread just stationary. The lambda critical is going to be given by 2 pi divided by k, okay and it is going to be given by 2 pi divided by k which he tells me is 0.697. So, I will just go with that 0.7 but remember this all being done in dimensionless. So, actual length is going to be a times that, okay. So, it is going to be 2 pi actual wavelength is going to be 2 pi divided by 0.7 times a. I think the jet is unstable when the wavelengths larger than 2 pi a because k goes from 0 to 1, k is less than 1, lambda must be greater than 2 pi by k. So, that is basically this tells you that the wavelengths which are actually unstable are larger than 2 pi by a. Any disturbance which is having a wavelength which is lower than 2 pi by a is going to be stable, okay. But that we cannot really conclude from this. We can only conclude about the unstable portion. What is I want to say? What we will do in the next class is try to get this upper bound using another method, okay. This range of wavelengths where you can decide the stability the threshold value. We are going to derive this using what is called the energy work principle, okay to derive the critical condition on wavelengths. That argument is slightly different from that approach is slightly different from what we have done now. In the sense that is more of a static argument. We do not use the dynamics, okay. So, idea is that what we have done the linear stability analysis is we are beginning with the actual governing equations and we are getting the condition of our stability. We are finding what the growth rate of the disturbances, okay. The time dependent factor is actually captured in the linear stability analysis. In this approach, the time dependency will not be captured. But you are going to use some energy argument and you are going to find out your critical wavelength for stable and stable behavior. And then we will compare these two approaches, okay. That is the idea. So, we will answer this question about this integration constant being 0.