 So, this lecture is on computational heat conduction. So, we started with introduction then moved on to grid generation and then we will start with heat conduction problem. So, pure diffusion problem and I will start with finite volume method for heat conduction problem. So, there I will show you how to obtain linear algebraic equation from control volume using certain approximation. Then I will talk about the two different methodologies solution methodology explicit method and implicit method, thereafter I will discuss implementation details how to convert it into a code. I will be little more detail as far as code development is considered in this topic as compared to the grid generation topic. Finally, I will discuss the solution algorithm and then I will show you some example problems. So, this slide shows that the finite volume method which is being discussed uses an approach which is different from what you see in what you see most of the CFD books. This is what you teach in your undergraduate fluid mechanics for heat transfer, how to obtain governing partial differential equation which you call as the Navistok equation which includes continuity as well as momentum equation. So, from control volume in undergraduate courses typically student learns how to get partial differential equation. Now, when you go for CFD you are interested into algebraic equation. Now, what you see that in most of the CFD books in fact all the CFD books they start from this governing partial differential equation and they show you how to obtain the algebraic equation, but the same set of algebraic equation can be also obtained from control volume. This is what will be shown here specially I had found that for the undergraduate students as we had discussed yesterday that many colleges you say that they have you have CFD course for undergraduate, but the students are not taking it up and one of the major reasons you feel that they feel it is quite mathematical they find it specially they are it is offered as an elective in the last year and the students go for easier option. So, in this we will follow an approach which is not that mathematical although some mathematics will definitely be involved. So, we start with control volume. So, rather than starting. So, the difference is rather than starting from differential equation we start with some picture some control volume and reach to the differential equation. So, this is we discussed in this course. So, this is the plate of length l 1 in the x direction l 2 in the y direction and this is the equispaced vertical lines equispaced or simplest form of the grid generation for heat conduction in plate and this is the centroid of the plate which is the cell center for finite volume method. So, there are 25 yellow circles which are the interior points for the plate. So, let us suppose you are interested you have taken this plate from a furnace at let us say 200 degree centigrade and it is subjected to certain thermal boundary condition. So, let us suppose left wall is subjected to 100, bottom wall is to 200, right wall to 300 and top wall to 400 degree centigrade or there is a boundary condition. So, you want to determine the temperature distribution. If you are doing an unsteady simulation you can create even a movie of the heat transfer. This is a typical computational stencil which is involved in the electric equation. This is the boundary control volume, border control volume this I had discussed in my previous lecture an interior control volume. In the last lecture we had shown you how to obtain the x and y coordinates of the vertices of the control volume. Once you know that you can calculate certain geometrical parameters which will be coming in finite volume method. What are those geometrical parameters? Width of the control volume in the x direction in the y direction the phase center. We have four phases earlier I was saying left phase, right phase, bottom phase, top phase. In CFD we call them as west phase, east phase, north phase and south phase. Then you have certain neighbors. So, a small letter this small w represents phase center capital W represents west neighboring cell center. This is south cell center, east cell center and north cell center. Other than the width of the control volume you also need the distance of a cell center from its neighbor which is used in calculating the temperature gradient at the phase center. So, once you know the coordinate of the vertices from the grid generation you can obtain the coordinates of the cell center. Then you can obtain the width of the control volume the distance between the cell center and its neighbors. Width of the control volume you can calculate through this expression surface area you can calculate through this expression volume you can calculate this expression. These are two-dimensional control volume. So, this geometrical parameters are later on I will show you that they will come as a coefficient in the algebraic equation. In grid generation after we obtain the coordinates of the vertices we obtain this parameter because this parameter will come into the algebraic equation as coefficient other than the thermo physical property. So, in my first lecture I showed you how from a control volume how to obtain differential equation. Now, in this second lecture what we will do is that we will not start from differential equation. We start from the same control volume and learn how to obtain an algebraic equation. What is the law of conservation of energy? Rate of change of internal energy or enthalpy of solid inside the control volume is equals to heat gained by conduction which is total heat inflow minus outflow plus heat gained by volumetric heat generation. So, total heat gained is equal to rate of change of internal energy. So, we start with certain words now we have to convert this words into algebraic equation earlier we had converted into differential equation. Now, let us suppose take this room and if you want to represent the rate of change of internal energy of the fluid inside this room. Rate of change of internal energy of the fluid inside this room let us suppose temperature inside this room is varying with respect to time. So, if you take one point temperature is varying with respect to time. How do you represent rate of change of internal energy? How do you represent internal energy? Mass into specific heat into temperature this is internal energy. If you want rate of change time rate of change and you do the temporal derivative. So, del by del t. So, rate of change of internal energy is equals to del by del t of rho into volume this is mass specific heat and temperature. If the fluid is incompressible density will not vary with respect to time specific heat does not vary with respect to time. If you take a fixed volume this will let me put it like this strictly speaking this can be expressed like this rho C p volume integral of del capital T by del small t. Do you agree with this? This rate of change of internal energy rather than writing like this I can write it as a volume integral. Because this rate of change of temperature varies from point to point in a volume is correct. In this room you have infinite number of points and at each point this term is varying. So, you can write it as a volume integral. So, in this you have to use mathematics to calculate this term. What is the variation of temperature gradient in space in a volume? Let us say in this room you have you put some temperature measuring instruments and let us suppose wherever you are sitting just at that location measure the temperature with respect to time. So, consider each one of you as a control volume. So, it will vary from point to point, but if I say that out of let us say 50 people I want to take only one value which value I should take who is sitting where in this room that corner that corner of front centroid. Anybody anyone would like to say with a different answer to this? Corner anybody would like to prefer or front or back or centroid is good enough. So, who is sitting at centroid? The cameraman is sitting at no. Why centroid? Let me ask you that question that approximation will be good that will depend on the size of the room or not. It would be good approximation on large size room or small size room that is the mathematics which is involved. If you understand that this understanding is very much similar to I will give you an example. Let us suppose you see a curve when I said volume the variation is with respect to volume, but for simplicity I am just showing you temperature or some function is a function of x. Let us suppose this function is varying in a very non-linear fashion. This centroid which you are saying is basically a linear approximation what you are using assuming is in this room in this volume the variation is linear in all the direction x, y, z direction. So, the idea which we are using is very simple. Instead of that volume if you take certain distance element in this region what is although let us suppose this function is varying from here to here, but if I want to represent by one point value what is that best one point value? Suppose this is x 1 this is x 2 I want one value in between x 1 to x 2. So, which location temperature you would like to prefer? At the middle. So, this middle or that centroid which you said is the same idea mathematically. So, you can if you can understand this curve. So, here what is varying is temperature gradient in a volume here this is a just a function. In CFD note that we use some approximation, but that approximation is good if your size of the control volume approaches 0 it is very small. That is why you want to have large number of grid points. You want to have very small volumes because smaller is the volume better is the approximation. And this idea is very much similar to what you study in your schools that a function may be very non-linear, but if you take very small part of it its variation is almost linear. This is the level of mathematics which will be used in the approximation. So, the first level of approximation is the average value of this temperature gradient in a volume is equal to the value at the centroid. So, the first level of approximation is in a volume I showed you earlier I had shown you here. In a volume the change of this what is this term rate of change of temperature in this volume. So, we need to do volume averaging of the rate of change of velocity in case of, but here it will be only temperature. Volume averaging of for the rate of change of temperature in a control volume as a value at the centroid. This is what we discussed just now. This is what is order of this approximation. This is mathematically it is called as second order approximation. Now there is a second level of approximation. The first level of approximation is volume averaging of rate of change of temperature. The second level of approximation is how I can obtain this del capital T by del small t in algebraic form. If you know finite difference method and time is much in the forward direction. Always whenever you solve a unsteady problem you have an initial condition. From initial condition so typically let us suppose you want to create a movie of flow or heat transfer always you have a situation you know the initial picture from initial condition and boundary condition. From initial condition you obtain the picture for let us say if your time resolution is 1 second using the picture of 0 second you obtain the picture for 1 second. Using the picture of 1 second when I say picture it is not only pictorial information which is used in the codes which is the fluid dynamic information like velocity distribution, pressure distribution. So, time always marches in the forward direction. So, here one first order approximation is that we do forward differencing in time. This is derivative with respect to time we use superscript for as running indices as i, j are running indices for x and y direction n is the running indices in the time direction. So, we do forward differencing and represent this temperature gradient as T at 10 plus 1 minus T at 10 divided by delta t at the centroid. So, this centroid representation first we this is the first level approximation volume averaging in the second level we do first order forward differencing. These are the two approximations which are used to represent the unsteady term as this. Any question on this? Mathematically how it can be said that the first level of approximation is the second order approximation. Later on when I will be discussing the advection schemes there I will show you mathematically that it is second order accurate. When I will be discussing the central difference scheme because this is like a central difference scheme where you are basically type doing some type of averaging. There mathematically I will show what we will do is that we will take a quadratic curve and then show that it is a second order accurate. Any other question? So, this is the you have a question finite volume method discretization of one term which term is this unsteady term. Right now I am showing it for rate of change of internal energy, but we will use the same thing for rate of change of x momentum rate of change of y momentum. Same approximations are used there also. Now let us go to the diffusion term. What is diffusion in conduction? Conduction heat transfer. So, what you do in heat transfer course you have q x, q x plus delta x, q y, q y plus delta y. But here now when I am doing when I am in CFD I will not write this as q x or q x plus delta x because I want to write this q on this surface as one point value. What is the best point? Earlier we talked about volume, now it is surface area. Answer remains and is there answer different or same. If I earlier my question was that I want to represent in a volume some function. You said that the best point is centroid. Here again in CFD as CFD we use a discrete method. So, we want to represent the variation as a one point value. What is the best point? Centroid of the surface area which is shown here. This small e is the centroid of the phase, east phase. So, the approximation which we are using is surface averaging of heat flux at the control volume phases. How do we represent in general the rate of heat transfer? Let us suppose perpendicular to this wall the length is unity. So, earlier I had shown you this is volume integral. Here it is actually surface integral, but in the z direction the length is unity. So, that is why I am because this q is varying with respect to y in general. As in earlier case your temperature gradient was varying in volume point to point. So, here point to point this heat flux is varying is the function of y and in the z direction length is unity. So, total rate of heat transfer can be written as integral of q d y where what will be the limit of this sorry from this y y s e to y n e. This can be expressed as average value multiplied by delta y and what is this average value we are taking centroid value. Similar thing we do at other phase centers. So, this is expressed as q n where any centroid of the north phase, east is centroid of the east phase and so on. This derivation is analogous to del by del x of q x plus del by del y of q y. Do you agree with this? When you teach conduction you also start with this q x, q x plus delta x, q y, q y plus delta y. You do a balance. How do you write this q? You write as q x into delta y. This is q x plus delta x into delta y and when you do the balance you write del by del x divided by volume anyway del by del x of q x plus del by del y of q y negative of that. That was differential form of the diffusion term and this is algebraic form. There is no differential term here. What is this term? This is q e into delta y. Here we are doing a balance where is in minus out. So, it is q w into delta y plus q s into delta y minus of q n into delta x minus of q e into delta y. So, this is an algebraic equation. It is not a differential equation. This completes the first level of derivation, first step of derivation analogous to del by del x of q x plus del by del y of q y. What is the second level of derivation in undergraduate courses of it transfer? How to obtain del square t by del x square plus del square t by del y square from del by del x of q x plus del by del y of q y. So, you use in the second step basically you substitute Fourier law of it conduction into that first derivative and then that first derivative gets converted into second derivative. Here the second level of approximation is we will do the same thing. Here also we will use a Fourier law of it conduction, but the difference is Fourier law of it conduction is in a differential term. We have to express this in algebraic term. This is the second level of approximation. So, the second level of approximation is discrete representation of Fourier law of it conduction. Let us suppose this is q e. This q e is actually in which direction x direction normal direction. So, it will involve d t by d x. So, to obtain d t by d x what you do? You draw a horizontal line. If you draw a horizontal line what you are finding? That horizontal line is intersecting this point and this point. So, your face center is lying it in between two yellow circle and you have to calculate gradient in that direction only. So, this becomes d t by d x equals to t capital E minus t small p divided by this delta s e is a distance between this two yellow circle. By finite difference method which is this what is this expression? You want to calculate gradient at a point using one point behind one point ahead. This is called as what central difference is. What is the order of accuracy? Second order Taylor series expansion I had shown in the previous slide. So, this is a second order finite difference discretization. Actually, when you want to calculate gradient here you do the same thing at all the phases. So, here if you want to calculate you take temperature here minus temperature here divided by the distance between the two points. Similarly, here but this type of variation basically what you are assuming is that between two yellow circle the variation is linear. That is why you are using this expression. So, but you are not making sure that there is a continuity between the linear relation. So, that is why it is called as piecewise linear approximation. If you follow Patankar books he uses this approximation because there is we are assuming linear variation between two yellow circle. But we are not sure that the linear variation between this two yellow circle and the next two yellow circle the slope is same. That is why we say that the variation is piecewise linear. There is not continuity in the linear slope of the variation. There is a discontinuity. So, this is the piecewise linear second order approximation. So, we have two levels of approximation. By using first level of approximation we get equation similar to del by del x of q x plus del by del y of q y which I shown in the previous slide which is this equation. What is the approximation we do surface averaging of the heat fluxes. And this is the second level of approximation where which is a discrete representation of normal gradient of temperature at the phase center. Note that this is a discrete representation of normal gradient of temperature. Then what you do in heat transfer you substitute this Fourier law of heat conduction into del by del x of q x plus del by del y of q y. We will do the same thing here. We will substitute this into this expression. Once you substitute you get this equation. What are the different terms? This is the rate of change of internal energy. What is this n plus 1 is the new time level. This is the old time level. And what is this term? This is the rate of this is total heat gained by conduction. And what is this term? This had not used any approximation. This is volumetric. It is per unit volume. You basically multiply this term with volume. Now, there are two types of method one which is called as an explicit and second is an implicit method which basically depends upon whether you are calculating this total heat gained by conduction using the temperature of the previous time level or the new time level. I will give you an example. Let us suppose you took a plate from a furnace at 250 degree centigrade. Let us say you have subjected that plate to a boundary condition. Let us say left wall 100, bottom wall 200, right wall 300, top wall 400 degree centigrade. So, you know the initial temperature distribution. So, you know the initial picture of let us say if you want to create a movie of temperature variation. Now, when you use this temperature in information, initial and boundary condition information, there are two approaches. One if let us suppose in this case I use a temperature of the previous time level. So, in this how many unknowns I have? So, all the temperature values which have superscript n I know from the previous time level. When I am just starting I know from initial condition. So, then in this how many unknowns are there? 1 sorry which is this one. Let us suppose in a problem there at 25 yellow circle you will get 25 equations, but in each equation there will be only one unknown. This is easier to solve because each equation can be solved independently. This is explicit method. So, I am writing separately this is the unknown which you want to do and these are the neighbors. I will show later on an animation stencil. So, you get here the equation is simple, but there is a price which you have to pay for. I will discuss that little later. There is a second method where this total heat gained by conduction is calculated using temperature of the same time level. Now, the problem here is let us suppose if you have 25 points you have 25 equation, but in each equation there will be not more than how many unknowns? How many neighbor comes into the expression? 4 neighbor and the center point and if all the neighbors are yellow circle and the center point is also a yellow circle. Then you will have maximum number of unknowns which is 5 unknowns, but if some neighbor is blue circle it is a boundary point and that point let us suppose it is a boundary condition then there will be less number of unknowns. So, in an implicit method in a two dimensional problem you get set of equations which are coupled. You cannot solve each equation independently. So, this has to be solved in CFD one way of mathematically one way of solving is that you take the inverse of the coefficient matrix multiplied by the constant vector, but we do not follow that method is called as direct method of solving matrices. Here what we follow is that we solve it iteratively. So, we start with the initial guess of the temperature at new time level of the neighbors and then keep iterating till it converges. For each time level there is a loop for when you do an unsteady simulation there is one loop for time that will anyway be there, but in an explicit method you get let us say 25 equation in each equation there is only one unknown. So, you can solve each equation independently, but in an implicit method at each time step you have to solve set of algebraic equation iteratively. So, there is an inner loop involved when you use an explicit method there is a topic called as stability analysis. If you look into book of JD Anderson he has shown the stability analysis using von Neumann stability criteria this what is mean by stability? Stability means that when you are solving for let us say steady state problem the numbers what do you expect in the real world situation the temperature should the difference of the temperature between two consecutive time step should die down if the temperature should reach to steady state, but computationally if your time step is not small enough depending upon grid size and thermo physical property. So, in an explicit method you need a certain you cannot have time step greater than sorry greater than certain value. So, what I mean is that let us suppose you have particular plate made of particular material particular grid size then let us suppose using stability criteria you found that the minimum time step should be let us say 1 second then you need to make sure that if you use a time step greater than 1 second then what will happen is that your code will not converge physically you want that the temperature difference should die down, but numerically it will happen in this method that it will not die down, but it will grow. So, the difference of the temperature between two consecutive time step approaches towards infinity. So, your numerics is not behaving what it should behave physically. So, from stability criteria one expression is derived alpha delta t 1 by delta x square plus 1 by delta y square should be less than half. So, this delta what is alpha k by rho c p this delta t now becomes function of thermo physical property and delta x and delta y. So, there is a time constraint you cannot use a time you have to use a time less than a particular value there is a limit which you can use like if you want let us suppose in a real world situation you took that plate from furnace and it reached to steady state in let us say 1 hour. If you are using an implicit method you can take a time step of 1 minute and you can do 60 time marching and reach to the steady state, but let us suppose you are using an explicit method and it demands that you have time step less than equal to 1 second and if you use 1 second as the time step you have to do 3600 times of time margin. In implicit method there is no such time constraint whether you are able to understand and appreciate. So, let us suppose you have a plate which reaches to steady state in 1 hour and if the for the same plate you are doing a simulation in an implicit method you can use a large time step, but it is not that you can use 1 hour directly as time step and then from 0 initial condition 1 time I think you should reach to the steady state. You cannot use a very large time step in implicit method, but you can use a reasonably large time step may be 1 minute or 5 minutes, but in an explicit method demands that let us suppose if delta t should be less than 1 second then you have to ensure that it is less than 1 second. But the problem is although let us say explicit method takes 3600 time marching whereas let us suppose implicit take 60 time marching, but for time marching which takes more computational time implicit because explicit you get 25 equation in each equation there is only 1 unknown you can solve it straight forward, but in implicit you have 25 equation with each equation not more than 5 and 1. You have to solve it iteratively per iteration time computational time which is required for implicit method for time step computational time required with implicit method is more. So, let us suppose per time step explicit method takes 1 second and implicit method takes 1 minute then what will happen there will be a break even both will take same computational time. So, it is difficult to say which one is better, but I will discuss the advantages and disadvantages. I will just like to show you some animation of the how the numbers evolve in an explicit method. This is a type of instantial in an explicit method. So, to calculate the temperature at new time level this is function of 5 different values of the previous time step. So, if this point is p then the same point at previous time level and the neighbors. So, let us suppose this is a plate which you have taken from furnace all the yellow circles are at 250 degree centigrade. This blue circles are at 100, this bottom blue circles are at 200, this is at 300, this is at 400. So, this is the initial temperature distribution which you know. Now, the way we have done finite volume discretization that is suppose there is no volumetric heat generation. So, this value in explicit method we get an expression like this where there is only one unknown. For the new time level there is only one point which is coming in this tensile and this is function of the same point it is not south east western 5 point values of the previous time step used to calculate the temperature of the new time level. This you do for, but in this case as I said if there are boundary values also which are used just see that. So, this equation is used how many times? There are 9 yellow circles it is used 9 times, but what will happen only the neighbor should change. So, in programming you can do it very easily how you do? Using loop. So, you vary i from let us say 2 to 5 j from 2 to 5 2 comma 2 it will solve here 2 comma 3 it will go here. So, this animation which has shown this is for 2 comma 2 2 comma 3 and so on. So, this way I am showing you that this set of algebraic equation in a loop how it evolves with respect to time 5 values are picked up from previous time level and the new time level value is obtained. This happened behind your computer screen after this what happens? So, let us suppose from 0 second you have obtained the temperature distribution for 1 second. Then from using the temperature of 1 second you go to the 2 second using the 2 second value you go to 3 second. So, you have pictures for 0 second 1 second 2 second 3 second till 3600 second. Then what you can do? You can club all this picture in a sequence and create a movie. This I showed for heat transfer the same thing happens for fluid flow also. In an implicit method what happens? You have a stencil like this what you are finding here is that from the previous time level only 1 value is coming that is the same point value. And in this the stencil is such that there are 2 bonded there you can see that in this stencil 2 blue circles are coming this is 100 sorry. Yeah boundary condition does not changes with respect to time. So, this blue point values will not change with respect to time. So, as I said that this boundary condition is such that this left wall is at 100 and this is at 200. So, 2 values are known from out of this 5 values. So, there are 3 unknowns here 3 yellow circles. If you go to the next point 4 yellow circle 4 unknowns 3 yellow circle 3 unknowns. So, you get self-algebraic equations you have not more than this here you have 5 unknowns. So, you get some equations where are there are 3 unknowns 4 unknowns. Basically the corner yellow circles have 3 unknowns central yellow circle has 5 unknown corner yellow circle will have sorry this yellow circle has 4 unknowns this is a corner 1 it has 3 unknowns. So, you get here set of a yellow equation where there are in this point you have maximum unknown 5 unknown where all neighbors are yellow circle. So, this has to be solved iteratively. So, other than the time marching there is a second loop to solve the set of a yellow equation. Any question into this? Is it occur only in the explicit method? Yeah. So, the question is the stability occurs only in the explicit method the answer is yes and first clarify more I think you need more understanding into what stability is what happens is that there are certain truncation error which occur whenever you solve in computer there is certain precision with which you solve. So, if those truncation error grow with respect to time then your numbers approaches towards infinity. So, with respect to time how this numbers grow if you follow a book on stability analysis you get all these details. So, from there this type of criteria are obtained this type of stability criteria you can mathematically it has been done for pure diffusion equation pure advection equation. But when you go to Ney-Stokes equation it is a combined convection diffusion equation with a source term. So, for that there is no stability criteria. So, when we use the Ney-Stokes equation we do not have any stability criteria for explicit method in Ney-Stokes equation. But what we do is that we take the stability criteria from pure convection pure advection pure diffusion and we use that as a guideline at least we make sure that. So, this is a necessary condition but it is not sufficient for Ney-Stokes equation. Yeah, I will just add to that for your information. There is a very good discussion based on sort of physical arguments in Patankar's book if you want to get some sort of an idea of why a numerical scheme should behave whether it is stable or unstable. It is a very nice discussion and purely based on physical arguments and it is worthwhile reading that it is I think in the chapter number 3 and 4 initial when he talks about heat conduction itself. What Professor Sharma is mentioning is actually what we what is called as a mathematical exercise called the von Neumann stability analysis as he mentioned. Is this a false diffusion? No, no, it is just the maintaining of the positivity of coefficients. Okay, okay, correct. It is a very nicely written set of paragraphs and I think that should help you a little Yeah, in fact I am planning to discuss that when I teach advection the coefficients become peccalite number relationship there is there in the advection maybe I will it is good that you pointed out I will emphasize that. So, let us go to the advantages and disadvantage of the two approaches. The explicit method is relatively easy to set up and program and that is what we are giving you because we want to initially give you programs which are easier to understand. So, the codes which will be there in tomorrow's lap session or conduction are based on explicit method. The implicit method advantage is that it takes lesser step to reach when you want to do as I said if you want to do calculation for 1 hour you can take lesser number of time step you can take let us say 1 minute then you need only 60 time step. The disadvantage of this explicit method is that it may take it is not always necessary that it may result in long computational time make calculations over a given interval of time. So, let us suppose if you want to do steady state computation of 1 hour then if your time step restriction is 1 second then it may take more time, but it is not always true. For a time step calculation computation time is more in case of implicit method as I discussed earlier and for transient problem if you are solving steady state problem probably it makes much sense to use an implicit method, but if you are solving a transient problem where you want to capture true transient in fluid mechanics or heat transfer there are 2 class of problems. First class of problems which reaches to a steady state situation and there are second class of problem which reaches to maybe let us say a periodic flow like if you talk of flow across a circular cylinder there is one vortex shedding phenomena if you have heard of. So, what happens the flow does not never reaches to a steady state, but it reaches to a periodic state. So, with respect to time the vortices break away from the rear side of the cylinder alternatively from the top and bottom. So, for inherently transient problems if you want to capture the transients accurately it is very important that you should take then the time step play important role and then perhaps explicit method not always is more beneficial. And second thing is that nowadays we talk of parallel computing. So, when you go to develop parallel algorithms it is easier to do in explicit app using an explicit app. Now, let us go to the implementation details. So, we have discussed the theory now how to develop code. Now, to develop code I will show give you some figures pictures to give you an idea how to develop code in easier way this is the way the code has been developed which is given to you. This is the simplest form of the grid generation equispaced horizontal and vertical lines. What I do here for to make you explain this in a easier way is that let us say this is the grid point of temperature. Where do we have heat fluxes at the faces? Where do we have heat fluxes in the x direction? Which faces? Vertical faces. This is the boundary points for the temperature. These are the running indices in the x and y direction. So, I am defining green square here represent heat flux in the x direction. Red inverted triangle represent heat flux in the y direction. Now, what happens to a face? A face is common to two control volume, but we want to represent this heat fluxes only by one indices. So, we need to follow up convention. So, we say that now positive faces in a control volume what are positive faces? East and north are positive faces. So, there are you can what I am trying to emphasize here is that we have three variables here temperature q x q y. Now, I will store this in a matrix form. When I store in matrix form there I need a running indices to tag a particular value of the any of this variable. You can have three different running indices one for temperature one for q x one for q y, but that is not good way of programming. Good way of programming is that I use only one running indices for one of the variable and I follow a convention. So, the convention which I am following here is the any green square I have to see that it what is the running indices of the control volume in which it is in the east face. So, the running indices of this green square corresponds to the running indices of this in running indices is j comma i. Similarly, why the running indices of this is corresponds to j this is the convention I am following. If I ask you what is the running indices of the green square which is here this will be at the east face center of yellow circle which is here what is the what is the running indices of that yellow circle here i minus 1 comma j. So, the green square which is sitting here will have q x j comma i minus 1 what will be the if there is a red inverted triangle here also this will be at the north face of j minus 1. So, in a cell there will be two green square on the positive face it will have the same running indices on the negative face here it will be i minus 1 and here it will be j minus 1. So, these are the green squares which are coming now let me ask you if I talk of the points which are inside the yellow circles which are inside the plate how many are they 5 in the x direction 5 in the y direction or these are the green squares how many are they in the x direction in the y direction. Here at this basically green square they are vertical face center right now we are having 5 control volumes in the x direction 5 control volumes in the y direction. So, we have 25 control volume in this 25 control volume what I have shown as green square are the number of vertical face centers how many are they how many in the x direction 6 in the x direction you have to calculate heat flux at the boundary also note that because when you apply the law of conservation of energy you will need this heat flux here. So, how many green squares are there in the x direction 6 how many in the y direction 5 ok just note this because these are some of the numbers which we will be using it in the pseudo code in the next slide. How many red inverted triangles are there in the x direction 5 and how many in the y direction 6. So, there is a interchange green square is 6 in x direction 5 in y direction red inverted triangle is 5 in the x direction 6 in the y direction. How many total circles are there if you take the boundary points also 7 by 7 these are the 5 yellow circles which are seeing and one is at the left boundary one is at the right boundary where you apply the boundary condition. So, including boundary there are 7 by 7 yellow circles and the green squares are 6 in the x direction 5 in the y direction ok. Now, let me ask you a question that when I was showing the discretization I did the discretization first for analogous to del y del x of q x and del y del x of q y. So, I said that in the coding also for you to explain this in a easier way what I the code which I will show and what has been given to you first we calculate heat fluxes. So, to calculate heat fluxes I have to scroll through all this green square. Now, to scroll through green squares I should know what is the running indices what is the running indices of you can know the running indices of any if I want to know the running indices of this green square I need to know the running indices of which yellow circle this yellow circle. What is the running indices of this yellow circle this is 1 this is 2 this is 3 this is j is equals to 3 and i is equals to 1 2 3 4. So, the q x here is i is equals to 4 j is equals to 3 q y here is same i is equals to 4 j is equals to 3 and here it will be i minus 1 here it will be j minus 1 and so on. Now, when I want to scroll through all this green squares this is my first point this corresponds to which control volume this is at the boundary. The way we are applying boundary condition the volume of the boundary control volume is 0 you can assume that this point is centroid of a volume which is having a 0 volume is this ok. The way we are applying the boundary condition we have this is a 0 control volume that is why the centroid of that control volume in the phase center is coinciding this is only possible if the volume of that control volume is 0. So, it is there is only one point there is no if the volume is 0 there will be no east point west point there will be only one point and the running in this is for this point will be i is equals to 1 j is equals to 1 what will be for this it will corresponds to this is i is equals to 2 j is equals to 1 ok. So, i is varying from 1 2 3 4 5 and 6 what is your i max i max here is 7 let me show you ok the way I am showing i max is this is my i is equals to 1 this is my i is equals to 2 3 4 5 6 and 7. So, my running in this is for temperature let me show you for temperature i is equals to 1 and i is equals to 7 is boundary j is equals to 1 and j is equals to 7 is boundary interior points are 2 to 6 in both the direction when I want to calculate q x it is varying from 1 to 6 this is the i for this will be i for this i for this is 6. So, it will vary from i is equals to 1 to 6 or i is equals to 1 to i max minus 1 and in the y direction it will vary from j for this will be equal to 2 2 3 4 5 6. So, in the i direction it varies from 1 to i max minus 1 and in the j direction it varies from 2 to j max minus 1 ok and i max j max are 7. So, it becomes 1 to 6 and 2 to 5 and if you go to red inverted triangle they are 5 in the x direction 6 in the y direction. So, here it varies from 2 to i max minus 1 and 1 to j max minus 1 if you understand this clearly through this picture you can develop a code. So, this is the pseudo code to calculate heat flux in at the x direction i max is 7. So, it is from 1 to 6 j max is also 7 this is 2 to 5. So, this way you scroll through all the green squares this is a to calculate the distance between two cell center what is this x l x l is the coordinate of the cell center. So, this is the distance between the cell center this is the previous time level temperature between the two neighbors this is at the east phase center ok. I will show in the earlier slide this q x how you can calculate temperature of i plus i is in this direction do not get confused. So, this will be p at i plus 1 minus p at i and how you will calculate q i j is in vertical direction note that temperature at j plus 1 minus temperature of j. You will not only use temperature, but in the denominator you will need to have the coordinates of the cell center. This will be x cell center and this will be y cell center. So, if you do use that you can calculate this is the heat flux in the x direction this is resistance between the cell center in the x direction. And this is note that here it is j is varying from 1 to 6 and i is varying from 2 to 5. So, this way we calculate the q x and q i on all green circle and red inverted triangles respectively. So, this completes the you know first aid of differential equation is del vadal x of q x and del vadal y of q y. In algebraic what is the discrete representation summation on west and south phases q into surface area minus summation of heat fluxes into surface area on the north and the east phases. This is the code of that. So, with this you have calculated the heat fluxes. What do you do in sorry this is basically the discrete representation of Fourier law of heat conduction. Once you know this what you do you do a balance. So, for a control volume if you want to do a balance I had told you that on the west and south phases negative phases heat is going in and in the positive phases it is coming out. So, what will be expression for balance q x here it will be i minus 1. This will be q y j minus 1 you add this 2 because here it is coming in you want to get gain of the heat transfer and this 2 you subtract multiplying by surface area. This is what is done here for the x direction surface area is delta y for the y direction surface area is delta. So, this is just the balance and also note that what is the loop here. This you are where you are calculating the balance at yellow circles. This loop is corresponds to green squares. This loop corresponds to red inverted triangle and this corresponds to yellow circle. So, note the loops which you are using. So, you calculate the appropriate variables at appropriate locations. So, you balance it once you have balance it got the total heat gain by conduction. Then you multiply by delta t rho c p v. This is the expression which you have got from the explicit method. Any question into this? In the program you will see let us say in psi lab language which is basically a c type of language. Same loop will be there. So, whatever is shown here you will can correlate directly in the code go back and see once we give the code as well as the slides.