 In today's class, we will discuss about hydraulic flood routing. In earlier classes, you have already been told about the flood routing by way of hydrologic methods, but in today's class, we will discuss about the hydraulic way of flood routing. So let us say this is a channel, the flood wave is moving in this direction and here I know what is the hydrograph. You all know what is a hydrograph, right? This is T, this is the discharge. So the time history of the discharge will be called the hydrograph. So in this channel system, I know the inflow hydrograph and here I know the rating curve. Let us say this is the downstream end of the channel. So here I know the rating curve, the stage discharge relationship. Given these conditions, I want to know the flow depth and the velocity or discharge at any other location. For example, here I want to know for a given time what is the discharge or flow velocity and what is the flow depth. So how we will go about this? When you say hydraulic routing, it can be done in many ways. So to start with, let us begin with dynamic routing. In dynamic flood routing, we use the governing equations as continuity equation for water flow and then momentum equation for water flow. If I consider the flow to be one dimensional, then what should I do? I should use these equations. The continuity equation will be partial h over partial t plus partial u h over partial x is equal to 0. Please remember that we have assumed here the flow to be one dimensional. That means suppose this is the cross section of a river. At any point, this will have three different components. The flow will be in three different directions. Let us say this is along the flow, this is along the width and this is along the depth. So when I say the flow is one dimensional, that means flow along the main direction or the longitudinal direction is prominent. Other two directions are not so significant. So I consider that the effect of these two components of the flow velocity are minimal. So I consider this direction only. So when I say this is one dimensional flow, this h is basically the flow depth and u is the flow velocity. Do not forget that we have also assumed that the bed is rigid. You know that all natural flow channels are with sediment flow but here for the sake of simplicity, we have assumed that the bed is rigid. That means no sediment flow. So this is the continuity equation for water flow and the momentum equation for water flow will be partial u h over partial t plus partial over partial x u square h g h square by 2 is equal to g h s o s f. Coming to different terms in these equations, as we said earlier u and h are respectively the flow velocity and flow depth. Here g is the acceleration due to gravity, s naught is the bed slope and as we are assuming one dimensional flow, s naught or the bed slope will be the slope along the flow. That means if this is the channel, this slope or in other words if I consider this to be x, this to be z, slope will be minus del z del x. So this is the slope or the bed slope. S f is the friction slope. What is friction slope? Suppose this is the water surface. So we consider three different types of slopes. One is the bed slope, which is the slope of this. One is the water surface slope. The slope of this line. Another is the friction slope or the slope of the energy line. I can draw the energy line here like this. So slope of the energy line is called the S f. If you consider the equations mathematically, these equations are hyperbolic partial differential equations. These are non-linear by nature. So a closed form analytical solution is very difficult. We can find closed form analytical solutions only for idealized cases. Therefore, one should go for the numerical solution. So we will study here the numerical solution to these equations, but before that let me explain you this S f. How one will estimate the S f? S f is the friction slope and it should be estimated by some resistance relationship. For example, Manning's equation. So one can use Stege's equation or Manning's equation, some sort of resistance relationship and you must know that these resistance relationships are valid for uniform flow. For example, Manning's equation is valid for uniform flow. That means we assume in our analysis that for unsteady flow, even the uniform flow resistance relationship is valid. This is an assumption. Strictly speaking, this will not be the case. In unsteady flow, the resistance relationship will be different. However, we assume that Manning's equation is equally valid in case of unsteady flow also. So coming back to the numerical solution, if I have to write some of the names here, then one will be method of characteristics. Then there is finite difference method. Then we have finite element method, also spectral methods. There are so many variations of these methods also. However, our discussion is limited to this finite difference method. You know this method of characteristics. This is generally used in pipe flows, but in channel flows, we generally use finite difference method. So let us consider about finite difference method. I will explain here one explicit method and if time permits, one implicit method. What is the difference between explicit and implicit method? In explicit method, suppose this is a partial derivative term. In the finite difference approach, what I need to do is I replace this partial derivative term with equivalent algebraic term and by doing so, I convert these partial differential equations to algebraic equations. Now, I have set of algebraic equations and unknowns. I try to solve the unknowns. Suppose I have a partial derivative term like this. So the spatial derivative will be expressed in terms of unknown time level values. For example, suppose I use a forward finite difference. So this is the forward finite difference, but look at the time levels. The time levels here are for unknown time level. This is called implicit and in explicit, we use by using same forward finite difference, but the time level will be for known time level. Here, these symbols k and k plus 1, they represent the known time level and the unknown time level. It will be more clear when we draw the grid. So let us say this is the finite difference grid. This is the distance and this is the time. To represent distance, I use the nodes as i's. So these are different nodes. That means different i values. They represent the distances and these are different time levels. Let us say I represent them through k. So this grid indicates the spatial variation and the temporal variation of our problem document. This is the spatial direction. This is the temporal direction. So this x means, let us say this is 0 kilometer. This might be 2 kilometers. This might be 4 kilometers like this, but when I say i, this designation i means some number. For example, this might be 1, 2, 3, 4, 5. Some nodes, I designate these nodes through these i values. Similarly, the time nodes are designated by k values. For example, k is 1 here, 2 here, 3 here, 4 here. In absolute sense, this k is equal to 1 may represent, let us say time is equal to 2 hours. So if I am interested to know the variables throughout the channel, do not forget this. This is the channel. These are the different nodes which I represent in this grid as i's and these nodes on these nodes both h and u. If you remember the governing equations as we wrote for the continuity equation and the momentum equation, the independent variables are x and t, the space and the time nodes. But the dependent variables, the unknowns are the flow variables u and h, the velocity and the flow depth. So here, at all the nodes, I need to find out what is h and what is u. On all the nodes, here I need to know what is h and what is u, what is h and what is u, what is h and what is u like this. So for all the nodes and also at all the time levels, I need to know how h and u vary. So h is a function of x and t, h varies with distance as well as with time, u also varies with distance and time. So how do I define the problem? So this flood routing problem I will define by this. Suppose this is the channel, I know at a given time the flow depth and the velocities at all the nodes. So this is called the initial condition. At t is equal to 0 for all x, I know u for all x is known and h for all x is known. So initial condition says at time t is equal to 0, I know the variables values at all the nodes. What is boundary condition? Boundary condition is at unknown time levels, it is greater than 0. I know the values at the boundaries or at the end nodes. As you know, the channel is divided into number of nodes. So this first node and the last node, the inflow section and the outflow section on these nodes, I should know the variables values. So this is one section, this is another section, this is called the inflow or you can say the starting point of my problem domain, this is the outflow section. So I know some condition here, I know some other condition here. From these given conditions called the boundary conditions, my results should be affected by these conditions. So let us once again formulate the problem and we are discussing one numerical method which is explicit and the name of the method is Lax Diffusive Scheme. So we are discussing about Lax Diffusive Scheme. I am not repeating the governing equations. You are supposed to know these governing equations. Few minutes back we discussed that. So here the main philosophy is to write the time derivative. Let us say this is a general type of time derivative. It should be represented through, I will come to this star business this. Look at this expression very carefully. This partial derivative term is made equivalent to this algebraic term and this is for node i. That means if you remember my finite difference grid, let us say this is a particular i and I am trying to evaluate the time derivative term at this node. For this particular node, I can write these equivalent algebraic terms and this f is the functional value. If this is h, this will be h. If this is velocity, this will be velocity at i at an unknown time level. That means if this is let us say this is k, this is k plus 1. If this is i, so this is unknown and this is known. Everything here is known. At this level I do not know anything. So my objective is to find out u and h values at these nodes. So at node i del f over del t is made equivalent like this. What is f star? This is the average of the neighborhood values. Mathematically it can be proved that if I use the value f i k here, the scheme will be unstable. What is an unstable scheme? An unstable scheme is when the numerical errors, total numerical errors, they grow with time stepping. As time increases, the total numerical error, they increase. So we do not want that. Our scheme should be stable. That means the error should be within some limits. To fulfill that condition, we use here a averaging parameter which is because ours is a case of one dimensional problem. So it will be f i minus 1 plus f i plus 1. And do not forget this is for the known time level. That means f i star here will be average of these two values. This is i minus 1, this is i plus 1. In case of a two dimensional problem, this should be made average of all the four neighborhoods. However, this is a one dimensional problem and we are writing this star parameter with two neighborhood parameters. Now this is the substitution for the time derivative. What about the space derivative? If this is a general statement for a space derivative, then I can write the equivalent algebraic expressions. This will be f i plus 1, f i minus 1 divided by 2 delta x. We use a central finite difference in case of spatial derivative. Do not forget to use the known time level values here. This is for node i. Once again I repeat, when we look the governing equations, we see that in the governing equations, there are spatial derivatives and there are time derivatives. The time derivative terms are made equivalent to algebraic expressions by way of time stepping in which we use the averaging of the parameters between the neighborhood values. However, in case of spatial derivative, we use central finite difference and because lax diffusive scheme is an explicit scheme, so spatial derivative will be expressed in terms of known time level values. That means, here we use known time level values. So, I will just give an example how to discretize, how to use these principles in the governing equations. If I rewrite my governing equation, let us say continuity equation. This is the continuity equation. So, in this equation, I have to substitute this for the time derivative which will be h k plus 1 i minus h i star h i star divided by delta t and for this, I will have h u i plus 1 h i plus 1 minus u i h i. This should be at known time levels divided by 2 delta x. I substitute the spatial derivative here and as we discussed for the time derivative here, we will be equal to 0. So, if I look at this discretized equation, this is called the discretized equation. This is the partial differential equation. This is the equivalent algebraic equation which is made, which is in the discretized form. This is called discretization. So, here again h i star, if I replace h i star as we discussed, it will be average of the neighborhood values. So, it will be half of h i minus 1 k and h i plus 1 k. So, if you look at this equation, everywhere whenever there is the superscript k, I know the values because in the numerical grid, k indicates known time level and k plus 1 is the unknown time level. That means, if I look very carefully this expression, I have only one unknown time level expression. So, I can write h i k plus 1 equal to minus here half. This will be u i plus 1 h i plus 1 u i h i. This is with k. This is with k divided by 2 delta x. Then delta t will come here multiplied by delta t and this will be plus half plus h i minus 1 h i plus 1. This is k. This is 2. So, on the right hand side, we have all the values with superscript k. That means, I know these values. So, the beauty of this explicit scheme is that we can know the unknown time levels value with the known time level values and I do not have to solve any matrix or any difficult thing. Explicitly, I can have the terms for each node. I can write this term for all the i values, but just hold on. When I say all the i values, if I draw the grid, this is i and this is k and if you relook the expression we wrote there, it is involving i plus 1 and i minus 1. That means, if this is my i, let us say this is i, it will involve evaluation of h here will involve this and this. So, here we have a problem and what is the problem? The problem is these two values will not exist for the end nodes. For example, if I am here, I can have this, but there is no node here. Similarly, for the downstream end, I do not have the node here. So, what we do is this lag scheme can be used. Just remember, lag scheme can be used for the interior nodes and for the end nodes or the boundary nodes, we use the boundary condition. What will be the boundary condition? The boundary conditions will depend on the problem statement. When we say problem statement, in some cases it will be the inflow hydrograph for the inflow section. So, you remember these statements, let us say inflow. At inflow section, I must use, of course we have not discussed the method of characteristics, but for the time being, just take it for granted that there is something called method of characteristics and in method of characteristics, we get one negative characteristic equation. So, for inflow section, one should use the negative characteristic equation and the condition from the problem statement. So, when we say from the problem statement, for example, if this is the channel and you have been provided with the inflow hydrograph, then this is your problem statement. So, one equation is this, the other equation is this. We have two equations and two unknowns in H and U. So, one should be able to find out both uniquely for the inflow section. One simplified method one can use in place of this characteristic equation is, you use one equation from this and you use. So, I am writing here 1 a in place of this negative characteristic equation. Suppose you do not know what is negative characteristic equation, you can go for this method. This method says use the extrapolation. What is extrapolation? Extrapolation is suppose these are the nodes. I know the values, let us say H, I know the H value here, I know the H value here. So, you extend this line to find out H value here. That means H, let us say this is 1, this is 2, this is 3. H 2 k plus 1 is known. How it will be known? This will be known by the Lax Diffusive Scheme. When we use the Lax Diffusive Scheme for the interior nodes, I can find out H and U for these nodes. So, this is known. Similarly, H 3 at unknown time level is known and what about H 1 k plus 1? As we said, we cannot use the Lax Scheme for this. So, we will use the boundary condition. So, this will be expressed in terms of these two. So, if I use a linear curve, that means I assume that the relationship between these three is a straight line. I can fit a straight line between the H values at these three nodes. So, it will become this plus this divided by 2 will become this or 2 times H 2 minus H 3 k plus 1. This is called the method of extrapolation. So, you use extrapolation method in place of method of characteristics and do not forget the condition from the problem statement. For example, if you have been provided with the inflow hydrograph, then use that condition and also this condition. These two will be sufficient to find out what is the H and what is U at node number 1. Similarly, for the end node or we can say the downstream end. Generally, we have the rating curve and we have the either extrapolation or positive characteristic equation. Now, we will discuss the flow chart of the Lax diffusive scheme so that the sequences, the sequential steps to be used in Lax scheme will be clear. So, we start the program here, then all the inflow values sorry input values. When I say input, what are the different input values? It will be different channel parameters such as slope length, its cross section at different places and then we have the numerical parameters then you may have the flow parameters and you may have the constants. For example, we have G in our expression. So, acceleration due to gravity G is equal to 9.81 that information you should put as input. Then you set the initial condition that means at t is equal to 0, you set U for all i's and H for all i's. Then from this now you have to find out what is your delta t. I will write your stability. I have not explained what is stability criteria. I will explain after this because this will be a breakage in the continuity. So, first let us complete this flow chart and then we will discuss about what is stability. For the time being just know that stability means we find out what is delta t. Unlike in case of implicit schemes, we have to follow certain conditions for delta t. In explicit schemes, although it looks very simple, we are paying the price here. For example, we have to follow the current condition to find out what is the delta t value. We cannot use delta t values at our own sweet will. So, after stability, one should go for the unsteady flow condition and here you can solve the H equation or the continuity equation. Solve continuity equation and solve momentum equation. When you solve continuity equation, you will find H k plus 1 for all i's except 1 and the last node. Similarly, here you will find u for all i values, but you cannot get these values for node number or i number 1 and i number last. So, what we do is we use here the boundary condition. Please remember that these two lines are from the lack scheme. The last line the boundary condition is not solved by lack scheme. You have to use the problem statement and common sense and general rules of the numerical method to find out what is the boundary condition. Then we should go for the next time step, but before that I should verify this diamond shape indicates a logical statement. So, I should verify whether I have reached the required time level or not. So, t here is verified. If I have already reached the required time level, then I should stop. If I have not reached the required time level, I should repeat, but before repeating I should swap the I should exchange because if this is the grid, let us say this is my k, this is my k plus 1. That means this is the known time level, this is the unknown time level. So, using these values, I calculate these values. So, here at the end of this, I am here. Now all the values are known. So, I should treat this as k. That means now this becomes k and this becomes k plus 1. So, I should swap the variables. So, here I should do h k is equal to h k plus 1 and u k is equal to u k plus 1. Then, once again it becomes the same thing. That means now I am at time level k, which is known for example here and I move forward, but before going there I should know what is stability, what is delta t. Then I solve the lax diffusive scheme, put the boundary conditions and get the answers. So, I time step this. That means I go to one time level, another time level, another time level like this and I will stop when I reach the required time level. For example, suppose there is a flood in a river and I know that the flood takes, let us say 4 hours to reach the downstream end. Then probably I should do the flood routing for 4 hours or maybe 10 hours, something like that. So, the user, the programmer has to give the time of computation. That means at which time level I want to stop my program. So, that has to be given. That will be given in this numerical parameters in the input file. So, you use the time of computation here and as soon as the program reaches the time of computation, it will stop. If I want to write the results, then where should I write, I should use, it depends. For example, let us say this is the channel. One sort of output will be, suppose I am at a place for a given x, I want the time variation of the hydrograph. Let us say, suppose this is h, this is t. Let us say this is the flow depth, this is time and this is like this. To get the result like this, what should I do? Where should I write the output? I should write the output here. Before I move here, this should be my output. So, here I should write t and also h k for i. This i is fixed. This is for a particular given. This is for some value, let us say 50. For a given distance, I have to know for each time level, what will be the variation of depth. If my objective is that for the whole domain, I should know how the surface look like. For example, suppose this is the water surface, this is the bed. So, here the objective is different. This is x, this is h. You compare this and this. Here the objective is different, here the objective is different. So, for this, what I need to do? After the end of all the computations, before stopping, I should write the output here. So, depending on the objective of the user, depending on the objective of the user, one should use the output statement. Now, let us go to the stability, because we have not discussed the stability. So, we will discuss the stability here or in other words, if we are using an explicit scheme, then we have to use a limited value of delta t. We cannot use any value of delta t. So, how to decide that? This is called Courant-Frederick-Lewis condition, C-F-L condition. This condition says delta t should be less than equal to delta x divided by u i plus minus square root g h i. So, we multiply here. How the computer will know that this is less than this? Delta t is less than this. So, to ensure that, we use a factor which is less than 1. So, we call that as the Courant number. We use a Courant number. That means delta t is equal to Courant number times delta x divided by u i plus. Let us use plus, because this will give you a stringent value square root g h i. Let me explain this. This Courant number C n is a factor which is less than 1 and this delta x will be given by the user. u i and h i are known. g is acceleration due to gravity. Based on this, delta t will be calculated. A typical value of Courant number for lag scheme is in the range of 0.6 to 0.8. Even higher values, one can numerically experiment what is the value of Courant number and accordingly, you can decide what value of Courant number is appropriate for the given problem. Once you find out delta t, please check that this is dependent on i. That means, if this is the numerical grid for each time level, let us say this is the time level k plus 1. So, to determine what will be this value of delta t, I need to find out delta t for each i and I should select the minimum of all the delta t's. So, you find out minimum of delta t i. Let us say this is i. So, for each i, I will find one delta t value and based on the minimum criteria, I will choose the minimum value of all the delta t's. So, I go for the minimum of the delta t to select what will be my next time level. So, once delta t is calculated, then we can follow this flow chart to find out this is the stability step. So, we can move further. Please remember that this is not uniform. That means, this delta t and this delta t, this delta t, these might be different. Delta t may not be equal for all the time steps. Now, when you do the lag scheme, what is expected? Sometimes, we find in flood routing some discontinuities like let us say something like this. Let us say this is x, this is h. In flood routing in hydraulics, this is called a bore. Mathematically, you can say that this is a discontinuity. So, lag scheme is not very good to capture the discontinuity like this. That means, strictly speaking, the numerical scheme should produce these as the result, but lag scheme will give what? It will give you something like this. It will give you some sort of the nature. It will indicate you the nature, but step like things will come. So, when we get step like things, we know that this is not a very good thing. So, what we do? We do something called the, this is called a numerical error. So, because lag scheme is only second order accurate in space and first order accurate in time, we find these type of errors. So, there is another scheme which is called MacCormack scheme. This, if we use this scheme, then the step like things will not come and it can capture the bore very correctly. So, as I mentioned in the MacCormack scheme, the philosophy will be slightly different, but it is similar to lag scheme. So, this is a two-step scheme. One is the predictor scheme. The other is the corrector step. In the predictor step, this is calculated by f, let us say P fk i divided by delta T and we also have the spatial derivative terms, which will be in terms of known time level divided by delta x. This is backward finite difference and this is time stepping. This is the predictor step and in the corrector step, we will have, let us say c is for corrector step. Now, this is known the predictor step. So, this is used and for the spatial derivative, we will use the predicted values, but here we have used the backward finite difference scheme. So, here in the corrector step, we will use the forward finite scheme. This is in terms of predicted values and the thing is not over here. There will be a final step, which will say what is k plus 1. That means, is half of f. Remember this very carefully. This is corrector step. This is not predictor. This is at time level k. So, in the scheme of the parameters, the variables are defined in terms of, let us say h and u at time level. I have h k and u k. These two are known. First, I will find out what is predicted value of h. Similarly, what is predicted values, value of u. From this, I will find what is corrector value of h and then, I will find what is this. This is the sequence. So, in McCormack scheme also, it is not free from numerical error. Any numerical scheme you use, it will be with some error. So, here we find oscillation. For example, if this is the shock, we will find some oscillations here. So, this is also not very good. So, what is done is there is something called artificial viscosity. This is used to kill the numerical oscillations. So, by using artificial viscosity and McCormack scheme, one can capture the shocks very well. So, to conclude today's lecture, I will say that this hydraulic way of doing the flood routine is very good. In this, we use the continuity equation and the momentum equation, which are also called the Saint-Venant equations. We solve these equations by some numerical method. We have discussed two explicit methods. One is the Lax diffusive scheme. The other is the McCormack scheme. As I told, McCormack scheme is better because it is fourth order accurate and this scheme captures the shock very well. So, we will stop here.