 Now, after formulation, let us go that how you will develop into code. This is a very important thing, maybe you will not find this in any book, but this helps a lot in code development. So, again here I will use lot of symbols to explain, to give you a picture, give you an animation and the same idea is exactly being implemented in the code which will be there in the afternoon session. So, we will take a simple problem, a Cartesian coordinate problem and this problem is called as the benchmark problem. This is probably the most widely used validation problem in Cartesian coordinate system. This is called as lead driven cavity flow. So, the idea is there is let us say square box, the left wall, bottom wall, right wall is stationary, the top wall is like a conveyor belt which is moving with a constant velocity. So, let us suppose the top wall is fitting tightly over this box and this is moving with a constant velocity u naught. This is called as the lead driven cavity and this is called as top wall is called as the lead, because that is driving the flow. If this is stationary, then the fluid will be stationary. Now, this is the way we draw the stagger grid. These are the grid points for pressure and temperature. This is for the u velocity. This is for the v velocity. This is for the boundary points for all the three variables. These are the running indices. Now, let us start that what is the first step? How to predict u velocity? First, we will do prediction using the pressure from the previous time level, mass flux prediction, then mass imbalance from the mass imbalance, pressure correction and so on. Let us first predict. When you want to predict u velocity, which control volume you have to use? u control volume. This will be the, now when you want to have u control volume, it will not be near to the boundary. It will not be equispaced vertical lines. You will have half size control volume on the left side and the right side for if you want to have u control volume. This is the u control volume. Compare this with the pressure control volume. This is staggered. This is at the boundary. Let me ask you a question. If you want to predict, how many points you will predict? How many green squares you have right now? Let us say, total you can say in the x direction, 6. y direction? y direction you have 7. Let me ask you that the same control volume, if you see in pressure, how many grid points were there? Do you remember? In the pressure, in the case of pressure, it was having, it is the same division which I am showing you earlier, I had shown earlier in conduction also. So, there how many points we had maximum? i max, j max, I was talking repeatedly 7, 7 by 7 and inside they were 5 by 5, correct? I am taking the same division. So, here pressure or mass conservation control volume will be 5 by 5, 25 control volume for mass conservation. Now, by looking into this, how many points you will predict for u velocity? Can you tell me? There is no need to predict at the boundary because we have a u is equals to 0, u is equals to 1 on top wall and u is equals to 0 on all the wall. So, there is no need to predict at the boundary. You need to predict only on that green squares which are sitting inside. So, how many are there? 4 in the x direction, 5 in the y direction. I am telling this because in CFD it is very important that your looping should be proper. So, to give you a feel of the numbers, I am showing you like this. Just note that continuity 5 by 5 x momentum 4 by 5. Now, what I am drawing this? What are this? These are, these are horizontal phase center. So, when you are applying the, at the phases, when you apply law of conservation of x momentum, you will need mass flux at this phase center to calculate advection. You will need to use advection scheme at this phase center. Basically, you calculate mass flux, advection flux and diffusion flux at this phase center. What is the diffusion flux in x momentum equation? Discuss stresses in x direction. What is the advection flux in this case? x momentum flux. So, in this yellow square, these are the vertical phase center of u control volume. Note that. And these are vertical, these are horizontal phase center of u control volume. At this phase center, you will calculate the viscous stresses acting in y direction. Sorry, viscous, actually in the yellow squares, you will get normal viscous stress acting in the x direction. In this pink square, you will get shear viscous stress acting in the x direction. Note that. So, this way I am showing you the horizontal and vertical phase center for the u control volume. And this is your u control volume, where v is sitting at the corner. This is the convention we had been following consistently. For the x direction on the positive phase, we take j, i. At this phase, it will be i minus 1. And the fluxes which you will encounter is mass flux, advection flux and diffusion flux. And this is for the y direction. We have this fluxes in two directions, x direction, y direction. For x direction, we use vertical phase center. For y direction, we use horizontal phase center. Note that I am showing you the loopings also. See here, it is 2 to 7 minus 2, interior control volume, green squares. Let me go back and show this to you. Green square, what is the loop? Just note that. 4. And what is this? 2 to 7 minus 2. 7 minus 2 is 5. 2 to is 5 is how much? 4. So, it will cross, scroll through this 4. And what is in the j direction? 2 to 7 minus 1. 2 to 6. That is 5. These are 5. Now, when you look into this vertical phase center, j is same. j is in the j direction, you have 5 yellow squares. But in the x direction, you can see it is not 4. It is 5. Green square is 4 in the x direction. Yellow square is 5. So, instead of 2 to i minus 2, this is 1 to i max minus 2. If you look into this pink square in the j direction, now, in the j direction, that green was 5. Whereas, pink in the, this is 6. 1, 2, 3, 4, 5, 6. This varies from 1 to 7 minus 1. Note that in CFD, first you need to understand infinite paper. You have to do this thing exactly. And in your code also, you need to have exact. There is no approximation. If you do one mistake, one i, you make i plus 1. You will not get the solution. The solutions will diverge. So, here nothing less than perfection is expected. But if you do this homework properly, before doing sitting in front of computer, you will do less mistakes. Before coding, it is suggested that you do, try to create proper looping. So, this type of effort is just shown to you, so that you do not struggle and get frustrated in coding. Let me tell you. It is very easy to get frustrated by coding that is very difficult. But the thing is that it requires a lot of focus in understanding, writing it out, because everything is expected perfectly. So, this looping basically we used here. So, this we used for that yellow square location. This is for that pink square location loop, because this where the vertical phase center, horizontal phase center. And that way we calculate mass fluxes. And note that this mass fluxes I had mentioned earlier that for u control volume, you have to do linear interpolation. So, for mass flux linear interpolation I had already mentioned. This is the advection flux using first order point scheme. And this is the diffusion fluxes, which is viscous stresses in the x direction. So, once you have calculated the fluxes, those were the loops to calculate the fluxes. This we did earlier when we had done conduction heat fluxes calculation, advection heat fluxes calculation. Yeah. Is there any question reason is there to use j M i in the programming? Many times when you are plotting, it is just for plotting convenience. There are some software which take this type of data structure. So, it is just for that nothing more than that. So, once you have calculated the fluxes for the u control volume, then you do the balance. This is the total balance, total advection. Physically what is in that rate at which exponent term is lost by the fluid in the control volume? This is the total viscous force acting in the x direction. This is the total pressure force acting in the x direction. So, using this total summed values finally, we predict velocity u velocity. Now, we after predicting u, we will predict v. And these are the control volumes for where you are getting half control volume. There we were getting in u left and right wall near to the left. And here you get bottom and top. And these are your v velocity control volume. 5 in the x direction, 4 in the y direction. So, 5 by 5 for pressure, 4 by 5 for u velocity, 5 by 4 for v velocity. These are the boundary point. So, these are your interior control volume. Note the loop for interior control volume. 2 to, in the i direction, it is 2 to 7 minus 1, 2 to 6, that is 5. In the j direction, 2 to 7 minus 2, that is 5. 7 minus 2 is 5. 2 to 5 is 4. You have 1, 2, 3, 4 in the y direction. So, what I mean is that you may be developing a complex code. But what I want to highlight here is that by doing some small simple diagrams, you can generate proper looping. And this looping is very important as far as coding is concerned. So, in pen and paper, you can create some, right now I have taken this many control volume. This you can even do with a 3 by 3 control volume also. That way, I started today morning. But this exercise is very help to decide what will be the proper loop. So, to create a final and navier stroke solver, there is lot of exercise or homework to be done. Lot of surgery needs to be done. So, this is some of those. I am trying to see most of those. In fact, what are these? These are vertical face centers, where we will calculate the fluxes in the x direction. These are at the boundary. Vertical face center proper looping has been in the j direction, it is same, but in the i direction, it increases by 1. You have 6 yellow inverted triangles and 5 red inverted triangles. Note that. This is for the horizontal face center. And this is a pseudo code. Note that I am marking this g m x minus 2, because this is different from what we were doing in conduction and advection. There it used to be g m x minus 1. This is to calculate in the x direction. This is to calculate in the y direction. Fluxes to calculate the v velocity. So, for v control volume, we have calculated the fluxes in the x and y direction with proper looping. Once you know the fluxes, you do the balance. So, in this slide, the balance is shown. Once you do the balance, then you can predict the velocity for the new time. You can predict the velocity. So, velocity prediction is done. Now, what is the next thing? We have to solve that. So, we have predicted velocity. Then, what you do? You do calculate mass imbalance. Where you calculate mass imbalance? Continuity equation. So, let us now go to continuity equation. So, this is your continuity equation. 5 by 5 control volume. 5 by 5 for pressure, 4 by 5 for u velocity, 5 by 4 for v velocity. You can call this as grid point for u velocity also, but you can also call it as a grid point for vertical face center, where you will calculate the fluxes. So, now, this is the code for character. So, we are done with the prediction. Now, this is the correction. So, based on the predicted velocity, we calculate the mass imbalance and check whether the mass imbalance is less than practical 0 value. Let us say 10 to power minus 3. If it is not, if it is satisfied, then go for computation of temperature. If it is a force connection problem, else continue. So, if mass imbalance is not less than 10 to power minus 3 or practical is 0 or convergence criteria, then you calculate pressure correction as a function of mass imbalance. Note that here you are only doing one iteration and you get an updated value of this value of pressure. With an updated value of pressure, once you get a pressure correction for one iteration, that iterative, that updated value of pressure correction you use to calculate updated value of velocity correction. Note that in this loop, you are always updating the predicted value. That is what he was pointing out that it is a simultaneous solution. It is not that we are solving the pressure correction equation till convergence and then we are adding the velocity correction to the predicted value. This we are doing in the same loop. So, we do one iteration of pressure correction, then update the predicted value. Using the updated value, we again calculate the new mass imbalance. This is done inside this loop. To make it more clear, what I mean is that this pressure correction poison equation, the right hand side changes after every iteration because this velocity prediction after each iteration of pressure correction we are updating. So, I talked about the velocity prediction, x momentum, equation, y momentum. So, u control volume, v control volume, then the pressure control volume where we do mass balance. Once you do all this, what you get? You get the correct velocity for a particular time step. Once you have got the correct velocity for a particular time step, you can, if it is a heat transfer problem also, then you solve the energy equation. Now, I am discussing the implementation details for energy equation. So, the pressure control volume is used for the calculation of temperature also. Note that for the symbol which I was using here, for u velocity I was using this as A x 1 and for y momentum this was A x 2. Here I am writing as A x 3. This M x, M y I am using for continuity and the energy equation because they are same. I am using M x 1 for x momentum and M x 2 for y momentum. This is the symbol which I am using. Actually, this loop is very much similar to what you have seen in conduction problem. So, for solving energy equation, first you calculate mass flux. Note that to calculate the mass flux while solving the energy equation, I am not doing any interpolation for the normal velocity because in that temperature control volume faces, the velocity grid point is already sitting which I was doing. That is why it was M x 1 in for x momentum and M x 2. There I had to do interpolation, but not for this. Then, this is the first order apparent scheme and this is the diffusion flux which is heat gained by conduction. Then finally, the balance is done to calculate the total enthalpy lost by the fluid control volume, total heat gained by conduction, total heat gained by volumetric heat generation and finally, obtained the temperature for the new time state. So, with this I complete the implementation detail where I showed you through some picture, animations, various how basically main thing message was how to do proper looping and I had also shown you pseudo code. Sudo code means I am not using the syntax of a particular programming language, but you can take this and use the proper syntax of a programming language and convert into a code. Anyway, I am not discussing the pre-processing and post-processing steps because I believe that those are straight forward, straight forward. Any question on this? Let us solve this problem. In boundary conditions, we specify only no slip condition and at the top velocity. At the boundary, the pressure need not be specified. You are solving pressure correction as a Poisson equation. You will not need boundary condition for pressure in staggered grid. I will let you point out, but you will need boundary condition for pressure correction. So, the next topic is boundary condition. Now, velocity boundary condition is very straight forward. So, you know it there is no slip condition. So, u is equals to 1 on the top wall and 0 in all the other wall, v is 0 at every wall. So, let us discuss about the pressure and pressure correction boundary. Actually, why you will need pressure boundary condition? If you are where pressure is coming in your formulation in the x momentum and y momentum equation. So, you have to look into the border u control volume and see whether the pressure which are involved are at the boundary or away from the boundary. So, let us look into that. Let us look into this. Do we need pressure boundary condition? The question is this. Do we need a pressure boundary condition? Let us try to. So, out of this, let us how many we have? 4 in the x direction, 5 in the y direction. The roll of boundary condition comes when you look into the border control volume. So, this are interior. Let us look into the border. Border means this 4, this 5, this 4 and this 5. Let us try to understand what pressure they will demand. This column of grid points will demand pressure from here and here and this pressure is not at the boundary. When you go on this side also, this will demand pressure from here and here. This note that this is u control volume. You need pressure only on the vertical phase center not on the horizontal phase center. So, do not say that I will need pressure here for this control volume. I will not need pressure here because this is a u control volume. In the x momentum equation, you need pressure force in the x direction where you need pressure only at the vertical phase center not on the horizontal phase center. For that I have to go to v control volume. There the picture will change entirely. So, what you see is that when you are taking this interior control volume for u velocity, you do not need pressure of boundary. Any problem you take, it will be that this is convenient only. Any problem you take in Cartesian coordinate. We will not put the accuracy. No. In fact, that is a good thing. If we avoid pressure boundary condition, many times pressure boundary condition causes lot of convergence issues. In fact, if you look into the more complex problem like multi phase flow problem or moving boundary problem, people find this as better as compared to the collocated view. So, pressure boundary condition will not come into picture because your calculation of pressure force in the x momentum equation and y momentum equation, the expression is such that it does not demand pressure at the boundary. Boundary condition comes because when you discretize for the border control volume, it demands a boundary value. Similarly, if you go to the y momentum equation, I can show it to you. So, in y momentum, when you want to calculate pressure force, pressure acts in which phase centers? Horizontal phase center. Pressure force in the y direction acts in the horizontal phase center. Now, here you see at horizontal phase center, it will be this to control pressure force in this control volume. It will demand pressure from here and here and this is away from the boundary. Here, here, here and here. You will not need pressure in the horizontal phase center, note that. So, that way in your discretization, discretize the equation, pressure does not come at the boundary. So, you do not need pressure at the boundary. Just wait. That is in the equation itself. It is del p by del x in x direction momentum equation and del p by del y only. In the y direction equation, it does not require del p by del x. Correct, correct, correct. Yeah, I mean, do you have a question? Sir, in this case, the problem does not specify pressure. There is one velocity and three stationary walls. Correct. But if problem involves pressure, then how things will change? If pressure is specified on one of the boundaries, say inlet or outlet. Then we use appropriate boundary condition, the pressure correction boundary condition, which takes care of it. What is basically saying is that there are problems like if the fluid is coming out from a pipe under ambient condition and the pressure is atmospheric. How does this take care of it? There also we do not apply the pressure boundary condition, but we apply pressure correction boundary condition, which takes care of it, which I will discuss next. That is a good question. Now, with this I had shown that we do not need pressure boundary condition. Now, let me and I had also mentioned that we need pressure correction boundary condition, because we are solving a Poisson equation. So, let us go to the pressure correction boundary condition. The first what I am showing you is what you said just now. If the pressure is given at some boundary, if something is given, there is no need to correct. So, the correction should be 0. So, if we know the pressure at some boundary, we apply a pressure correction is equal to 0, that automatically takes care of that prescribed pressure boundary condition. Is that clear? Boundary condition for u velocity. Boundary condition for u velocity will be required, because when you want to calculate, let us say mass flux here, you will need u velocity here. So, you will need mean of this and this value. Your first unknown we will be there on the second page somewhere and in the, no we are quite staggered in the text direction. Yeah, no you need here mass flux. So, at the phase center, you need u velocity on the horizontal vertical phase center and v velocity on the horizontal phase. So, for that you have to do interpolation. So, in that phase that interpolation. Yeah, u velocity boundary condition will definitely be used. Strictly speaking, this boundary condition needs to be used, because boundary condition acts as a jury to the problem. They dictate the solution. So, the pressure is prescribed, pressure correction boundary condition is 0. Now, most of the time, what you see is that the normal velocity is always prescribed. And actually, we have an equation from which we can, if we know the normal velocity, we can derive or we can degenerate the boundary condition for pressure correction. Is that clear? What I am saying is something like this. Let me show you that actually, the way we have written here our velocity correction. If you want to express it in discrete differential form, it is something like this. I would like you note this equation. I will write this in paper and I will show it to you. Just see this equation, after writing I will show it to you. So, this equation I am showing you here. This is up prime is equal to delta t divided by rho delta vp. This can also be written as up prime n plus 1 rho delta vp divided by delta t is equals to pp prime minus pe prime into delta y. What will be the differential form of this equation? Can anybody tell me? What will be the differential form of this discrete equation? What will be the differential form of this discrete equation? This is a, let us say, finite volume discretized of a, let us say, differential term. So, what will be the differential part of this? Yes, rho del u prime by del t is equals to del p by del x. So, suppose u prime is prescribed, u prime is given and it is not also changing with respect to time and is not function of time. Then what will happen? This will become 0. Is that clear? So, then what we say is that this equation is there in a semi-explicit method or in semi-implicit method also. So, if u prime is prescribed del p prime by del x is equals to 0 and if v prime is prescribed del p prime by del y is equals to 0. In general, I can say if the normal velocity is prescribed in a boundary del p prime by del n is equals to 0. This is specific. This is u is the velocity in the x direction, gradient in x direction, velocity in the y direction, gradient in the y direction. This is the velocity in the normal direction, this is the gradient in the normal direction. So, we obtain pressure correction boundary condition indirectly. Indirectly means we first see what is the normal velocity at the boundary because in almost all the problem either normal velocity is known or pressure is known. Actually, normal velocity you know even in case of inviscid flow. What is the boundary condition in case of inviscid flow on a solid surface? For an inviscid flow on a solid surface, what is the boundary condition we have? Normal velocity. Normal velocity. What is no penetration? So, by knowing the normal velocity if at the boundary it is prescribed we take and in all the problems this thing happens. Either pressure is known or normal velocity is known. So, pressure is known we take p prime is equals to 0. If normal velocity is known we take the normal gradient of pressure correction is equals to 0. I would like to point out that although I am not talking of pressure boundary condition here, but later on we will use it in let us say collocated grid which is the next topic after t session. Pressure boundary condition is an approximation. This is not an approximation. We use boundary layer theory and come up with pressure boundary condition and which is an approximation. Is it that in this case the pressure boundary condition, pressure correction boundary conditions are always Nyman boundary condition? Not always. If you have an outlaid boundary where pressure is prescribed it is not Nyman. This is this is not Nyman. If u n is prescribed then it is Nyman. Yeah, then it is Nyman. Otherwise it is not. Otherwise it is initially. I would like to point out this that is staggered grid for methodology. This is more robust. After p session I will talk of another method which I believe as per as formulation is concerned it is not that robust. But there is some other advantage due to which people jumped into it because it is easy to use that approach in complex geometry. This staggered grid if you have to do in complex geometry or unstructured grid your life is horrible. Just imagine if you have unstructured grid and this type of staggering how you can manage. For structured grid this is good. Even in case of complex geometry where your control volume shapes are changing neighbor by neighbor this formulation. Staggering creates I will say staggering is easy as far as the way I showed you but staggering in case of unstructured grid or complex geometry becomes quite complicated. So to solve complex geometry problem people came up with collocated grid method. But nowadays we have personally realized that if you want to use for we are developing codes for moving boundary problems and the multi phase flow simulation. This collocated grid as I believe that it is not mathematically that rigorous it creates lot of convergence issues. So for the lid driven cavity problem the boundary condition is this. You do not have prescribed pressure at any boundary. You have all four solid wall. So all at all four solid wall you have a Neumann boundary condition. It is like an insulated boundary condition del p by del n is equals to 0. And this is the discreet form of that pressure correction boundary condition. And this is at the velocity boundary condition. Any question on this? Is it clear? This is basically the discreet form of normal gradient. On the bottom wall this is equal to that border value. This is at the top wall sorry left wall i is equals to 1 is equals to i is equals to 2. I max equals to I max minus 1. j is equals to 1 is equals to j is equals to 2 and so on. So finally solution algorithm. So these are basically the user input, thermo physical property, geometrical parameter, number of grid points, boundary condition inputs, grid generation, stability criteria using stability criteria for time step. I would agree like to point out that you do not have any stability criteria for the Neumann stoke equation. We have stability criteria for pure diffusion and pure advection. And we get two different time step and we take the minimum of the two but it is not guaranteed to convert. That has been implemented in the code also. Then you set the initial condition because we are fine using a transient approach. Set the boundary condition. So this is the pre-processing stage. Then you predict you predict u velocity. You calculate fluxes and the source term which is the pressure force in the x direction. Calculate total advection, total diffusion and then finally predict the velocities both u and v velocities. First you do this for u then repeat it for v. So in this slide what I am showing you is the prediction of u velocity and v velocity. So the previous slide was on pre-processing. This slide is on velocity prediction. And the next slide will be on solving the continuity equation and updating the predicted equations iteratively. So the previous slide was momentum equation. Here we check for the mass imbalance. If there is some mass imbalance we get the pressure correction. Then calculate the velocity correction as a function of pressure correction and get the updated value of the mass flux predicted prediction and then this continues. I had shown you implementation details separately for each of these predictions, corrections and this is the solution algorithm step by step procedure. So in the afternoon session you will have a lab where you will have a code which you can see where you can see this has been implemented exactly. So with this I had come to the end of this session. If you have any question I will be happy to answer. In the initial condition we usually prescribe u and v velocity. So do we need to prescribe the pressure correction also in the initial condition? Yes we are following an iterative method. I would not say it is an initial condition. You can call it as an initial guess. Initial guess of the pressure correction. Because when you are solving this pressure correction Poisson equation you need an initial guess not an initial condition. Any other question? For pressure and velocity. So this question is whether we need to use any under relaxation factor for velocity, pressure and velocity. The answer is that if you are following a semi explicit method so your prediction is not iterative. So you do not need any under relaxation for that. But in case you are using a semi implicit method you need it whereas pressure correction is concerned we are solving this equation and in this case depending upon a problem you may have to use it. Sir in this case you divide the domain first into pressure control volumes. Yes. And then staggered the u and v controllers. Correct. Whether this will be the case for all problems or yes. Although I have shown this staggered grid in a Cartesian coordinate system but this can be used in cylindrical coordinate also spherical coordinate also. But if you want to use for complex geometry people do not use for complex geometry because there it becomes too complicated. For standard coordinate system this staggered grid staggering you can do and neighboring it is possible even in non-uniform. Right now I am showing you for a uniform grid. I am taking equispaced horizontal and vertical lines that way it is uniform. But you can do it for non-uniform but in standard coordinate system. Even for Cartesian always first pressure control volume and then staggered. Even if velocity is given suppose on the left side. But anyway we are using that velocity boundary condition we are using anyway because I had mentioned earlier that you calculate mass flux you need to interpolate the boundary value. So the boundary velocity boundary condition will be used definitely. Right now I am not discussing physical problems the results from the let us say lead driven cavity flow. In the next topic I will be discussing three types of in fact four types of problems. One is isothermal flow where is there is no heat transfer. Second is force connection. Third is mixed connection and fourth is natural connection. But I will take the same closed cavity problem and in the afternoon you have a lab session and you have the problems which are given are all these four types of. I would like to point out that in the afternoon session we have given a code where you can solve let us say isothermal force connection, mixed connection, natural connection. We have developed that code in a non-dimensional form. We have a non-number, Grashev number, Prandtl number, Rayleigh number are the governing parameter. So the user input is in non-dimensional form and we have included also the benchmark results for all this problem. For all this four types of problem we have taken from the research paper the published data and when you run the code we have we are showing you that overlap results also. However as you have a limited time this code if you want to get a grade independent solution it may take 2 hours, 3 hours or 4 hours. We are giving you telling you to run the code on a coarser grade. You will see lot of inaccuracies but do not get bothered about it. If you run it we are also saying take it as a homework and run on fine grade. We had told you in that problem sheet. So I just wanted to inform that. So we will beg for tea and come back and start the new topic. Thank you for that.