 Alright, so let us look at the first sort of a numerical exercise that we will go through. The last problem that we looked at when we just discuss these exact solutions that diffusion equation with a source is what we will take for the discussion of this finite difference approach. So before that, let me just quickly outline the overall sort of numerical solution methodology that we follow and this is true for whether you are doing a finite difference or whether you are doing finite volume. So what we normally do is we will start with governing equations and then we carry out a process called discretization on a measure of grid which will result in a set of algebraic equations which then are solved using some solution algorithm after which we obtain the numerical solution. So the numerical solution then is essentially obtained on the discrete points which are the grid points or the mesh points as you say. After you obtain that numerical solution, we end up doing post processing of the solution data and then that involves essentially the analysis of the problem or the solution. So for example the standard post processing you can think of one example I can give here is let us say you solve the fluid flow equations and the energy equation. Thereby you obtain at the end of it a temperature field. So from the temperature field you can compute the temperature gradient for example at the wall from which you can calculate the data for a heat transfer coefficient. So that is how you can think about as an example of what I mean by a post process. So this sort of a procedure is very general it does not have to be restricted to only finite difference essentially the same thing that you will follow for the finite volume which will be taken up later. So let us quickly look at what we are talking about here. As I said this is from the point of view of introducing to you one problem where the solution is going to be worked out using finite differences. In the sense that the procedure that I am going to talk about is actually restricted to this particular problem. So I am not going to talk about general finite difference philosophy much and so on. With respect to this specific problem that we just looked at we will solve it. So what was that? It was basically this transient flow between infinite parallel plates and the governing equation was this diffusion equation with the source term with that source equal to minus 1 over rho dp dx and which was a constant. And there were these boundary conditions the lower plate is stationary upper plate is moving at constant velocity of capital U and the initial condition is everywhere at time equal to 0 the velocity is 0. So the idea is that this is what we will call a time evolution problem or a time marching problem in the sense that starting from the initial condition we want to find out how the solution evolves as a function of time. So in any of these numerical simulation problems what you do is you first come up with a grid which is basically division of your domain of interest into a number of these mesh points. So in this particular case again as you can see we have a y domain only going from 0 to capital H. So that is what I have shown here through this sketch. On the left I have y equal to 0 corresponding to which I am placing one node called let us say the left boundary node if you want and then I am placing these successive nodes at equal distance in this case I am calling that as delta y till I end the domain with final node called capital N as the right boundary node at y equal to H let us say. So these are the points 1, 2, 3 etc all the way up to N is where the numerical solution will be obtained. We will have no information whatsoever about the U data because we are solving for U as you can see as a function of both y and t. So we will have no information in this region between the grid point. So one thing that I would like to point out even though perhaps that is not how it has been done in the Sylab code is I will prefer if you actually plot your results as a set of discrete data point because in principle what we are doing as a part of CFD simulation is that we are actually obtaining the data at these discrete points 1, 2, 3, 4 up to N. So it does not make really sense to come up with a continuous distribution it is better if you plot this as a separate discrete data point. I am not too sure whether it has been done that way or not in the Sylab the RA is saying no which is fine but for your own sake I think it may be useful to do it later in that sense. So then 1 and capital N are the boundary nodes and delta y in this case is the grid spacing. We are going with a uniform grid spacing it does not have to be uniform but for one representative case of this finite difference I am just going with a standard uniform distribution. So the finite difference approach involves basically utilizing the Taylor series expansion to replace the partial derivatives with respect to time and with respect to space using corresponding algebraic expressions let us say. So as I said I am going to restrict myself only to this particular problem. So let us see what is done in that case then. So I just consider one general interior node point I and the neighbors for these on the left are I-1 and I-1 as most of you would know. So what I do is I express the variable U at let us say I plus first node in terms of its value at ith node and the usual Taylor series standard expansion. So here I am including du dy at the ith node times delta y plus d2u dy squared at the ith node times the second term and so on. You do not have to restrict yourself only to second term but I am going to do it in this case. Similarly I am going to express the value of U at the I-1 node in terms of the value of U at I and all I am doing then is now I am moving by a distance of minus delta y in that sense to the left and accordingly the Taylor series is written in that fashion. So this is let us say I am calling 1 and 2. If you add 1 and 2 what you will see is the term involving this du dy will simply go away. This guy will get doubled up and if you rearrange what you obtain is an approximation for the second derivative of y sorry second derivative of U with respect to y at the node point I which is given as U at the I plus first node minus 2 times U at the ith node plus U by U of I-1 that is the U value at the previous node the whole thing divided by delta y squared and many of you would know this obviously and many of you would say that this is a second order central difference and so on. I am not even bothered about getting into those details. Primarily the reason is because we are not going to follow the finite difference methodology in the next three days. We are going to focus on the finite volume. So it is just to introduce how to obtain approximations for derivative expressions using Taylor series expansion is what I am doing here. So this is one way or one approximation for the second order derivative with respect to space of the variable at the node point I. Now that is what I want because what I have is the second order derivative with respect to y of U. I also have a first order derivative with respect to time of U. So let us see if we can employ the same thing. So I have two options here either I express the value of U at the node I at a later time instant which I am calling by superscript t plus delta t using the value of U at the same node I at the previous time step plus a Taylor expansion in time now. So I will have du dt at the ith node times delta t not even bothered with adding more term. If you simply rearrange this you will obtain an approximation for du dt at the node I as the value of U at I later time minus value of U at I previous time divided by or current time I should say previous by divided by delta t. And normally for the sake of representations this t and t plus delta t perhaps can be shown using a superscript n and n plus 1 etc. So n is the current time level n plus 1 is the next time level which will be essentially t plus delta t this is one option or the other option is going backwards in time we went forward in time in the case on the left hand side we can likewise go backward in time. So then in this particular case we will write that U I at t minus delta t will be then expressed in terms of U I at t minus then du dt times delta t. And if you rearrange this you will see that this comes out as U I at t minus U I at t minus delta t over delta t and again using the same superscripts that t level is n t minus delta t level is n minus 1 and delta t. So then I have two possibilities in this case as far as the time derivative is concerned and I have decided that I will use this one particular way of representing the second order differential with respect to Y. So then let me complete the nomenclature then what we have come up with here for the second order difference is what is popularly called as a central difference as most of you would know. The one on the left hand side here for the time derivative is what is popularly called as a forward difference and likewise the one on the right here is what we will be called as a backward time. So then I can choose whether I want to go with a central difference here and a forward difference for the time or a central difference here and a backward difference in time. Accordingly the finite difference equation corresponding to the original differential equation which we started with for each interior node I can be written using either a forward time and a central space approximation or backward time and central space approximation. So as you can see the space derivative is evaluated at the current time level always here whereas the time derivative can be either a forward derivative or a backward derivative and the source as it is which is a constant that remains the way. So let us see the first option here which is what will be called as the forward time central space form or scheme which is also called as an explicit method and we will see why. So if you simply rearrange this expression you can write it as the u value at the ith node for the later time instant is equal to the u value at the ith node at the current time instant plus all this business which is also at the current time instant. This is written for each interior node I let us say 1 by 1. So what you realize is that since I am beginning with an initial condition at time equal to 0 I am going to assign the u value at time equal to 0 everywhere. So the solution in that sense at time equal to 0 is known and then I want to go from this solution equal to whatever it is at time equal to 0 to the solution at time equal to 0 plus delta t let us say which is at delta t. Utilizing the solution at the zeroth level I can explicitly write the solution at the delta t level without really doing anything other than just visiting each grid point in a sequence utilizing this right hand side which is all known. So plugging the values of the known variables on the right hand side for the previous time step I can find out 1 by 1 the value of the solution at the later time and that is why this is called as an explicit method that is we know the solution explicitly in terms of the known solution on the previous side all right. So this is all done for the interior nodes remember we have the boundary conditions also. So in this particular case in this particular problem as we have noted the boundary conditions are particularly simple all that is is that the u value at y equal to 0 is 0 and u value at y equal to h is equal to capital U. So y equal to 0 corresponds to the i value of 1 and similarly y value of h corresponds to i of capital N. So I simply write that u of 1 this really is immaterial in that sense because at all times it is the same however just for completeness I am writing it at the new time level equal to 0 similarly u at the right boundary point is equal to capital U for all. The N plus 1 time level marking is required if you have time dependent boundary conditions. If there are time dependent boundary conditions which are known then you have to appropriately incorporate those at every time step and that is why it is okay to keep utilizing this time superscript always just to remind ourselves that there is a possibility of a time dependent boundary condition as well. So this entire equation which is basically your algorithm for the explicit method that is it. You are going to code only this in the sense that you will decide how many nodes you want to place in your entire y domain in that sense you are going to decide your own delta y you decide your own delta t how what sort of time steps I want to take to go from the 0th level solution to the next level and so on and so forth. Nu is going to be provided as part of the problem statement source is going to be provided as part of the problem statement. So entire RHS is known as I said and graphically if I want to show how the solution is progressing I am showing it at the bottom here where I have drawn the nth level and the n plus first level domain let us say at any grid point I what I am doing is I want to find out the solution at this grid point at the n plus first level. I am utilizing its value at the previous time level and the neighbor values at the previous time level if you see that that is exactly what is written on the right hand side and these arrows are supposed to simply say that this information is getting used to compute the value of u at this grid point for the next time level and so on for the entire set of interior grid points one by one. Very simple actually the explicit method is by far the most simple way of coding it. The only problem is that there is a so-called stability criterion and clearly I am not going to get into the discussion and derivation of what the stability ideas are. Very briefly unless you restrict the time step to be smaller than what is given by a criterion such as this what you will see is that the solution will not remain stable the numerical solution it will simply blow up and this is something that we normally ask our students to verify. It is very easy there is a criterion that is worked out based on some stability theory so as to say. So we will say that you know make sure that you violate that criterion when you code it and see what happens and clearly the solution is going all over the place and they get the feel that yeah it is really not working all right. So that is all about really the explicit scheme or the forward time central space scheme which is very very easy and very easy to program also there is hardly any programming effort involved all right. On the other hand we had another choice this backward time difference so let us see here backward time and central space. Let us see if we take that what happens. So this is the expression that I had written earlier u i n minus u i n minus 1 etc. So this n is really a dummy index in that sense so what I am doing is I am replacing n by n plus 1 if I do that I will simply write this as this will become n plus 1 this will become n this will become n plus 1 n plus 1 n plus 1 since anyway time is some in some way a dummy index. So I rewrite this in the form of u i n plus 1 minus u i n equal to all this in the n plus first time and then if I want to rewrite this equation in the form that we had written for the explicit scheme I realize that I have a slight problem to handle here because now I am realizing that I am interested in finding the solution at this grid point at the n plus first time level look how the information is going there for the calculation of the solution this is fine this is no problem this we have available from our previous time step what about these guys these are not available so the solution at this grid point is expressed in terms of the unknowns which are their neighbors at the same level so this is also not a problem so that is a slightly troublesome situation in the sense that we do not know this value we do not know this value so how do I get this no big deal what is done is you rewrite that equation again this one which is which is written out here and by rearranging the equation in the following form I write it as some coefficient times the u at the i minus first grid new time level plus some other coefficient times u at the ith grid point new level plus some other coefficient u i plus 1 new level equal to u i at the old level plus s times delta t which is something that we know from the previous time step so now let me not carry this constant in the form that is written out here you see that new delta t over delta y square it is coming same with the minus sign here and here so I am just going to call that as plus a for the sake of brevity so I essentially write this as a times u i minus 1 n plus 1 this entire constant I will call as b u i n plus 1 again plus a u i plus 1 n plus 1 equal to c i just for the sake of combining these two together at the nth level and that is it so what is this do this is what I would like to point out most of you would perhaps know this if you write this out for i equal to 2 3 4 5 6 7 8 9 all the way for each i you will see that you are writing this as an equation which involves three unknowns essentially the value at the node of interest the value of one neighbor to the left the value of one neighbor to the right and then there is a known constant let us see on to the right hand side and this way you can actually write the entire system if you want for all grid points one by one in other words if you write it you basically come to the conclusion that in order for all these u values to be found in this manner at the next time level I basically have to deal with the system of linear algebraic equations for the unknown interior grid points and basically the grid points of interest here which will constitute this linear system will be from i equal to 2 to n minus 1 why because the i equal to 1 and i equal to capital N are our boundary points where the values are known you can always include those as additional equations to complete this which is what I have done later but fundamentally this is a system of linear equations for the interior grid points. So explicitly if I want to write this I basically will realize that it can be written as a coefficient matrix times the u value that we want at all grid points at the new time level equal to this right hand side at the old time level which is all known. So in that sense whatever was written with this index notation can be written now as a matrix equation a times u equal to c where u is the column matrix you can you can say which is of unknowns. Now I am dropping the superscripts now I know that what I am looking for is the solution at the new time step having known the solution at the previous time step the solution at the new time step is basically going to be given by all these guys in the u column matrix. The right hand side column matrix is something that we know which is u i at the previous u at the previous time step plus s times delta t and this turns out to be a special form of the coefficient matrix where if you see I have added in the first row the boundary condition corresponding to the left node point which is u of 1 equal to 0 and at the last in the last row I have added the boundary condition corresponding to the right grid point which is u of capital N equal to capital U. So do you see the equivalence between these two the index form as it is written and the matrix form all that you have to do is really you actually have to explicitly evaluate each of these so 1 multiplied by u 1 plus 0 multiplied by u 2 plus 0 multiplied by u 3 so these are all 0 so what is left is nothing except 1 multiplied by u 1 equal to 0 which is our boundary condition fine. Then I will have for the second row a multiplied by u 1 plus b multiplied by u 2 plus a multiplied by u 3 equal to c 2 which is same as this written for i equal to 2 and so on is that fine and the last one again is all 0s not giving anything up to u of capital N minus 1 only final equation would be 1 times u n giving you capital U which is your other boundary condition. So for the purpose of clarity once I have shown this as a matrix representation where the matrix is in a situation that barring let us say the first row and the last row which I have just forcefully added as the boundary conditions everything comes in the so called tridiagonal form so there is only the main diagonal which has all those b coefficients the one which is to the left and the one which is to the right has coefficients a and this is quite obvious from your index equation which has been written out here. So to solve this what we do is every so the way to solve this then using this so called backward difference or the implicit method is that at time equal to 0 everything will be available as the initial condition. So initialize that at time equal to delta t I basically need to solve this linear system then I will obtain the solution at time equal to delta t fine then I go to time equal to delta t plus delta t which is say 2 delta t again I have to go back and solve this linear system. So for every time step that you take you have to solve a linear system like this where the right hand side will always correspond to the time level where the solution is known at the current time level and what you are solving for is this u1 sorry u1 and u1 are known u2 to un minus 1 at the new time. So how do you solve this most of you would know there is this so called TDMA tri diagonal matrix algorithm or the Thomas algorithm which is nothing but a Gauss elimination algorithm only for a tri diagonal system. So what ends up happening in what hand so most of you perhaps know Gauss elimination I am very sure in fact most of you probably know this also. So in a Gauss elimination what happens there are two steps one is what we call a forward elimination and then there is a back substitution. So the forward elimination does what what sort of a matrix form is brought about either an upper triangular or lower if you are doing it in a slightly different manner but let us say we are doing it in the standard way using which we will get an upper triangular form. Now here what will be the form of the upper triangular matrix in this case I would call it an upper bi diagonal right because what will be left at the end of the forward elimination process is the main diagonal and only this one which is above it that is it. So I will call that in this special case where the coefficient matrix is tri diagonal the Gauss elimination process actually will end up with a coefficient matrix which will be bi diagonal rather than calling it as formally upper triangular and that is what I have outlined here the algorithm really the forward elimination which will lead to an upper bi diagonal matrix and I do not want to explain the algorithm I want you to read and make sure that this is how it will be. So what you end up doing is actually since everything below the main diagonal is going to be made to zero you do not care about it all that you need to do is you need to store the values of the diagonal elements as they are changed and that is precisely how they are changed in a statement such as this. Correspondingly you will have to change the right hand side entries as well and those are also written in this form. Having done that the back substitution simply goes one after the other in the sense that the boundary condition is known anyway having known the boundary condition you simply work out one by one a previous equation to obtain the solution for all those previous grid points. So this is what the let us say the code statement that you may want to write in the code which will give you the back substitution. Having done this you basically have u of all i's that is going from n minus 1 to 2 at the new level and that is what your solution that you are after at the new time level. The you can see that there is lot more work involved in an implicit method because every time step you are required to solve that tri diagonal system fine but the advantage of the implicit method is that it is what is called as unconditionally stable. So it does not matter what sort of a time step that you take it will always give you a bounded solution. The solution will never explode just the way it did for the case of explicit situation and that is more or less it really what has been done in the in the program that you will run now is the problem that we talked about that transient problem in the infinite plate is solved using both the explicit method and the implicit method as I have been outlined here. So very quickly though it turns out that some of these books nowadays include some of this numerical methods discussion in particular this Potter's book actually has a fairly large CFD introduction where actually something very similar has been done a diffusion equation type situation in chapter 15 and Fox's book has a very small discussion which is actually not that that useful. There are obviously numerical methods book which will talk in far more detail so perhaps you should not pay any attention to these references much. What I will do is let me just project the problem statement that we have and I will just point out what has been done here. So as you can see it is the same problem that we talked about the transient problem in the infinite plate parallel plate etc the boundary conditions are as they are initial condition is known okay. Now the data is given so obviously for a numerical solution to work you have to provide numbers otherwise it will not work right. So I am providing some numbers this is kinematic viscosity density the capital U of 40 h of 4 centimeters and I am already giving looks like delta y of 0.001 has it been hard coded okay. So that is fine since the problem statement ask for it that is perfectly fine. Once you have been asked to look at the explicit method with a delta t of 0.002 seconds to compute the velocity profile for three different cases one is that no imposed axial pressure gradient then a adverse pressure gradient and later a favorable pressure gradient what you have asked you have been asked is plot the velocity profile at the time levels given by this I think if I remember you have done for more time levels than these and this is important superimpose the analytical solution which we worked out in the session before lunch at the same time level which essentially means that series which we are talking about earlier has been programmed to convergence so to say. So now that convergence criterion is what something that he has incorporated in terms of comparing successive terms you keep on adding until the difference between successive terms becomes less than some small tolerance value. And the same type of exercise is given as part two where the implicit scheme is required to be coded for essentially similar type of numbers is just that here you can see that since the delta t values immaterial since the implicit scheme is unconditionally stable you can actually utilize whichever way you want. The accuracy of course will depend on how fine a time step you choose but the solution is not going to explode even if you choose a fairly large time step. All right so specifically what we are looking for is second page of this third page of this document have already been given to you so just give me a few minutes do you have this open with you oh very good okay. So you can see figure 1.1 where what has been asked is numerically and analytically calculated velocity profiles at these various times using the fully explicit scheme for a equal to dpdx0 b dpdx of 20000 plus and 3 for minus 30000 and similarly for 1.2 for the implicit.