 Let us quickly go through what has been discussed before the T break. I had started with the computational heat conduction topic with grid generation and showed that in a simple example a square plate I have shown you that the grid generation simplest form of the grid generation can be done by drawing equi-space horizontal and vertical lines getting certain equal size of the control volume. Then the grid points at the centroids of the control volume are shown by yellow circles and those are the points for which we want to get a linear algebraic equations. I had also shown you the boundary grid point shown by blue circles. So and thereafter there are not only do this geometrical division process but we calculate certain values also which are used which comes as a coefficient in the linear algebraic equations which we create later on. Those are the geometrical parameters like width of the control volume, the surface area of the control volume, the volume of the control volume, the distance of a cell center from a neighboring cell center. So this completed the grid generation part then went into the second level of discretization which is obtaining the algebraic equations from the conservation law. So here I had mentioned that here the procedure which is being followed is that we do not start from the differential equation in a CFD course. We start from a same control volume which is used in the undergraduate heat transfer of fluid mechanics course and obtain the same algebraic equations. So in the towards the end of the last lecture I had shown that term by term starting with the unsteady term then moving on to diffusion term and source term. For each of this term I had shown that we use two levels of discretization. The first level for the unsteady term is the volume averaging. The first level of discretization for the diffusion term is the surface averaging. The second level of discretization is the discrete representation of rate of change of temperature for unsteady term. The second level of discretization for the diffusion term is discrete representation of Fourier law of heat conduction. For the discrete representation of the temporal gradient we have used forward difference, first order forward difference. For the discrete representation of Fourier law of heat conduction we have used second order central difference or piecewise linear approximation. I had also shown you that the approximation here used are mathematically very simple which you have studied since your school days. So I believe that this procedure of obtaining the set of algebraic equation is not does not involve complicated mathematics like Gauss divergence theorem. And so the traditional method of finite volume method where we start with the differential equation. So when you teach CFD course you have to start the first or the second lecture with the set of differential equation. And the objective of here showing you the traditional finite volume method which you see most of the CFD book present CFD book is to show that the algebraic equation which I had shown you here. So in the next slide using starting from the governing partial differential equation I want to show you that indeed we get the same algebraic equation. So the end result is same whether you follow this mathematics based approach or the physics based approach which was which is proposed here. So in a differential approach or traditional finite volume method we start with differential equation. In the present method earlier I was starting with the control volume here we start with the differential equation. So this is the unstudied term this is the diffusion term and this is the source term. This equations are basically obtained from law of conservation of energy and applying delta volume tends volume of the control volume tending to 0. So there we have assumed that volume is tending to 0. Now in the finite volume method as here volume is tending to 0 as this equation is per unit volume basis we take the volume integral to convert this per unit volume expression into an expression which corresponds directly to the law of conservation of energy when the volume integral is taken there are certain terms which are converted into surface integral. Now you should understand that which are the terms which are converted into surface integral. So there is a mathematical procedure called as Gauss divergence theorem using which let us suppose I will show you that this volume integral k del square t which is represent total heat gained by conduction this volume integral is converted into a surface integral. Now in a Gauss divergence theorem when you apply the Gauss divergence theorem this del square t gets converted since double derivative gets converted into first derivative. However you have to take the product of this del t dot n. Now this del t dot n I would like to mention that what is this n? This n is a unit normal vector I will show you through a in whiteboard what is del t dot n. So let us take a control volume this control volume has four phases centroid of the east phase this is a centroid of the north phase west phase and south phase. Now what is the unit normal vector let us obtain what is the unit normal vector on each surface. In this surface unit normal vector is plus i in this surface it is plus j in this surface unit normal vector is minus j unit normal vector here is minus i. So on the positive surface is the unit normal have positive sign on the negative surface it is a negative sign because it is in the negative x and negative y direction. So when you want to calculate del t dot n what is del t? Del t is del t by del x i plus del t by del y j dot n. So when you do that what happens is you get for the positive phases you get a positive sign for the negative phases you get a negative sign. So it comes out that for this n is positive plus i for east phase for west phase it is minus i. So this expression I would like to highlight that the this surface area which we use we use the magnitude as delta x y for horizontal surface delta y for vertical surface and for the positive surface you have a positive sign for negative surface you have a negative sign and this way also finally you end up with the same algebraic equations for the diffusion term and here I would like to mention that unsteady term and the source term we do not apply the Gauss divergence theorem. So we apply the Gauss divergence theorem only to this term not to this term because this term unsteady term and the source term they are not flux term. So we apply this Gauss divergence theorem only to those terms in the conservation law which have flux term like for mass flow rate you have flux term mass flux similarly for enthalpy flux similarly for force per unit area. So those term in the conservation law which can be expressed as flux for those only we apply this Gauss divergence theorem the other terms are what you call as a volumetric terms like this is directly proportional to volume this volumetric heat generation is also directly proportional to volume. So we do not apply the Gauss divergence theorem to the unsteady term and the source term. So note that only those terms in the governing equation or conservation laws which consist of flux term are converted from volume to surface integral using Gauss divergence theorem in the traditional finite volume method. Thereafter the approximations which are discussed in the previous lecture the same approximations are used and we end up with same linear algebraic equation. So in both the methods the approximations are same surface averaging, volume averaging, use wise linear approximation, first order discretization in the time direction. So the approximations are same in the governing equation based traditional finite volume method and the control volume based novel finite volume method which is proposed here. This results in the same linear algebraic equation. However I had felt that when you want to teach specially under graduate students it is better rather than starting from this differential equations and reaching to this algebraic equation you start from a control volume which they see in their first class of let us say under graduate heat transfer course and reach to the same algebraic equation. So this way I had shown you that whether you follow the physics based or mathematics based finite volume method you get the same end product which is the algebraic equations. Now I will discuss that once you get an algebraic equation what is the procedure to solve the set of algebraic equations. Note that I am not if you look into the expression for temperature there is certain superscript let me show through here for temperature temperature terms are here and temperature terms are here. Now on the left hand side I am having superscript for temperature on the right hand side right now I am not showing you any superscript what does superscript represent? Super script represents time level. So this is the temperature of the new time level this is the temperature of the old time level. So what does this term represent is the total heat gained by conduction. Now if I use a superscript N in this case then I am using temperature of the old time level to obtain the total heat gained by conduction. If I use a superscript N plus 1 on all this temperature then I am calculating total heat gained by conduction using the temperature of the new time level. So in a solution methodology if you are using the temperature of the previous time level and then calculating total heat gained by conduction then this is called as an explicit method and in that case I just wanted to show you that the way we have discretized you calculate temperature at the new time level as a function of temperature at the old time level how many neighbors are involved in this expression east, west, north, south and temperature of the same point of the previous time level. So note that if there are 25 points if you are using this method you will get 25 equations and in each equation there will be only one unknown and that one unknown will use 5 values from the previous time level 4 from the neighbors and one from the same point. This is what is called as computational instantial later on I will show you an animation also for this. So this is the algebraic equation. So the nature of the algebraic equation is there are 25 yellow circles you get 25 equations in each equation there is only one unknown it is easy to solve this set of equations each equation can be solved independently they are not coupled because each equation there is only one unknown. Whereas if you use an implicit method the total heat gained by conduction is calculated using the temperature of the same time level in this case you there are 25 yellow circles you get 25 equations and in each equation there will be not more than 5 unknowns what I am why I am saying not more than because if it is a border cell then if there is some boundary conditions then one or more than one neighbor can be coming from boundary as a boundary condition if it is a constant temperature then that value will be act as a constant and it will come to be. So in this case if there are 20 yellow circles you get 25 equations 25 linear algebraic equations but this set of equations are coupled why because each equation cannot be solved independently because each equation has more than one unknown but the system is closed it is just that you have to solve system of linear algebraic equations which are coupled but note that all are linear equations. So you may be feeling that explicit method is good because in each time step if you are using an implicit method you get a you have to solve system of algebraic equations and in computational fluid dynamics how do you solve a system of algebraic equation one way you might have studied since your school days is that there is a coefficient matrix there is an unknown vector you take the inverse of the coefficient matrix multiplied by the on the right hand side the vector which you have and then you get the unknown vector. However in computational fluid dynamics we do not use that method which is called as a direct method of solving the system of algebraic equations because one of the main reason is as I said that there are 25 equations there in each equation there is not more than 5 unknowns. So if there are 25 so let us suppose t1 to t25 are the temperatures which you want to calculate out of the 25 temperature values in each equation there will be only 5 non-zero coefficients other 20 coefficients will be 0 in each row of the coefficient matrix. So this in this case it creates what we call as a pentadagonal matrix of the coefficient matrix so it is called as a sparse matrix because right now we are talking of 25 points in actual CFD calculations you can have millions of points. So when there are more number of points there are more number of zeros also and if there are more number of zeros taking the inverse of that coefficient matrix is a very costly affair. So we do not use that direct method of solving the system of algebraic equations we solve it iteratively using Gauss-Riedel method or Jacobi iterations or a dry diagonal matrix algorithm. So in each time step so as far as unsteady computation is concerned there will always be one time loop so your solution will mark picture by picture because let us suppose you want to create a movie of the temperature distribution. So you are to create each picture you have to move ahead in the time direction you start with one initial condition initial picture and you move in the time directions creating pictures for different time instants and finally you can create a movie out of it when you use an explicit method you do not there is no iteration involved because in each equation there is only one unknown however if you use an implicit method there is an interior inner loop also to obtain the temperature at each time step. So however it is easy to program the explicit method the explicit method takes less computational time per time step but there is a time step restrictions there is a stability analysis right now where I am not showing the detailed stability analysis if you want to refer you can refer a book on computational fluid dynamics by J. D. Anderson where it is shown that the stability analysis is used and the idea is to quantify how the error go with respect to time and this equation comes which is called as the equation for stability requirement. So the time step cannot be there is a restriction as far as time step in an explicit method is concerned. So once you have decided delta x and delta y let us suppose the example which we have taken delta x plus 0.2 delta y was 0.2. So in this case it will be 1 divided by 0.2 1 divided by 0.2 is how much 5 5 sorry it is a 0.2 square so it will be 0.04. So 1 divided by 0.04 is how much 25 25 plus 25 is 50. So 1 divided by so this alpha delta t should be less than equal to 1 divided by 100 delta c should be less than 1 divided by 100 alpha. So if your delta t is not less than that value then what will happen to if you use an explicit method is your temperature will grow with respect to time you will get but it will not reach to a steady state situation and second thing your number will approach towards infinity each picture which you will generate will be just a numerical artifact it will not represent the real world temperature distribution. So let us suppose if you are doing a problem where which is a steady state in 1 hour but if you decide a grid size let us suppose you have a time step restriction that you have to use a time step less than 1 second then you will need 3600 time marching to reach to steady state in an explicit method. Whereas if you use an implicit method there is no such time step restrictions you can use 1 minute as let us say as a time step and you can do 60 time marching. So here you need 3600 time marching here you need 60 time marching however per time marching this takes less time explicit method takes less computational time implicit method takes more computational time and in a computational procedure ultimately any method whenever we propose we have to say that what is the price involved with a method and where the which is called as computational cost like when you have to run a computer you have to pay electricity bill you have to pay not only for let us say your computers but also for running ACs if it is a big system there is a computational cost which is involved. So any simulation we are doing by any method we want it should be done in least time. So which is better that we that varies from problem to problem I will discuss advantage and disadvantage in the next slide but before I go into that I would like to discuss what we call as computational stencil I would like to do this with an animation before I do that I will show you in this slide in an explicit method all the neighbors are from the previous time level in an implicit method all the neighbors are from the new time level. Now what I am doing here is that these are the 5 values which are let us go to the animation I will show you through an animation so that you will get a better feel of this computational stencil what I would basically like to draw your attention is that in this figure there are 9 yellow circles. So we want to calculate the temperature at this 9 yellow circles. So what I will show is that the 5 values of the previous time steps are used to calculate the value of the next time level for this. So this stencil will move 9 times so that so we solve 9 algebraic equation in this case note that and in those 9 algebraic equation there is only one unknown. So I will show you an animation where you will get an idea. So let us take an example let us suppose there is a plate which we have taken from a furnace where let us suppose that this is your left hand side this is at 100 degree centigrade so all these blue circles are at 100 degree centigrade all these are at 200 degree centigrade. Let us assume that the boundary condition is such that these are at 300 degree centigrade. Let us assume that these are at 400 degree centigrade. So we are just saying that the left wall 100 bottom wall 200 right wall 300 top wall 400 degree centigrade these are coming from a boundary condition. Let us suppose this plate is taken from a furnace at 250 degree centigrade so all these yellow circle we will use an initial condition so that will be 250 degree centigrade. So this will form the initial picture of your movie. Now when you are using explicit method let us assume that you use an explicit stability criteria and let us assume that you got a time step of one second. So now the idea is how to obtain the value of temperature at one second how to obtain the picture of temperature distribution in one second. Now to do that I had shown you those points for the time at one second the boundary values as it is a Dirichlet boundary condition. So the blue circle values will be 100 200 300 400 on different walls. Now we have to obtain the temperature on 9 yellow circles. So here I am showing you a computational stencil for the first point okay now so for the temperature at this point of the new time level it uses temperature of the same point of the previous time level it is north neighbor, east neighbor, south neighbor and west neighbor. But note that south neighbor and west neighbor are blue circles where this has a temperature of 100 from Dirichlet boundary condition this has a temperature of 200. So in this case two values are coming from the boundary. However you know the value of the yellow circles it is coming from initial condition this is 250 this is 250 this is 250 three values 250 one boundary value 100 another boundary value 200. So you know all these five values you solve this algebraic equation using the five neighboring values and obtain the temperature of this point. But this you will get temperature of only one point you have to obtain for various point. So this way what I am showing you is that five values are picked from the previous time level and certain neighbors are used north south east west neighbor are used and we use an expression corresponding to this linear algebraic equation. So here you have nine equations note that in each equation you have only one unknown this arrow is pointing for the new time level at only one yellow circle. So this is what is called as computational stencils. So you get feel that how nine different values in this case are obtained through this computational stencils. So this is just an algebraic equations but with this you get a feel that how iterations proceeds. Now as far as the scrolling is concerned in programming language how do you scroll you scroll by using loops. So let us say in this case i is varying from 1 to 5. So you have to scroll how many yellow circles you have three in the x direction three in the y direction. So what loop you will use i is equal to 2 to 4 j is equal to 2 to 4. And for neighbors when it is a east neighbor you use i plus 1 for west neighbor i minus 1 north j plus 1 south j minus 1. So that way you use running indices. So there is a loop which you use you use an expression which is in terms of the running indices of the loop and solve it. In an implicit method this is a computational stencil. Now here the expression is search that there are five values from the new time step which are involved. Now here this yellow circle is of the new time level which I do not know. So in this case I get an algebraic equations where there are three unknown. This boundary value does not change with respect to time. So in the previous time also it was hundred next time also it is hundred. This is also two hundred but three yellow circle values in previous time step they were two fifty but in the new time step they are unknown. So in this equation you have note that three unknowns. Why? Maximum except there will be five unknowns. Why? Three is coming because there are two boundary values from the boundary conditions you know that. Here there are four unknown values, three unknown values, four unknown values, five unknown values and so on. I am just seeing the number of yellow circle as this is moving and that are the number of unknowns. I will repeat it again. So just see as this stencil moves how many yellow circles are there in that red finger at new time level. Three yellow circles are there. So three unknown values in this equation. Four unknown values in this equation. Three unknown values in this equation. Four and so on. So here you get nine algebraic equation. In some equation there is three unknown. In some equation there are four unknown. In one equation here in this case there are five unknowns. So here you get again I would like to point out that here you have a system of algebraic equation where each equation is coupled with another algebraic equation because there are more than one unknowns. And this we solve in COD in an iterative manner. Let us come back to the lecture. With this I had shown you that I had shown you you can say a graphical or a pictorial form of the equations, how the equations are used in calculating the values at the new time levels. As I mentioned that explicit approach is relatively simple to set up and program. On the other hand it has a certain disadvantage that it results in long computational time to make calculations over a given interval of time because there is a time step restriction. So you may have to do more time marching. Whereas an implicit method you can use a larger time step. So you can do a lesser number of time marching. However it has a disadvantage that the computational time per time step is much larger. For an inherently unsteady problem this method is more accurate because in those cases you need a temporarily accurate solutions for which you need a very small time step. So after showing you the talking about the explicit implicit method discussing this tensile. Now I am showing you the solution algorithm for explicit method. We enter the inputs like thermo physical property. Note that you can have a various type of induction heat problem depending upon the material of the plate which you are using. You can have different problems depending upon the size of the plate and as a user you have to enter the number of control volumes in the x and y direction, the boundary conditions and this is what is called as convergence criteria. Because you want to stop this if let us suppose if you take a video camera and you want to shoot a movie. But typically let us say in a fluid flow and heat transfer situation what happens is that initially the picture keeps on changing but finally a state reaches when then the picture does not change with respect to time. This is what is called as a steady state condition. So to check that the data in this case the temperature which we are getting at different time steps whether if the it is approaching towards a steady state. To check for steady state the temperature difference should die down between the two consecutive time step which should reach 0 but the computer does not understand what is 0. So in any numerical method we have to use what we call as convergence criteria. So as a user you have to define convergence criteria we define typically by symbol epsilon. Second step is a grid generation where we calculate all the geometrical parameters of the all the control volume. Third step there is a mistake here this should not be initial guess it should be initial condition initial guess we use in a steady state formulation. So here we have to give initial condition for temperature as in the stencil I was saying it is 250 degree centigrade and Dirichlet boundary condition for temperature. Note that here I am showing you Dirichlet boundary condition and a non Dirichlet boundary condition has two different steps. There is a reason behind it later on I will show you that when you use a Dirichlet boundary condition temperature value is a constant at the boundary which does not change with respect to time. But when you use a non Dirichlet boundary condition like constant heat flux or a convective boundary condition when you discretize the temperature value is a function of border value the border value is an interior value which changes with respect to time. So the boundary point value changes with respect to time. So this step 4 non Dirichlet boundary condition is inside the loop of time marching. The step 3 does not change with respect to time. But in this step 4 the boundary temperature changes with respect to time and it is put in a loop for when we once we reach the steady state condition. Now once you solve for want to solve for you start with an initial condition that is once you know the temperature when you want to calculate temperature for the new time level at each time step you have two matrices. One is the temperature of the new time level and second the temperature of the old time level. Now what you do is that otherwise what may happen is that maybe you will say that if a problem reaches to a steady state in let us say 1 hour and in explicit method if your time step is 1 second then to create an animation let us suppose you want to store the data of each of this time step result then you will need 3600 matrices rather than doing that what typically is done is that at each time step we have two matrices one from the previous time level second from the new time level and we compare for convergence for steady state condition. If the temperature has not reached to a steady state condition whatever temperature we get we take it as a old value. So before we move into time marching whatever temperature we have we take it as a old volume. Now here the solution algorithm or the implementation details which will be shown later is such that we program what we have understood from the conservation law. What we have learned from the conservation law is first at the phases of the control volume we should calculate the heat fluxes and once you calculate the heat fluxes then you balance the heat fluxes in the control volume and calculate the temperature. So the temperature calculation here we are not solving the direct linear algebraic equations but here it is shown in a more intuitive or more connected to the conservation law where the idea is first we will calculate the heat fluxes at the phases of the control volume using the temperature of the old time level and in the second step we will do a balance and calculate the total heat. So with this step we get a total heat gain by conduction which is used to calculate the temperature. This will be clear in the following slides where I discussed the implementation actually in mode. Then we check for convergence what is convergence this is just to check whether the temperature has reached to a steady state. One of the simplest way to check is that we have got two pictures. First picture let us say temperature of the previous time level, second picture of the temperature at the next time level. Now to check for the steady state condition what you do is that the simplest way of is to take the difference of the temperature of the two time level. So, you get another matrix which is let us say t minus t volt and in that matrix pick the maximum value that is the maximum difference value. And check that maximum difference whether it is reached to let us say 10 to power minus 3 or 10 to power minus 4 which is which you are defining as practically 0 which is your convergence criteria epsilon. But rather than using that condition typically in CFD calculation we calculate what we call as the root mean square value of the temperature difference whose expression is given by this. So, we calculate a root mean square value of temperature and that value we check for converges. If it has converged then we have reached to steady state otherwise we go back to step 4 and continue this various steps. So, this loop goes on this time marching loop goes on till we reach to steady state condition. Now I was saying that first we calculate the heat fluxes then we calculate the temperature. So, to show this I will show you some one animation where it will be clear to you. Now let us take what I am showing you here is that there is this is a one dimensional problem. So, what is a one dimensional problem typically it comes into picture when you have a thin sheet heat transfer across a thin sheet. So, when you take a sheet the dimensions in the direction of thickness is very small as compared to other two directions. In this case we take the heat transfer as one dimension. So, let us say that this length corresponds to the thickness of the sheet. Now what we are doing in this case we are defining the sheet let us say the thickness of the sheet into 5 equal parts this is the first part, second part, third part, fourth part and fifth part. Let us go to the animation which will be coming in this window where I will show this. So, let us suppose this is the thickness of the sheet and you are dividing into 5 equal parts and then you get 5 control values in the x direction and these are the boundary gate points. Let us suppose as the first let us say 1D heat conduction problem typically you might have seen is left wall is at 0 degree centigrade right wall is at 100 degree centigrade. So, that let us suppose here temperature is 0 degree centigrade, here temperature is 100 degree centigrade from boundary conditions and let us assume that here we have to take an initial condition. So, let us assume that this plate initially is at ambient temperature. So, let us assume that all this yellow circle have a temperature of 30 degree centigrade. Now what I am saying is that we will calculate the heat fluxes on the phases of the control volume. Now tell me how many phases are there of the control volume? How many control volumes you have? 5 control volume. A phase is common to 2 control volume which is the first control volume. First control volume is from here to here. So, this is your first phase center, this is your second phase center. This phase center is common to this. So, this will be the third phase center. This is the fourth phase center. This is the fifth phase center. This is the sixth phase center. So, note that in this case there are 7 grid points. There are 5 yellow circle, 5 control volumes, interior control volumes and there are 6 phase center. What we do is that at this 6 phase center we calculate the heat gained by conduction and the stencil for that is given by. So, using the temperature of the old time levels where Tb1 is 0 degree centigrade, Tb2 is 100 degree centigrade. We will calculate the heat flux in the x direction. What is the expression for heat flux in a discrete form? What is the discrete representation of Fourier Laff heat conduction? To calculate the dT by dx here, here we use the one sided forward differencing. Temperature here minus temperature here divided by the distance between the two. What is temperature here? 250 degree centigrade. What is the temperature here? 0 degree centigrade. What is this distance? Let us take the total length as 1. So, the coordinate here is 0.1, 0.3, 0.5, 0.7, 0.9 and 1. So, the distance here is 0.1. So, 250 minus 100 divided by 0.1 multiplied by minus k will give you heat flux at this phase center. If you want to calculate heat flux at this phase center, it will be temperature here minus temperature here divided by the distance between the two. What is temperature here? 250. Here also it is 250. This will be, heat flux will be 0. So, this way we calculate the heat flux. So, this is the computational stencil to calculate the heat flux. What I am showing you in this animation as green square is the, you can call it grid point for heat flux in the x direction. Where is it located? It is located at the phase center of each control volume. How many green squares are there in this figure? Six. There are six phase center, six green square. So, this is going as per law of conservation of energy. That if you know a temperature distribution, you can calculate heat fluxes at the phases of the control volume. We can use this heat fluxes and obtain temperature distribution like in a control volume. Unstudy term is rho C P T N plus 1 minus T N divided by delta T is equals to Q W minus Q E. On the west phase, it is entering. On the east phase, it is leaving. So, we calculate temperature at the old time level using the heat fluxes of the previous time level. So, two heat fluxes which corresponds to the east and west phases are used to calculate the temperature at the cell center. So, this is a computational stencil. So, with this stencil or this animation, you get a feel how you can write a program. So, here what you will see which you do not find in most of the CFD books is that in the CFD books, you get expressions. Here, I am showing you figures, animations to give you a feel so that you get you could understand it in a better way those expressions and you can start programming. Later on, I will show you implementation detail. I will even give you pseudo codes which is basically a program but not in a programming language with proper looping which you can directly take to a programming language and use proper syntax and convert it into a program. Let us come back to the slide. So, this is the computational stencil. I had earlier shown you the computational stencil in an explicit method and implicit method but here I am showing you computational stencil for the implementation detail to calculate the heat fluxes and finally the temperature in case of 1D and steady state conduction problem. So, this is the expression which I was talking of. This is the first step which I was saying where this expression is in terms of heat fluxes and this heat fluxes are expressed in terms of the Fourier law discrete form of the Fourier law of heat conduction. Just to give you an idea, let us take let us do a simple calculation. Let us assume that the surface areas are unity. Let us assume that the density thermo physical property density specific conductivity is unity. Let us take delta x and delta t as 1. So, with this values if you substitute into this expression. So, this is 1, k is 1, surface area is 1, volume is 1, rho is 1, Cp is 1. So, with this in the sample calculation you will get an expression as Tpn plus 1 is equals to Tpn plus qw minus qe. So, let us take a case where left wall is at 0 degree centigrade, right wall is at 100 degree centigrade and let us suppose initially the plate is at 50 degree centigrade. So, the all yellow circle have a temperature 50 degree centigrade. So, what happens in the first time step is shown here. So, if you want to calculate dt by dx here, this will be temperature here 50 minus 0 divided by this distance. So, what is this distance? This is 0, this is how much is the delta x in this case. This is divided into this is 0, here I am taking it as a 0.05, this is 0, 50 minus 0. Here in the previous slide let me show you that delta x I had taken it to be 1, what will be delta x by 2? 0.5, okay. So, what is delta x? This much is delta x, this distance is delta x, this distance is delta x. You can see that at the end, the distance between the two circle is half. So, this is delta x by 2, delta x is 1. So, it will be 0.5. So, when you calculate d t by d x, it will be 50 minus 0 divided by 0.5. So, d t by d x comes out to be 50 divided by 0.5, which is 100. But, you have a minus key also, k is 1. So, you get minus 100 here. Similarly, when you calculate d t by d x here, it is d t by d x will be 0, because 50 minus 50 will be coming here. Here also d t by d x is 0. Here also d t by d x is 0. On the right hand side, on this phase, d t by d x is 100 minus 50 divided by 0.5. So, d t by d x will come out to be plus 100. But, when you multiply by minus k, it will come out to be minus 100. So, with this, you have calculated q x on all the phase center. Note that there are 1, 2, 3, 4, 5, 6, 7 values here. But, q x, you have 6 values, because there are 6 phase centers. Now, you do a balance. So, this is the heat flux on this phase. So, what is this phase for this cell center? This is the west phase. This is q w for this cell control volume. q e is 0. So, total heat, which is going in is minus 100 plus 0, which is minus 100. So, this is minus 100. This is minus 0. This is minus 100. And then, you have to add with t p n. What is the value of temperature here? 50. So, 50 minus 100 makes it minus 50. So, this way, you can calculate the temperature. So, in this slide, I am showing you a sample calculation. So, whatever has been discussed in the previous slide, using those equations, you can obtain this number distribution. Typically, in an undergraduate course or in a CFD course, to get the feel of how the number grows, this is the sample calculation. These are also used as a sample in assignment problems or in the examinations. So, this is a solution algorithm for one-dimensional situations. Due to shortage of time, I will not be able to go into this. I had already discussed in detail this for a two-dimensional unsteady state conduction. Now, I will go to the steady state formulation. In a steady state formulation, there are in computational fluid dynamics, there are two types of formulation. One which is called as unsteady state formulation, where you use an unsteady state equation. But many a times, in fact, most of the times, we have problems which reaches to a steady state. So, you feel that why to start with an unsteady equation? Why cannot we start with a steady state equation? This is what is called as a steady state formulation. So, in case of, let us say, unsteady state conduction, you have unsteady term also. But here, in this case, in a steady state case, you do not have del capital T by del small t. So, in this case, you just have to solve del square T by del x square is equals to 0. In fact, you know the analytical solution of this also in undergraduate, from undergraduate courses. And this is the boundary conditions for those. And this is the finite volume discretization of this govern equation. This had already shown in my previous slides. And this can be expressed into this form, where temperature at a particular point is a function of its west neighbor and its east neighbor. And these are the distance between the two consecutive cell center. This equation, when you apply for, let us say, note that this equation, discretized equation is applied only on this 5 yellow circles. Note that this equation is not applied on the boundary points. Note that we apply the conservation law only on this 5 yellow circles. At the boundary, we apply the boundary conditions. So, for this 5 yellow circle T 2 to T 6, we substitute this delta x delta x E delta x W. And we end up with this coefficient matrix. So, in a steady state conduction, I was mentioning that you start with a steady state equation. You do a finite volume discretization, where you do not have any time term. And here you end up with certain set of algebraic equations. Like in this case, there are 5 yellow circles. You get, this is the first algebraic equations, where the unknown is T 2. This is the second algebraic, sorry, in this first algebraic equation, the coefficient of T 2 and T 3 are non-zero. Note that the coefficient of T 4, T 5 and T 6 is 0. When you go to T 3, T 2 is what? T 2 is a border control volume. T 3 is an interior control volume. T 2 has one neighbor, which is boundary. So, for T 2, this boundary value comes from the boundary condition. So, it has only 2 non-zero coefficients. When you reach to the interior grid points like T 3, T 4, T 5, there are 3 unknown, 3 non-zero coefficients. When you reach to the border grid point T 6, then there is only non-zero coefficients are for 2 values. So, you see that the structure of the coefficient matrix is, this is the main diagonal, this is the upper diagonal, this is the lower diagonal. So, you get what are the nature of the coefficient matrix is what we call as a tri diagonal matrix. So, this is, there are lesser number of zeros as compared to non-zero values. When you only have, let us say 5 yellow circles, but instead just imagine, if you have millions such yellow circles or control volumes or grid points, then there will be, the percentage of zeros will be much more as compared to what is here. So, when there are more zeros in a coefficient matrix, we do not want to take the inverse of this matrix, because there will be lot of unnecessary computations. So, instead of taking the inverse and multiplying with this vector, which is called as a non-iterative, which is called as a non-iterative or a direct method, we solve by a non-iterative technique, what is called as a tri diagonal matrix algorithm or by an iterative technique, which is called as a Jacobi as Gauss-Siedel method. When you go to multidimensional case, even if you are using a tri diagonal matrix algorithm, still then it becomes iterative in nature. So, what I am trying to emphasize in a steady state formulation, you march iteration by iteration. What is mean by iteration by iteration? I will show in the next slide. Let us take a case where left hand side left wall is at 0 degree centigrade, right wall is at 100 degree centigrade, T b 1 is 0, T b 2 is 100. Now, from this I can write an equation that for i is equals to 2, which is the left border cell, T 2 is equals to 2, T b 1 plus T 3 divided by 3. From this, this is the first equation, which I am showing you here. The last equation, which is coming from here is T 6 is equals to 2 b 2 plus T 5 divided by T 3. For the interior grid points, you get an equation like this. So, instead of initial condition, in a steady state formulation, you start with an initial guess, which is 50 in this case. And if you are using a Jacobi iteration, then after first iteration you get these values. If you are using a Gauss-Siedel method, you get these values. You can check these values later on. So, in a Jacobi iteration, the idea is that the previous time level value of T 3 is used to calculate T 2. And in a Gauss-Siedel method, the updated value of T 3 is used to calculate T 2. So, this is the solution algorithm for a steady state formulation. Professor Puranik had taught finite difference method and he has taken a problem in a steady state formulation. And this is the solution algorithm corresponding to this. The only difference here is that, here in this case we march iteration by iteration instead of time by time. The steady state formulation for 2D heat conduction is shown here. This is the finite volume discretization of this steady state heat conduction equation in two dimension. If you take a problem like this, where left wall is at, let us suppose 100 degree centigrade, bottom wall is at 200, right wall is at 300, top wall is at 400 degree centigrade. Then you get this type of algebraic equations. Now, here you see that instead of tri-diagonal matrix, you get a pentadagonal matrix. Because this is the main diagonal. This is the upper diagonal. This is the lower diagonal. And you have one more diagonal where you have a non-zero value. Here on the upside and on the lower side. So, there are 5 diagonals and you have a non-zero coefficients. And this is a constant matrix. Note that in the constant matrix, you have a boundary conditions which are coming. All these values are coming from the boundary conditions. So, this algebraic equation is solved in a steady state formulation. Now, the solution algorithm is similar to the steady state formulation. Only thing here is that we move iteration by iteration. Now, I will move to the implementation details for an explicit method. I will show this implementation detail through an animation. So, let us go to the animation which will be coming in this window. So, this is a plate. We want to solve 2D conduction problem in this. We are drawing equi-spaced vertical lines, equi-spaced horizontal lines. And we are dividing into 25 control volumes. These are the 25 interior grid points. These are the boundary grid points on the left, bottom, right and top walls. These are the grid points for temperature. Just to give you a feel, a pictorial representation before you go on to the code development. I am just showing you that these are the points for temperature. So, how many yellow circles are shown in this figure? Seven in the x direction, seven in the y direction. Now, I will show you the grid point to calculate heat fluxes in the x direction. Now, a phase is common to two control volume. Let us suppose if you take a vertical phase, it is east phase of this control volume and it will be a west phase of neighboring control volume. Such as this phase center is east phase for this control volume. It is a west phase for this control volume. So, when you want to calculate heat flux at this phase, you want to calculate only ones. You want only one running indices. In a matrix, you need one number corresponding to this phase center. So, to avoid duplication, you need to use a convention that what should be the running indices for the calculation of heat fluxes in the x direction. Now, here the convention, right now if you look into the vertical phase center, there is a vertical phase center here also. There are two vertical phase centers. This is the west phase center for this control volume. This is the east phase center. Now, the convention which is being used here is the qx for a vertical phase center. The i, j will corresponds to the cell on which this square is lying on the east phase center. This is lying on us at the east phase center of this cell whose running indices is j, i. So, the running indices of qx will be j, i. Similarly, I can have a vertical phase center here also, but the j, i will corresponds to the j, i of this cell center. What is the running indices for this cell center? It is j, i minus 1. So, here you will get qx j, i minus 1. Similarly, so on the positive phases, we are defining the running indices corresponding to the positive phases. So, the qi at j, i, this is the north phase center of this temperature which has j, i. So, this will be your qi j, i. Here also you can have qi, but the qi here will be j minus 1, i. When you program, you have to create loop to calculate qx and you need to create loops in a proper way. So, in this control volumes, you have 5 control volumes in x direction, 5 control volumes in the y direction. Now, what I will show you is the vertical phase center or the grid points for the qx. So, how many green squares are there in the x direction? 6 green squares in the x direction, 5 in the y direction. What is the loop? What will be the loop which you use for scrolling through these green squares? The loop will be this is i is equals to 1, this qx is lying on the east phase of this control volume. What is i and j corresponding to this? i is 2. So, for this control, this square i is 2, i is 3 for this, i is 4 for this, 5 for this, 6 for this. So, to scroll through when you want to calculate qx, your loop will be 2 to 6 or 2 to i max minus 1 because i max is your 7. When you want to scroll in the y direction to calculate qx because you have to scroll through all this, 6 in the x direction, 5 in the y direction. So, when you move into this, this will be j is equals to 2, j is equals to 3, 4, 5, 6. So, the loop in the x direction is 1 to i max minus 1 and in the y direction it is 2 to j max minus 1. Just to give you an idea, the number of yellow circle is 7 by 7, the number of green square is 6 by 5, note that. And when you, these are the grid points for horizontal phase center, how many are they in the x direction, 5 in the x direction, 6 in the y direction. So, the loop to scroll through the qi will be 2 to i max minus 1 and 1 to j max minus 1. So, with this figure I am trying to give you a feel that if you want to develop a program, you have to, here we are proposing a, using a method where we first calculate qx on vertical phase center, we calculate qi on horizontal phase center. Now, to calculate this, you need to use a loop and for any, so to calculate qx at any point, let us say qx I want to calculate here. How do we calculate qx here? I use temperature of its east and west neighbor. The temperature here is i plus 1, here it is i, ti plus 1 minus ti divided by the distance between the two points. Let us go back to the slide. So, this is a, what I called as a pseudo code. This is not corresponding to a programming language, but this shows you how to calculate heat flux in the x direction, what is the proper loop to calculate heat flux. This line shows you how to calculate the distance between the cell center. This line shows you how to calculate the heat flux in the x direction. So, this part of the program is to calculate heat flux on all the green squares. This is the part of the program to calculate heat flux in the y direction at all red inverted triangles. Note that here it is 2 to j max minus 1, that is 2 to 6, 1 to 6. Here it is j is 1 to 6 and i is 2 to 5. You can correlate this with the figure which is given here. And this is the loop to scroll through all the yellow circles which are inside 2 to 6 because now you want to do a heat balance in all this 25 control volume. So, you scroll through 2 to i max minus 1, 2 to j max minus 1. This is the line which you have to program to calculate total heat gain by conduction. This is the heat flux from the best phase. This is the heat flux from the east phase. This is the surface area of the vertical phase delta y. So, with this much heat transfer in the x direction, this much heat transfer in the y direction, total heat gain by conduction. This is the total heat gain by volumetric heat generation. You add this to multiply by delta t divided by rho Cp delta v. So, this is the change of temperature with respect to time. You add it to previous time level value and obtain the temperature for the next time level. So, this is what I call as a pseudo code. You just have to convert this into a programming language to develop a program. However, we have given you a silap code where this has been implemented exactly. This is a pseudo code for one-dimensional steady state conduction. I will not discuss this. So, let me go to the next slide where I will discuss the example problems. So, this is all CFD development part. Now, there is a second part which is CFD application where I will take examples. I have taken two type of problems, one-dimensional problem and two-dimensional problems. Let us take, I will show you a movie for temperature distribution. Let us suppose left wall is at 0 degree centigrade, right wall is at 100 degree centigrade. So, left wall is at 0 degree centigrade, right wall is at 100 degree centigrade and initial condition is 30 degree centigrade. Let us say ambient temperature. So, all the cello circle are at 30 degree centigrade, left wall 0, right wall 100. So, you know that we have an excess solution for this and you know that the temperature variation is linear. So, this is a good problem to start with. If you want to develop a program, it is always good to start developing program for a problem which already has an analytical solution. So, you get a confidence that the results which you obtained from simulation is correct. So, in this case let me show you a movie corresponding to this. Now, the movie which I will show you is x is shown in the x axis, temperature is shown in the y axis. Let us come to this window and I will show you a movie corresponding to this case. So, here you can see that actually in this case I have taken the number of grid points as whatever is given here. So, corresponding to this yellow circles you can see on the x axis I have x coordinate and this on the y axis I have a temperature. So, you can see on the y axis how it is varying 0 to 100 degree centigrade. Note that this is a conduction problem without volumetric heat generation. So, if left wall is 0, right wall is 100. So, the intermediate temperatures will be limited between 0 and 100. So, the left wall is 0, right wall is 100 and the initial condition is 30 degree centigrade. So, this is your initial picture. So, this is the temperature as a function of x plotted at time is equals to 0. At different time step this picture will change and you can get a movie which corresponds to this. Finally, it reaches to a linear variations which you can let us come back to the slide. So, this is your linear variations. So, the movie continues and finally, it reaches to this steady state condition. Now, let me take a problem where left wall is at 100 degree centigrade and on the right hand side I have a non-dischlep boundary condition. Let us suppose right side is exposed to an ambient well let us say the h value is by force convection convective heat transfer coefficient is 100 watt per meter square Kelvin and the ambient temperature is 30 degree centigrade. So, for this case I will show you a movie. Let us come back to the small window. So, this is a movie corresponding to that. I would like to draw your attention that left wall is always at 100, but the right boundary temperature note that this is changing with respect to time. So, whenever you have a non-dischlep boundary condition the boundary value is a function of border value. Border value is a part of the inside the domain and which changes with time marching. So, the boundary point value changes with respect to time and whenever you have a non-dischlep boundary condition that non-dischlep boundary condition steps for non-dischlep boundary condition has to be inside the time marching loop. So, the steady state let us come back to the slide. So, the steady state temperature profile will still be linear in this case also. However, note that in this case the boundary value the right boundary where you have a non-dischlep boundary condition value varies with respect to time. You can vary the value of h and you can get different variation and this is an analytical solution corresponding to that. Now, let me take a problem with heat generation. So, I will show you an animation where let us suppose left wall is 0, right wall is 100, initial condition is 30. We are taking the same case. Now, in the earlier case heat generation was 0. Here volumetric heat generation is 30 kilowatt per meter cube. Let us go to the small window to see the animation corresponding to this case. It is a 1D conduction problem with a dischlep boundary condition at both the ends. However, it has certain volumetric heat generation. Now, it has a volumetric heat generation. The point here to note is that let us go back to the animation. This is an animation for left wall at 0 degree centigrade, right wall at 100 degree centigrade and there is a volumetric heat generation. So, the animation corresponding to this, note that whenever you have a volumetric heat generation depending upon the magnitude of volumetric heat generation, if it is large enough then you may get a maximum of temperature which is greater than 100 degree centigrade. That is what is being shown here. Now, let us come back to this slide and I can show you the study state temperature distribution. What I want you to note here is that the maximum value of temperature is how is greater than 100. So, it is not limited between 0 and 100 when you have a volumetric heat generation. Now, let us take a case where this is the same case, but with the different boundary conditions and you get a study state temperature distribution as this. You have been given programs using the methods which has been discussed in this lecture which will be helpful for you in better understanding. Now, let us go to the second problem which is a two-dimensional problem. Left wall is let us say 100 degree centigrade, bottom wall 200, right wall 300, top wall 400 degree centigrade and for this let us I will show you an animation here. This is an animation for temperature control. The line here represents where the temperature is constant. So, there are different lines. If the line shows 100 means on that line the temperature is 100. So, the different values corresponding to different lines are shown here. There is a color also and there is a color bar which gives you the magnitude of the temperature. Note that here the temperature is bounded between the minimum value and the maximum value at the boundary which is 100 and 400. Let us come back to the slide. Here what I want to point out is that whenever you are solving a problem in CFD or let us say in research you do not know what results it will come, but you should be clever enough to know what result will not come or should not come. Like in this case, if left wall is at 100 degree centigrade and you see that the isotherms are let us suppose here temperature is 150. So, the temperature here is 140, 130, 120. The temperature in this nearby region is greater than 100. So, the heat transfer should be outward. On the bottom wall is at 200. So, on this line temperature is 200. Just above this on the left hand side the temperature is less than 200. On the right hand side the temperature is greater than 200. So, in this case there is a change in the direction of heat fluxes. So, what I want to emphasize here is that I say many a times that you create pictures and movies, but this pictures of let us say in this case the picture of a temperature distribution gives you an understanding of the direction of heat transfer which is shown through this figure. So, like for flow across a car if you are a designer, if you create a movie for flow across a car fluid dynamic movie and calculate the drag force, a particular design is giving you less drag then it has it has a particular relationship with the type of flow structures which are generated. So, what I want to emphasize here is that this pictures or movies has a relationship with engineering parameter. So, why a particular design is better that depends upon the nature of the flow and from the physics of the flow which you understand from this pictures or movies. Using the number distribution near to the bottom plate I had shown here how you can calculate the heat fluxes and total rate of heat transfer. I will take the last problem which is left wall 100, bottom wall insulated, right wall constant heat flux and top wall convective boundary condition. This is the discretization of the boundary condition for the different boundaries and I will show you an animation. Let us go to the small window. So, here you have a steel plate the initial condition is 30 degree centigrade, boundary condition is left wall 100, sorry a left wall 100, bottom wall insulated, right wall constant heat flux, top wall convective heat transfer boundary condition. Note that the line which are hitting the bottom wall they are always perpendicular at the wall. Also note that the line which are hitting the right wall the slope of the line is constant. Why? Because you have a insulated boundary condition on the bottom wall. So, at the wall deltie by del y should be 0. So, this line in each picture should hit bottom wall where deltie by del y which is a slope of the temperature should be 0. Similarly, on the right wall you have a constant heat flux boundary condition. So, deltie by del x is minus q w by k. So, the slope is minus q w by k which is a constant quantity. So, if you have developed a code and in each picture this thing is not happening then you should be clever enough to know that this results are not correct. That is what I mean by that in research you should know what results should not come as we do not have control over what results will come. So, by looking into this picture you should be able to understand whether it is correct or not. Like for flow across a car if you are having a domain where near to the outlet if you have a let us come back to the slide. I was saying that if for flow across a car near to the exit if you are having a recirculation region it means your domain size is not correct. Here I am showing you how the heat flux varies on the different wall. So, I will end here.