 So, after the implementation details I will discuss the solution algorithm. So, the solution algorithm is as I said that the first part of any code is preprocessor and in preprocessor basically what we do is that there are certain user inputs like in the codes which we had given you will have to enter what is the density of the material specific heat of the material conductivity of the material. So, the material property is one of the input thermo physical property in general. Then there will be some geometrical parameter what is the length of the plate height of the plate what is the maximum point in the x and y direction what are the boundary conditions input. We have given you a code where which is made general like although it is in Cartesian coordinate system. So, when the code is in Cartesian coordinate system it is expected that the boundary will be having a left wall right wall bottom wall and top wall. However, there are switches like on the left wall whether you want Dirichlet boundary condition Neumann Neumann means insulated or convective boundary condition. Actually we have given four options if I remember constant wall temperature insulated boundary condition constant heat flux boundary condition and convective boundary condition because these are the four boundary conditions which you see in heat transfer. So, there are certain boundary conditions input corresponding to each wall and then there is a convergence criteria to check whether your code has reached to steady state. These are some of the inputs then there is a grid generation in this case we are following a simplest form of grid generation where we need to calculate the coordinates of the vertices of equispaced horizontal and intersection of the equispaced horizontal and vertical lines. But when you generate grids you calculate the width of the control volume distance between the cell centers surface areas and the volume. Then you set the initial guess the code which we are discussing here is an unsteady state code. So, for this you need an initial condition the problems which will be given you there will be an initial condition also. Note that here I am setting that the here you set the Dirichlet boundary condition I am not talking of non-Dirichlet boundary condition. Non-Dirichlet boundary condition I have put in separate step why? Note that the boundary condition I have I am in step 3 I am writing only Dirichlet boundary condition. In step 4 I am writing non-Dirichlet boundary condition why? In Dirichlet and non-Dirichlet what happens the boundary values what happens it changes with respect to time it is a function of border point. Boundary point is a function of border point whenever you have a gradient boundary condition and you know border point is a part of inside the domain it changes with respect to time. So, whenever you have a non-Dirichlet boundary condition what a non-Dirichlet Neumann or Robin type of which I discussed earlier for those type of boundary condition you have a gradient and you get when you discretize it you get boundary point as a function of border point which changes with respect to time. So, this has to be brought inside the loop after each time marking you have to update this boundary point value as a function of border point. Before you go to computation you whatever temperature you have now take it as a old value then I had done input step for each phases calculate q x and q y then do a balance and calculate total heat gain by conduction and then finally temperature. Once you have got temperature then you check for convergence when there are two ways of checking the convergence in fact there are many ways, but one of the common ways is that one of the simplest way is that between two consecutive time step at each point let us say there are 25 yellow circle. So, you will have 25 temperature from the previous time step and from the new time step you take calculate the 25 differences of those temperature and out of those 25 difference values what whichever one is the maximum you use a small program to obtain the maximum of a matrix that maximum you check for steady state criteria, but a better way it has been found that instead of taking the absolute difference if you take the RMS of the difference. So, you take the difference between the present and old time take it square at all the so at 25 points you take the difference take it square take mean of that and take the square root of that. So, this instead of taking the absolute different absolute maximum difference this is considered to be better. So, we take the root mean square value of the temperature difference and that is what in your code you see RMS it has been defined given a name as RMS residue if I am not wrong. So, this is used to check the unsteadiness, but many a time now many a time the question comes whether we should call delta t is it tending to 0 as steady state or whether we should have del capital T by del small t tending to 0 as steady state. So, if you ensure this you ensure try to ensure that delta t tends to 0, but if you divide by delta t time step you try to ensure del capital T by del small t tends to 0. However, if you are using explicit method and use this criteria for conduction this works properly there is not much problem, but when especially when you are solving we had seen that if you solve a natural convection problem and if the temperatures are very large then the velocity inside in natural convection increases tremendously then the time step varies a lot. So, what happens is this temperature difference between 2 let us say consecutive time step goes down, but if we are by stability criteria if your delta t is increasing then this value oscillates a lot. So, for I suggest that it is better to use this boundary condition because sometimes this cause lot of stability problem. So, in the code finally we have used this boundary condition, but I just wanted to mention here if it has converge means it has reached a steady state if not you go back to there is a mistake here it has to go back to step 4 not step 5, because in step 4 you have non Dirichlet boundary condition you have to go back to step 4 where you have a non Dirichlet boundary condition where boundary is a function of border value. So, this loop goes on and on till you reach to a steady state till this r m s t is less than equal to epsilon which is a convergence criteria may be 10 to power minus 3. So, this completes the formulation the implementation detail and solution algorithm any question on there it is very important that you understand this. I will repeat again in the formulation we have used two levels of approximation let us say in case of diffusion the first level of approximation was whenever you take a control volume to calculate total heat gain by conduction you express heat transfer as fluxes on the phases. So, the first level of approximation corresponds to how to represent the heat fluxes at the phases how to write q x q x plus delta x q y q y plus delta y on the surfaces as one point value because q x varies all over the vertical surface q y varies over all over the horizontal surface. So, this approximation is first level of approximation what is the second level of approximation how to write Fourier law of heat conduction in algebraic form how to convert the differential form of Fourier law of heat conduction into algebraic form. We basically use second order central difference scheme piecewise linear approximation to represent heat flux at phases we use centroid value which is also second order approximation these two approximations we use and we get algebraic equation. Then I talked about the solution methodology explicit method and implicit method then I talked about how when you want to develop the code let us do in stages because formulation also we did in stages let us develop code in stages. Let us first calculate q x q y for that we use a discrete form of Fourier law of heat conduction for that we use some convention add even shown some symbol some colorful squares and triangles. So, that you get a picture of in your mind it is not mathematics is some picture which you see and you know what physically you want to do you want to balance the heat and then once you have balance the heat to know the temperature evolution you basically apply rate of change of internal energy is equal to total heat gained by conduction and then finally I had shown you step by step procedure as a solution. Any question please. So, let us suppose there is a plate of length l 1 in the x direction l 2 in the y direction I will take two class of heat transfer problem one where I will take all Dirichlet boundary condition. So, all Dirichlet means I take let us say left wall 100 bottom wall 200 right wall 300 top wall 400 degree centigrade and for this if you use the formulation the implementation details and if you develop the code using the solution algorithm finally you create a movie where at each time step if you have stored the temperature for each time step you get one picture. So, for this case if you take a steel plate. So, in your code you there will be density specific heat and conductivity as a property and this is the way the initial condition is 100 degree centigrade in this problem this is isotherm note that left wall is always blue 100 degree centigrade bottom wall is always green 200 right wall orange 300 top wall close to white 400 boundary condition should not change with respect to time this is just one of the check let me tell you in research we do not know what results we will get, but we should be clever enough to know what results we should not get. So, if you develop a code and there where you are finding that your left wall is not 100 bottom wall is not 100. So, it is not obeying the boundary condition how can you expect that your results will be correct. So, the purpose one of the implicit purpose of this course is that to make you clever enough that when you use CFD you know that at least this result I should not get or if your student comes with some result you should say no you are wrong just by looking into this not going into the details of the result just looking into the result you can understand I will give you some more hint to know how good or bad a result is or whether it is a physical or it is not ok. Now, whenever many times you may feel that I can use a software and get engineering parameters and let us suppose for flow across a car or flow across an aeroplane I will calculate lift force and drag force, but what is the need of that creating that movie. Let us say for flow across a car fluid dynamic movie there is a need if you have a if you want to have if you want to understand the fluid dynamics phenomena it is a I would say that is more of a scientific investigation if you want to understand that if you have designed 10 different shapes of the car and let us suppose you do CFD simulation and you find that one of the design is giving least drag force, but the reason why that design is giving least drag force if you create a fluid dynamic movie and not only velocity contours or let us say velocity vectors or pressure contours you may have to draw a vorticity contour you may have to develop some other tool to understand because fluid dynamics is not that easy that you create one type of movie and you can understand everything. Many a times we do not as I mentioned say which cannot say with certainty that drag force is reducing due to this there are many reasons because of highly coupled problem. So similarly here I want to give you a message that once you get this type of temperature distribution it has some story what is that story bottom wall is at 200 degree centigrade can you talk about the rate of heat transfer in the bottom wall this is note that this is a steady state heat conduction with no volumetric heat generation. So, the solution has to be bounded and there cannot be local maxima in the temperature distribution what I mean is that let us suppose if this is at 200 this is at 250 so all the point between in this region will be in between 200 and 250 it cannot be less than 200 or greater than 250 there cannot be a local minima or a local maxima if you have volumetric heat generation then it may happen that you can have a local maxima, but there also it is not necessary it depends upon the magnitude of volumetric heat generation as well as the boundary condition. So, near to the bottom wall by looking into this temperature distribution bottom wall is at 200 degree centigrade can you make a statement on the direction of heat transfer at what which point heat transfer will be 0 on the bottom wall bottom wall is at 200 there is an isotherm 200 which is heating it. So, at this point rate of heat transfer will be 0 it is heating normally d t by d y will be 0 what about left part of this fall temperatures are what 150 160 170 180 190 200. So, temperature is here bottom wall 200 in this region temperature is less than 200. So, heat transfer will be upward here it will be here it is greater than 200 210 220 230 240 and so on here heat transfer will be downward. So, this picture they speak something ok. So, many times when you get if you plot let us suppose in this wall how q y varies with respect to x. So, this picture says that it should change this direction it should pass to 0 initially it should be upward and then it should go to downward crossing 0 smoothly. So, this I give wanted to give an idea that when you plot the engineering parameter such as in this case I have taken a simplest of the problem the variation of q y as a function of y this is called as local heat fluxes this local values when you do post processing in the next slide may be I will show you. So, on the top wall you see everywhere the temperature is less than 400. So, heat transfer is inward on the right wall what you see is left wall is 100 and near everywhere it is greater than 100. So, heat transfer is outward on the bottom wall there is change in the direction of the heat transfer and in the right wall also there is change in the direction of heat transfer ok. So, this picture gives the reason why there is change in the direction of heat flux on the right wall and bottom wall. Now, how do we calculate because once you get a run a CFD code you get a temperature distribution, but how do you calculate this heat fluxes at the boundary actually the way I had given you the formulation the codes we are calculating q x q y. So, it is automatically calculated in this formulation ok, but if you are solving let us suppose by finite difference method where you are obtaining temperature distribution. So, let us suppose your coding does not involve q x q y because you can develop code in that way also and if you want to calculate q because this q x q y strictly speaking as an engineer you need it only near to the boundary not interior of the domain. Many times you may feel that it is unnecessary, but I had shown you because I felt that this is easy way to understand when you want to develop into code. So, normally if temperature distribution is known then at this boundary what we do is that how do you calculate q y here you will have to calculate q y. So, q y is equal to minus k del t by del y. So, you need to calculate del t by del y at this boundary point using the border point note that here you have to use one sided difference like forward difference. So, let us suppose at the bottom wall left wall it is 100 bottom wall it is 200. So, if you use d t by d y this temperature minus this temperature divided by the distance between the two points you can calculate del t by del y multiplied by minus k and then see. So, what I wanted to show you here is that you can see that in this three values are less than 200 here the heat transfer will be 0 here that values are greater than 200. So, there is a change in the magnitude of temperature from the left part to the right part across 200. So, if you plot at the bottom wall which is this I think red inverted triangle what you see is that initially it is positive then it crosses the 0 then it becomes negative on the bottom wall. So, note that we do numerical differentiation to calculate the local value of heat fluxes. So, in the bottom wall the heat flux let us suppose you have five points. So, you will get five different value, but many times we are interested in one value then you have to do numerical integration. So, numerical differentiation is needed to calculate the local heat flux and you need to do numerical integration to calculate the total heat fluxes. Many times if you are teaching a course on numerical methods you can take this as an example for numerical differentiation and numerical integration. These are good examples from heat transfer. We do this in our first course which we call as the computational method in thermal and fluid engineering. This is a second class of boundary condition. Left wall now I have all types of boundary condition Dirichlet, Neumann and Robin-Horwitz type of boundary condition. Left wall 100, bottom wall insulated, right wall constant heat flux top wall convective boundary conditions. Now, in this case when it is a Dirichlet boundary condition it is easy to apply the boundary condition because it is just the value which you have to give to the blue circles, but in this case you need to discretize like this is insulated del t by del y is 0. So, the temperature here is equal to this, the temperature here is equals to this. This is written as p at all i j is equals to 1 is equals to p at all i j is equals to 2. This is the discretization of constant heat flux boundary condition. This is the discretization of Robin-Horwitz mixed type of boundary condition. Minus k d t by d y is equals to h temperature at the wall minus t infinity and this is the Dirichlet boundary condition. Note that whether you are applying finite volume method, but the discretization of the boundary condition is done by finite difference. That is why some amount of finite difference was very necessary before going to the finite volume method. What I would do in a first we calculate the temperature. Once it is convergent then once only once we calculate heat flux values. What is the need not to calculate every time? There is no need. As I said that there is no need. It is just that I feel that this is when you are coding you have the feel of what you are doing, but you can avoid it. It is not necessary. It is just to have a better understanding of the conservation laws and the formulation and the coding. I wanted to have a same relationship between all of them. That is why I followed this example. For a beginner it is good, but later on it is not necessary. It takes unnecessary computation time. Now for this problem let us suppose left wall is at 100 degree centigrade, bottom wall insulated. Right top wall let us suppose there is a force convection with h is equals to 100 watt per meter square Kelvin and the ambient temperature is 30 degree centigrade. Right wall suppose let us say wall heat flux is 10 kilowatt per meter square. So for this type of animations which you will get for temperature distribution. Now here again you by looking into this animation you can understand whether it is correct or not. How? You should know what is your boundary condition and by looking into each of this picture you can say it is correct or not. Because from the boundary condition you know the slope of the line. Bottom wall is insulated. So the isotherm should be in heating perpendicular. See this animation. In each picture just see at the wall not away from the wall. At the wall it is heating normally first thing. Second thing right wall you have a constant heat flux. So this will not heat normal that is okay. But what is dT by dx equals to minus qw by k. So the slope has to be constant. So note that in this animation on the right wall each isotherm is heating at same constant angle. So whenever you are doing CFD simulation whenever you get some result be clever enough to judge from the boundary conditions or the flow structures near to the boundary. That whether your reserve seems to be physical or it is just a numerical artifact. This is a steady state temperature distribution. Just I would like to point out for a beginner as far as this graphics is concerned. Once we get the temperature distribution you get in a matrix form and then you take it to a software. And what you need to do is that you need to give the x coordinate of the yellow circle points or blue circle points, the y coordinate of that and the temperature value. And then you can use a contour command and then you can get and there are two types of contours. Earlier what animation which I was showing it is a flooded contour and what is shown here is the line contours. Here you have an isotherms which are represented by line and in the previous animation it is a flooded contour. So although there is a line also there is a combination there is an overlap plot of the flooded contour and the line contour. So this contour here represents the isotherms. This is the steady state variation. Now here again the same thing happens left wall is at 100. There is a top half where the temperature nearby is less than 100. Bottom half which is greater. So there is change in the direction of heat flux on the left wall. Bottom wall is insulated so no heat transfer. Top wall you have a convective boundary condition. So there is a heat loss from the top wall. Right wall you have a constant heat flux. So it will have inverted heat flux. Here again if you plot here I have plotted the heat flux on the left wall. So you can see that on the bottom half bottom portion of the left wall the heat transfer is negative. It crosses through a zero and then it becomes positive. There is a change in the direction of heat flux on the left wall. So that completes the conduction in Cartian coordinate system. I was not planning in cylindrical coordinate system but I think in the Boranich lecture somebody mentioned that if we can incorporate cylindrical coordinate system also included but the idea which I had used in Cartian coordinate same idea is used here. It is just that the geometrical parameters which are involved are different. And moreover let us suppose if you are taking in a polar coordinate system use the same law of conservation of energy rate of change of internal energy equals to heat gain by conduction plus heat generation which you can represent as this. So there I was showing you qx, qy. Now if you take a polar control volume what happens is that you get heat transfer in the radial direction and angular direction. And in this case the surface area are different. The surface area here is in a 2D case delta r on both the surface and here it will be r plus delta r into delta phi this will be r delta phi. So the surface area in this case what you find is that the surface area changes along the radial direction. It does not happens in case of x and y direction. Other than that the way we do the discretization it is almost same. So if you want to calculate this there is one yellow circle lying here one yellow circle lying here. So all those ideas we use actually right now in this slide I think what I am trying to show you is that you can get the differential equation. How you get the differential equation in your heat transfer course and cylindrical coordinate system? You do a balance and divide by the volume use the limits. So this is the differential equation in cylindrical coordinate system. In the next slide I show you the algebraic form. This is the first level of approximation where what you do is that at all the four surfaces you represent the heat fluxes at the centroid of each of the surfaces. So you have centroid represented as east, west, north, south. So instead of writing q phi q phi plus delta phi q r q r plus delta r you write the values at centroid as subscript and then you calculate do the balance. So from the negative faces it is coming in from the positive faces it is going out. So in minus out is the total heat gain and the second level of approximation is you do the discretization. The Fourier law of heat conduction is slightly different in cylindrical coordinate system as compared to Cartesian coordinate system. So this is the discrete representation of Fourier law of heat conduction in cylindrical coordinate system. Delta phi is the angle between these two points. Rp is the radius at point p. These are the geometrical parameters which are involved. So with this you can calculate that this qe, qw, qn and qs and then you do a balance and then the procedure is very much similar to what I have shown in the conduction. Even you can develop a general purpose code both for Cartesian and cylindrical coordinate system and only what will vary is the geometrical parameters. Like the denominator may be in a Cartesian coordinate system this was delta x here it is rp delta phi. The geometrical parameters will only vary and this is the explicit when you use an explicit method this is the final. So from the second level of approximation if you substitute this q into the expression which you have got in the first level of approximation you get the final equation as this and as we are calculating the total heat gain by conduction using the temperature of the previous time level this is an explicit method. So for explicit method this is the equation discretized equation. I will not go into much detail in cylindrical coordinate system but I would like to mention here that in a cylindrical coordinate system you have boundaries let us say in this problem at least polar coordinate system let us suppose you have a plate with a circular hole circular plate with a circular hole. So there is an inner boundary and there is an outer boundary both are circle r is equals to r1 r is equals to r2. So in this case what you find is that you do not have a computational boundary in the angular direction. So if you are solving for because this is circular so what happens is that if you want to solve for this this is one of the neighbor this is another neighbor actually when we want to develop a code even if it is in a cylindrical coordinate system we want to map this point in Cartesian coordinate system like this because we need to scroll in some starting from some point I need to start my i from some point. Let us say I start my j from the center in this direction let us suppose my i is starting in the anticlockwise direction from here this is my i is equals to 1, 2, 3, 4, 5, 6, 7, 8 but when I am solving for let us say this point I need one neighbor behind. So what I do is that I do not put this as 1 I put this as 2 and the value of this one will be equal to the last so what I want to tell you is that what is blinking here as this green line actually I say that this is a 0 volume control volume here I have some actually if you want to draw the domain this when it has to pick its east boundary it has to pick this value. So for the i is equals to 1 what I do is that for all j i is equals to 1 when we apply the boundary condition it is what is called as periodic boundary condition. So the boundary value is equals to the last border value so the first boundary value left boundary value is the so what is the last border value this one if you move in this direction which is the last yellow circle which you will encounter this one. So the way we apply boundary condition is that the value at i is equals to 1 is equals to the last yellow circle and the boundary condition at i max is equals to the first yellow circle this is what is called as periodic boundary condition in CFD. So we need to make a branch cut so that you can map it like a cartesian because finally when you want to scroll you want to scroll in this fashion where at the i is equals to 1 and i is equals to i max you use a periodic boundary condition here again I had used the same instead of q x I have q in the theta direction and q in the radial direction instead of q y. Here again I am showing you a pseudo code but this pseudo code will be slightly different from the earlier pseudo code because of the periodicity. You need not calculate q x on that branch cut that is why there is slight difference you can go into the details little later but what I maybe I can go to the figure and point out in this case I need heat fluxes on this phase as a function of these two values and this I need only once. So there is a periodicity so that I can avoid using the type of loops which you see in the theta and the r direction are little different because in the theta direction it is starting from 2 but in the j direction j is starting from 1 but so in earlier case this q y j was starting from 1 but in q x I was starting from 1 but in this case there is a slight difference due to that periodicity. As far as calculated this is calculation of q in the theta direction this is the calculation of q in the radial direction and then once you have calculated q in the theta and the radial direction you do a balance and then finally calculate the temperature of the Newton levels. So I think with this I had come to the end of this lecture on heat conduction.