 This is probably the most important and the most difficult topic, so as a teacher there is a challenge. So I tried my best and it is up to you to see how far I had been successful. Solution of new stoke equations on a staggered grid. Let me go to the background solving heat conduction equation or heat advection or convection without fluid flow. All these are much easier as compared to when you go to solve the new stoke equation. So that is why this topic is much more challenging or much more difficult as compared to that topics. So there are a lot of details, concepts in this. So I will start with introduction. Second topic is finite volume method. So you may be thinking why finite volume method so many times. I did it for conduction, advection, convection, then for full name stoke equation. Then I am again talking of finite volume method. The reason being that on a staggered grid you will see that we use certain mathematical procedure when we use finite volume method. Then we will go to the formulation that is semi-explicit and semi-implicit method and finally I will discuss in lot of detail this semi-explicit method. Now let us go to the topic on solution of new stoke equations on a staggered grid. So in the end of this topic I will discuss the implementation detail for the semi-explicit method and solution algorithm. Actually we have given a code which is made up of the semi-explicit method. Note that in all these lectures I am basically mostly talking of the implementation detail for semi-explicit method because we believe that for a beginner it is easy to program in a semi-explicit or fully explicit way. So these are the governing equations which I had shown earlier. This is the general transport equations. This is the way we do the grid generation, the simplest form of the grid generation. Note that this is the first course. That is why we are starting with the simplest form of the grid generation. Although when you look into software when you are solving a practical problem you may have much more complex domain shape and size, the grid size. But for a beginner we believe that this is a good grid generation to begin with. Now let me discuss the issue or the difficulties in new stoke solver. Before that I would like to go back to the governing equations. Now I would like to draw your attention to the nature of the equation. What is your ultimate objective? Our ultimate objective is let us say that in a CFD course we try to develop a video camera like tool. So let us suppose I want to create a movie for u velocity. So when you want to create a movie for u velocity you should have a term del u by del t in some conservation law or some equation. So in this equation you see that you have a del u by del t term. So you can consider this equation as an equation from which you can create a movie for this equation as an equation from which you can create a movie for u velocity. You can consider this equation as an equation where you can create movie for v velocity which is the other variable which is left. If you look into this three equation other than u velocity and v velocity which is the other dependent variable pressure. Now is there any conservation law or is there any governing equations where you have del p by del t whether you have a pressure temporal gradient of pressure because if you want to create a movie you want picture with respect to time. Now if you want it from equation there should be a time there should be an equation or there should be term in equation which should have a temporal derivative. So what you realize is there is no conservation law or there is no governing equations which has del p by del t rate of change of pressure. We have rate of change of x momentum from which we get rate of change of u velocity. We have rate of change of y momentum from which we get rate of change of v velocity. We have rate of change of internal energy or enthalpy from which we get rate of change of temperature but we do not have a conservation law which says rate of change of pressure. This was one of the reason that there was lot of difficulty in solving this set of equations. There is no specific equation for pressure but there is one conservation law which we are not using right now. X momentum we are using for let us say to create movie for u velocity, y momentum to create a movie for v velocity and energy equation to create a movie for temperature. So we are left with one variable pressure and we are left with one conservation law or one equation which we have not used which is the mass conservation. So can I call this equation as an equation from which we will solve pressure. Right now you have 4 equations and right now we are trying to see that from which equation which so the x momentum is for u velocity, y momentum for v velocity, energy equation for temperature. So can we call continuity equation as an equation for pressure but we do not have any pressure in continuity equation. So how we can call this continuity equation as an equation for pressure. So this is one of the questions which needs to be understood to solve the Navier-Stokes equation. So we have to do something and we have to convert this mass conservation equation or continuity equation to a pressure equation that is what is done in CFD. That is why it took lot of time before people started solving full Navier-Stokes equations. There was lot of initial struggle. Initially people were solving the fluid flow by using another method called a stream function vorticity method in which what is done is that this x momentum equation is operated by del y del y. So when you operate this x momentum by del y del y you get del square p by del x y del x del y and the y momentum is operated by del y del x. Then you get del square p by del x del y. So we do cross differentiation and then subtract so that the pressure term is completely animated that is what is done in stream function vorticity formulation. As this pressure is not having any specific equation let us so if you are not able to handle let us eliminate that. So that is what is what gave rise to stream function vorticity formulation and this is what is called as a primitive variable formulation. These variables are called as primitive variables in fluid mechanics because those stream function and vorticity are derived variables. They are not fundamental variables. So after discussing that let me so even let us suppose people converted the continuity equation into pressure equation and then they started solving this. Now I would like to draw your attention to the second level of approximation which we have used in case of pressure term and in case of mass fluxes. So what I am showing you here is the application of second level of discretization which is the algebraic which is the right now when you do the finite volume discretization of the pressure force in the x direction you get p on pressure on the faces on the west face the pressure force in the positive x direction. So you have a plus sign on the east face pressure force is compressive. So it is in the negative x direction so there is a negative sign. Now this pressure force on the faces are calculated using linear interpolations and if we have a grid structure such that the face center is lying exactly in between the cell center then p small e is calculated as p capital E plus p capital P divided by 2. Similarly pw is calculated by this expression surface area for both the faces is delta y. So if you take this difference of the pressure of at the faces of the control volume you get an expression you get an algebraic expression which is in terms of cell center value that is fine but in this expression there is no pressure value of the same control volume for which you are solving. You can see that this p capital P subscript p cancels down. So let us suppose you have a pressure distribution like this. Now this is the algebraic representation which you have got actually this multiplied by delta y. So if I apply this equation here you will see that it is indeed obeying this equation. You have got a number distribution like this. So if I take if I apply this equation at this point then this is my p subscript capital P this is p subscript capital E this is p subscript capital W. So if I apply to this the equation which I have for this control volume or for this great point the equation is p capital W which is 10 minus p capital E which is again 10 and the pressure force is coming up to be 0. When you apply to this control volume this is p capital W and this is p capital E. So again it is coming up to be 0. What I wanted to tell you is that if we follow this type of second level of discretization for pressure we get a convert solution. It is not that the solution is diverging. We get a convert solution. We get a number distribution which is obeying let us say mass conservation momentum conservation but do you expect that the flow field will be like this zigzag with a zigzag pressure field? No. So you get a number distribution which is following the conservation laws but it is not correct. So this gives rise to zigzag pressure field. So the point is we should try to avoid this type of interpolation. So either we should come up with some different interpolation or extrapolation or we should not have this type of interpolation at all. We should have grid points for pressure at the phase center. This to give you an idea. That is the only way you can avoid this procedure completely. This is what is done in staggered grid. Similarly when you look into applying second level of discretization to the normal velocity on the phase center. This is the discrete representation of this is after first level of discretization you express the mass fluxes at the phase center. You end up with normal velocities at the phase center and for this again if you apply u small e is equal to u capital E plus up divided by 2 u small w is equal to up plus u w divided by 2. Again you end up with the continuity equation which is not using the 2 neighboring difference of velocity of the 2 neighbors but which is using the alternate. Like here if I apply this equation this is the east neighbor and this is the west neighbor. So here it is 10 and this difference is 0. Right now I am showing you in one direction you get a wavy field in the y direction also. So you get a number distribution which is wavy in nature. So alternate numbers are changing and there are only 2 numbers but that number distribution obeys this algebraic equation which is the finite volume representation of continuity equation. So this is just to show you that you get a number distribution which is obeying mass conservation but you know that this result is non-physical. So the issue is how to avoid the interpolation for normal velocity on the phase center as shown in this slide and how to avoid pressure interpolations at phase center and the first remedy which was proposed is what is called as staggered grid. I will show you this through an animation. So the what is being proposed is this. This is the control volume v. Let us say this is a lead driven cavity problem. The problem which will be solving in today's lab session the first problem. So what is that you have a square cavity and the top wall of the cavity is like a conveyor belt which is moving in the x direction with a velocity u0 and v is 0. All the other walls of this cavity is stationary. So let us suppose for that problem we obtain the grid generation by drawing equispaced horizontal and vertical lines and we get certain interior control volumes. Then I would like to point out that this yellow circle which is shown here in a staggered grid we will solve for pressure and temperature at this grid point and this green square which is shown here which I had shown earlier for calculating heat fluxes in the x direction or enthalpy flux in the x direction which is located at the phase center. These are the grid points we will use to solve for u velocity. These are the grid points which we will use for solving for v velocities. These are the boundary grid points. At boundary we will have all the three different types of grid points because we have to apply boundary condition for let us say u velocity v velocity temperature pressure. So this is your staggered grid. I will show you in the next slide this is the running indices. These are the boundary nodes. So in this slide I will be showing you the different control volume. So this is the control volume for pressure and temperature which is called as a scalar control volume. This is the control volume for u velocity. Now note that in this u velocity control volume on the phases there is a circle and that circle is for pressure. On which phases is this circle? On vertical phases and which are the pressures which acts to calculate pressure force in the x direction? We take pressure on which phases? We take pressure on vertical phases because pressure is always normal in this compressive. So note that in this u control volume if you want to calculate pressure force in the x direction then the pressure grid point is lying at the phase center. So if the value, if the grid point is lying where we want it so what there is no need for any interpolation or extrapolation. This is the way we avoid pressure interpolations at the phase center of u control volume. I am showing you here this is the pressure control volume, this is the u control volume and this is the v control volume. Now the way it has been done that this u velocity u control volume I am taking the grid point at the east phase center positive phase center. The v velocity control volume we are taking on the positive phase and phase center again. This is the positive phase of that vertical phase. This is the positive phase or north phase or horizontal phase and these are the three different control volumes and I am also showing you the grid points which are lying. Let us see this control volume. In this control volume note that this green square which has u velocity. This control volume as I said that in this control volume we used to apply law of conservation of mass. This control volume we applied to we used to apply law of conservation of x momentum. This control volume we used to apply law of conservation of y momentum because in this control volume pressure is sitting at the center and we are saying that pressure variable we will obtain from mass conservation. This I discussed in my previous slides. You can see that in this control volume u velocity is at the centroid so we will apply x momentum here. Similarly v is sitting at the centroid of this control volume so we will apply y momentum here. When you are applying mass conservation laws in this control volume then what is the flux which comes mass flux. Now on the vertical phases which velocity constitute the normal velocity u and you see that the grid point over u velocity is sitting on the vertical phases. When you go to the horizontal phases the normal velocity is v and here again v velocity grid point is sitting. Interpolation interpolations are needed only if at the phases the grid points are not there but here the way we have done staggering we are making sure that the interpolation of the normal velocity to calculate the mass flux is not needed. Once you ensure this then we say that the coupling pressure velocity coupling is maintained otherwise we say that there is a pressure velocity decoupling. What is you see this word in many of the CFD books what effectively it means is this pressure velocity decoupling is that if you do interpolations to calculate the velocity or pressure of a particular grid point the algebraic equations which you get as sure I will show go back to the earlier slide like to calculate pressure at this grid point the equation which you are using the expression is not having pressure of the same grid point that is this. There is it is said that there is a breakage in communication between the velocity of the same point and the pressure of the same point because in the next slide when you use linear interpolation for the normal velocity you end up with equation where the velocity of the same grid point is not used. So, in your calculation procedure mass conservation momentum conservation to calculate velocity and pressure of a particular grid point which is notated by p your expressions are such that it does not have subscript p. So, to calculate velocity and pressure of a same grid point your discretized algebraic equation does not have expression for velocities and pressures of the same grid point. So, there is a breakage in communication between the velocity and pressure of the same grid point and this is what is called as pressure velocity decoupling both this linear interpolation for normal velocity as well as pressure as you can convince yourself through this figure. So, here I am showing you the different grid points also I would like to draw your attention that how I had decided that u is this is u capital P u capital W how this is v capital P this is v capital S here again we are following a convention which I had discussed earlier. Our convention is that on the positive faces we have this is p represent i comma j. So, we are using a convention where you can see that on the positive faces is having same i j which is there for the centroid and this is i minus 1 and this is j minus 1. So, we are using the same thing which I had discussed earlier in case of conduction advection and convection. So, and you can see that to calculate the when we do the mass conservation then to calculate the mass flux here you take rho up into delta y to calculate the mass flux here or mass flow rate rho uw into delta y rho vs into delta x rho vp into delta x. So, this way you calculate the mass flow rates and note that you are not doing any interpolations. So, you have avoided the interpolations completely in the mass conservation and when you look into this x momentum applications now you see that the pressure you when you want to calculate pressure force in the x direction. Here we have eliminated the interpolation for the normal velocity here we have eliminated the interpolation of the pressure at the face center. Similarly, to calculate the pressure force in the y direction you can see that we have eliminated the interpolations which are needed for pressure. So, with this staggering of the grid points we have avoided two type of interpolations interpolation first interpolation for normal velocity second interpolation for pressure I will show this through an animation here this is a grid point for pressure this is for u velocity this is for v velocity. So, let us take this control volume 1 by 1 this is a control volume where we apply a law of conservation of mass. So, when we apply law of conservation of mass this is the velocity which is a u velocity which is the normal velocity on the vertical face is used to calculate the mass flow rate on the this face and the grid point is lying there. So, we have placed our grid point such that we do not have to do any interpolation now the next thing we have put the v velocities you can see in such a way that again no interpolations are needed all the grid points are placed appropriately. And the conventions which we are following based on that convention you can see that this is u p, u w, v p and v s. Now, let us take a u control volume now when you take a this as u control volume this is u p can you tell me where will be v p sorry pressure grid point you can see in the slide that if this is u p I am asking you that in this control volume where will be p p. So, you can see this and very easily answer that p p should be here. So, it should be here that is what I will show now and where should be p e you know that this is p p then here this there should be p e that is what is shown here. But, before that I would like to remind you but, before that I will show you if p p is this then v p is here. So, let me show you v p this is v p then this will become v s this is p then this becomes p e and this becomes v s e I would point out subscript p is i comma j subscript e is i plus 1 comma j subscript s is i comma j minus 1 subscript s e is i plus 1 j minus 1. So, in this control volume we apply law of conservation of x momentum in this control volume we apply law of conservation of y momentum. Now, where will be p p here now p p should be at the bottom if this is v p p p should be at the bottom of this control volume which is shown here. If this is p p you know on the east side you have u p here which is shown here. If this is u p this will be u w this will be if this is p p then this should be north neighbor and here again if this is n the east face of that this will be u n and this will be u n w n w means i minus 1 j plus 1. The reason I am showing all this grid point involved in a particular control volume there are three types of control volume and there are in each control volume I am showing all the different grid points which comes into those control volume. The reason being that when you go into the formulations you will realize that so note that although you have avoided the interpolation to calculate the mass flow rate in this control volume but note that mass flow rate we need to mass flow rate actually driver in the adduction equation. Mass flow rate is used for mass balance in mass conservation but it is also used in x momentum equation it is also used in y momentum equation. Now, in this control volume when you want to calculate the mass flow rate as you have done staggering now you see that to calculate mass flow rate you need u velocity here u velocity here v velocity here and v velocity here. Now, note that the there are three mass flow I mean say as far as special locations are concerned the special location to calculate the mass flow rate here is different the special location to calculate the mass flow rate for applying law of conservation of momentum is different special location not to calculate the mass flow rate here what we do is that we let us say I need v velocity here which are the nearest v velocity here v capital E and v capital P. So, I do a linear interpolation and I calculate v velocity on the north face center as v capital E plus v capital P divided by 2 which I will show you in the next slide. Here velocity I calculate as u capital P plus u capital W divided by 2. So, that way I want to tell you that mass flow rate is not only used in the mass conservation it is used in the advection term in the momentum equation and to calculate the mass flow rate we avoid any interpolation in the mass conservation but we have to do interpolation to calculate mass flow rate for momentum conservation which I will show you in detail in the next slides. So, this is a slide to show you that we are avoiding any interpolation we have when we apply law of conservation of mass. So, here we are applying only first level of approximation we are not applying second level of approximation because second level of approximation is not needed because the velocity grid points are staggered such that they are lying at the face center of the pressure control volume which is the control volume where we apply law of conservation of mass. Now, this is the way you can discretize the unsteady term this is this involves only the grid point of the cell center. So, this is done for rate of change of x momentum in this control volume rate of change of y momentum in this control volume and rate of change of temporary internal energy in this control note that I had told you that the yellow circle is control volume for scalars it is not only the control volume for where we apply the mass conservation, this is not the control volume where we apply mass conservation but in this is a scalar control volume, temperature is also a scalar so we apply energy conservation also. In case if you apply energy conservation note that in this case if you apply energy conservation here to calculate the enthalpy flux now when you are calculating the LP flux for this control volume the mass flux which are involved are such that for that you will not need the any interpolation. So, there are three special location of mass flux first in three special location of mass flux in three different types of control volume. For the law of conservation of mass as well as for the law of conservation of energy the mass flow rate which is involved we avoid any interpolations. But for law of conservation of x momentum and bi-momentum we have to do interpolation which I will show the expressions in the coming slides like shown here. You can see that the expression is note that I have written m subscript there are two subscript first subscript representing the phase and the second subscript represent the control volume for which this is the mass flow rate. Sorry this is small m is not mass flow rate it is mass flow rate per unit area this is mass flux. Here the symbol small letter is used for flux term small m for mass flux small a for advection flux small d for diffusion flux. So this is mass flux which is the product of density multiplied by normal velocity. So for east and west phase which are vertical phases the normal velocity is u. Now here we do not have green squares and here we do not have red inverted triangle. So we do not have great point for the normal velocity. So to calculate the mass flux we need u velocity here and which are the nearest neighbor u p and u e. So we take u p plus u e divided by 2 which is shown here. Similarly here it will be u p plus u capital W divided by 2 which is shown here. When you go to the north phase we need v velocity and which are the nearest v velocity grid point v capital E plus v capital P divided by 2 which is shown here. Similarly we do on the south phase. We do the same thing in the y momentum equation. Note that here we have taken in the for east and west phases we have taken the u neighbors in the horizontal direction. But in the next slide you will see to calculate u here you take the vertical neighbors because those are the nearest neighbors. To calculate u here you take u n plus u p divided by 2 which is shown here. To calculate here you take u w n plus u w divided by 2 and here now v is calculated using the vertical neighbors v p plus v capital N divided by 2. I started with unsteady term now I am actually first I started with the mass conservation where I said there is no interpolation involved. Then I went into the unsteady term and then I am going one by one to the different advection term x momentum. Now this is y momentum. Next slide will be energy equation advection term which is enthalpy or internal energy lost by the control volume. Now here when you want to calculate the total advection the normal velocity grid points are lying. So for that you do not need any interpolation. So in the previous two slides this slide and this slide we were doing interpolation for the normal velocity. For the energy equation it is not needed. Most of the CFD books which you see they show you all this expression. But what is different in this lecture or slide is that I had given a pictorial information because you as a teacher know that when you want to specially you have to develop a program using all this formulation. So whether you understand the formulation only through equations or equations supplemented by figures which is also supplemented by animations then the memory retention or understanding improves a lot that is what I had tried. And after advection next thing is diffusion. So this slide I am showing you now I have to calculate viscous forces in the x direction. Now this viscous forces for viscous forces you calculate normal gradient of u velocity. Only u velocities are involved in this calculation. So in this slide I am showing you only u control volume and its neighbor and the approximation second level approximation is very much similar to what I had shown you in the previous slides for conduction. Like here viscous stress is mu del u by del x. So it is u capital E minus u capital P divided by the distance between the two points which is this. Similarly on the other phases this is for this phase, this is for south phase, this is for west phase and this is for north phase. This is to calculate the total viscous force in the y direction. This is to calculate the total viscous force in the y direction and this is to calculate the total conduction heat fluxes or total heat gain by conduction. Note that here I am using plus k again I would like to point out and I had reversed the direction of the arrow because in our summation I would like to remind you that we do east north minus west south. Now when you have east north now it becomes inflow when it is a plus k. So do not get confused at how the direction has been changed. Note that here it is a plus k. Let us go to this window maybe I will show you some animations just to give you a better feel of what has been taught just now. So this is the application of x momentum equation in the u control volume. This is the way x momentum is entering and leaving the control volume. Note that the mass fluxes which are involved at the phases we do not have normal velocities. So for that we have to do this interpolations. Similarly when you go to apply momentum conservation to calculate the advection this is the advection momentum flux in the y direction entering and leaving from the vertical phases this is from the horizontal phases. And as again mass flux are not the normal velocity grid points are not at the center of the at the phase center we have to do this interpolations. Similarly this is for the energy equation. Here note that we do not have to do any interpolation. So here I had shown you the viscous forces also and in the viscous forces we do piece by linear approximation which we had already done for conduction heat fluxes. And this is the pressure force acting in the x direction. Now here note that we do not need any interpolation so that is why I am not showing you any expression here because pressure is lying at the phase center directly. Here again pressure when you apply want to calculate pressure force in the y direction pressure grid points are lying at the centroid. So this is not needed. Now let me show this. So this way I had completed the finite volume discretization in staggered grid. So earlier topic was also finite volume method but I had completed the first sub topic of finite volume method. So there was a repetition but there was a reason where I wanted to show you that in this staggered grid what are the important points which you have to remember where you have to do interpolation where we have avoided the interpolations. So let me show you this through an animation which will be coming into this window. So this is the for staggered grid rate of change of x momentum of the fluid inside the control volume. This is the rate at which x momentum is changes across the control volume. This is the viscous forces acting in the x direction in a staggered grid. This is the pressure force acting in the x direction. Similarly this is for the y momentum equation in a staggered grid and similarly this is for the energy equation. So we had completed the finite volume discretization but before we break there is one idea which is probably I would suggest the most important or most difficult concept to understand. I would like to end this lecture today through that concept and I had tried hard to create animations and to discuss philosophy rather than expression to convey that idea. So what I am trying basically to tell you is let me come back to that animation directly because we do not have much time. As I said that we mass conservation we convert into an equation for pressure. Now there is a mathematical procedure where we do it but in this slide I will show you discuss the philosophy of that. So let us come to this animation. So let us say there is a control volume. Actually the idea is something like this that if you follow a fully explicit method in case of solving the Navier-Stokes equations then if you follow a fully explicit method then if you are solving x momentum equation then you calculate total advection, total viscous forces acting in the x direction total pressure forces acting in the x direction and you get a u velocity. Then you solve a y momentum equation and when you solve it by a fully explicit method from there you get a v velocity. Now this u velocity v velocity which you have got by a fully explicit method if you calculate the mass fluxes you get this type of mass fluxes. So now the problem is is this physical whether it can happen for an incompressible flow in a real world situation the answer is no. So this is a numerical situation and this is what we call as a first iteration of when you want to obtain pressure. Pressure you do not know. So you have to start with by taking pressure of the old time level but when you do that you end up with this number distribution and this is not obeying mass balance. So we have to come up with a procedure so this picture or animation is for first iteration this is not a real world situation this is just a numerical situation. This is the first iteration of after you have predicted velocities you get this type of mass distribution. Now by looking into this picture and if I say that you have a control that you can control the pressure of this control volume then by looking into this picture you realize that 4 kg per second is going out and 2 kg per second is going in. So more mass is going out and let us suppose you have a radio type knob where by which you can control this pressure of this control volume. So when more mass is going out whether you would increase the pressure or you would decrease the pressure. When more mass is going out you would decrease the pressure. When you decrease the pressure more mass will come in so that you hope that mass balance will occur. So this is shown that you are decreasing the pressure from 1 atmosphere minus 0.5. So you are let us say rotating your radio knob in the negative direction and you are rotating by an amount such that 1 atmosphere becomes 0.5 atmosphere. This minus 0.5 is what is called as pressure correction but when you have done this you will get new mass fluxes. So when you do this you get a new mass fluxes and this is the second iteration. So when you change the 1 atmosphere to 0.5 atmosphere you got this new mass fluxes whether they are obeying mass conservation no. So this is not the final value this is the second iterative value. Now what is happening more mass is coming in now what you have to do you have to increase to resist the inflows. So now you increase by 1.5 so that the pressure becomes 2 atmosphere. Right now I am showing this philosophically. So this way the picture repeats and all this picture corresponds to different iteration and there is not one knob if there are 25 grid points you have 25 knobs. So this is the coupled system non-linear system. But fortunately we are fortunate that finally we get a picture where mass balance is obeyed in an approximate sense is approximately 0 mass conservation is approximately 0. This I am showing you philosophically in the coming lecture I will show you specific equations. We have equations to calculate that from 1 atmosphere how we obtain minus 0.5 and we get 0.5 we have expression from which we correct the mass fluxes from 1 kg on the west phase how we obtain 3 kg on the north phase 2 kg how it changes to 1 kg we have specific mathematical expressions. So here I am showing you philosophically through an animation in the coming lectures tomorrow I will show you mathematical expressions. But in books if you directly see those expressions you do not get feel so through this pictures and animation I feel you get a better understanding of what that equations are doing. So I will discuss specific equations which are written symbolically here I will discuss more tomorrow morning. Thank you for your attention.