Lecture 9: Conservation Laws / Analysis / Shocks

Flash and JavaScript are required for this feature.

Download the video from iTunes U or the Internet Archive.

The following content is provided by MIT OpenCourseWare under a Creative Commons license. Additional information about our license and MIT OpenCourseWare in general is available at ocw.mit.edu.

PROFESSOR: All right, so I'll start with this review board, which shows the conservation law. Scalar, one unknown, one dimension, x, one space variable, and the flux function is f of u. That's the differential form. We have seen that at some points that differential form produces a discontinuous solution, a discontinuity. In that case we are better off with the integrated intgeral form. So I just integrate that one from a point x left to a point x right from an a to b and that gives me-- I have then the integral. So if u was density for example, the integral of density will be the mass. That's the quantity that I'm seeing here, so this is the rate of change of mass in the interval, comes from flux going out minus flux coming in. So that's the conservation law in its pure statement.

And we learned that solution for this-- while we're looking at the differential problem-- we can follow these characteristics from the initial time zero. So at time zero we look at a point x and the value at that point just travels-- is sort of safely carried up the characteristic line, which is that equals x0. And the problem arose when characteristics either collided into each other or separated. Those are the two difficulties and they have different resolutions. So this is what we met in the traffic example when a light goes red.

When a light goes red cars are coming up to the light, in front of the light the density is low. The cars have to brake quickly. So something happens-- that's a question mark there. What do you do when characteristics collide? When the value here cannot be both the initial value that's coming up that characteristic and the initial value that's coming up there because they're different? So we need a rule for what to do there. And we also need a rule for what to do in case of an empty space. When no characteristics-- here are two characteristics are hitting, here are none are touching there. Characteristics are separating in this case and the question is, what solution do you fill in the middle when you haven't got a characteristic? It's not coming from the start. OK, what solution to choose?

The answers will be different. And let me take the first one, the collision question, first. So what will happen is a shock will form. There'll be a shock line so that these characteristics carry information into the shocked, by the way, that's a big deal-- until they hit the shock and these characteristics carry their information. And across the shock we have different values, u on the left and u on the right. So that's one possible form of the solution because that will solve the equation in this region where the characteristics from the right are following the formula. It'll solve it from the left, so everything comes down to the jump condition at the shock.

What's the shock speed? Where does this shock go? It's to find that shock in the xt-plane that I have to use this integral form. I must use the integral form because I've got a discontinuity and the differential form just doesn't make sense. OK, I wrote down the key idea here or the conclusion, but not the reasoning. I'm looking for this jump condition at the shock, which is going to locate the shock. And the key quantity that I have to find is the shock speed, which will tell me where that shock is going, how it's moving up in time. So time is going that way, x is going this way in this picture. So this has a speed, s of t. And here's the conclusion. There's the jump condition. So what does it mean?

This notation, this bracket notation is the jump. The shock equation is beautiful. It says that s, the shock speed, times the jump in u is the jump in f of u. And let me take this example that is sort of the Burgers' equation example, so this is the example that comes from Burgers' equation. That's the example where the flux is this nice parabolic function u squared over 2. Then the change in flux, the jump in the flux function is u right squared minus u left squared over 2 and I divide by the jump in u, which is u right minus u left and that cancels beautifully and leaves me with the 1/2 of u right plus u left. It tells me that the shock speed is exactly halfway between, it's exactly the average in this Burger's example of the characteristic speed from the right than the characteristics speed from the left. And notice, it's in between. That's crucial. That tells us that this thing here is good. It tells us that characteristics are entering the shock from both sides. Characteristic lines are entering the shock. The shock is somehow in this sort of nonlinear case, it's somehow kind of cutting out bits of the solution, reducing the energy, reducing the entropy, reducing the total variation as we'll see it in a minute. So that's the answer here and I guess I have to justify, where did this come from?

Well, it came from this integral form. Suppose I integrate from a little bit to the left of the shock to a little bit to the right of the shock. Of course, that's where I want to integrate in order to capture the shock. OK, so if I just take that equation and write it over here and I'll imagine-- because I'm not going to be x left and x right are going to be pretty near the shock-- u up to the shock will pretty much have a constant value, u left. So approximately that is the shock position. Let me say that shock is at captial X is the position of the shock. So I have x minus x left times u on the left of the shock. I'm doing that integral. Oh, I need d by dt of all this. So d by dt. It's d by dt of the integral of u dx. Let me just put here what I'm computing. I'm computing that from x left to x right.

Well, I'm going up to x, up to capital X where the shock is and then onward. So here I went up to x using the value u left. And now will be the onward part, which will be x right minus x, the half of the interval that's on the right of the shock times u right. And then copying this term is exactly our jump in f of u. Plus the jump in f of u is zero. Those words, total variation are a reminder to myself for what comes next. OK, so is that OK? That will be approximately zero because u wasn't exactly constant on the left and right. It didn't stay at the value ul, but we're squeezing the interval down where it's essentially constant. OK, so now just look and see what we've got. dx dt is the shock speed, so that's s. Multiplying ul and here I have the derivative of that is zero, the derivative of xr is zero. Here I have a minus, so it's minus ur. dx dt times that number is e plus the jump in f of u, left to right, is 0. And that's our equation. If I put this on the opposite side it becomes u right minus u left. Shall I do that?

So I'll keep this on this side, got a good sign, but I'll move this onto the other side so it will be the shock speed, d capital X dt times ur minus ul, which is the jump in u. So that's great. We now have a way to deal with shocks that collide. So that was a derivation of this jump condition, which first had a name, two authors [? Rankin Whogonio ?] in gas dynamics, but then it was seen as the right condition in all conservation laws. OK and now I'm ready to comment on the second issue. What do you do when characteristics separate?

Well, let me say what you don't do. You could, without wrecking the equation, you could try to put a shock in there and give this the value, u on the left even though a characteristic isn't telling us what's going on in this space or this one. You could try to put a shock in there. You could maybe satisfy the jump condition if you put the right shock in, but the characteristics would be coming out and that's wrong. And it's called the entropy condition, which says that this can't happen. That the derivative of entropy has a certain sign and it rules this out. It rules out this possibility of characteristics emerging from the shock and we need a totally different solution. And it turns out to be a fan in there. So we put in a fan of characteristics coming from here and that gives a solution that actually will depend on x/t. This solution will be a simple ratio, x/t enters that.

In other words, that gives us a smooth solution and I guess thinking about the traffic example, that's exactly what happens. So this is the shock part, now for the fan part, what happens when we have a bunch of cars sitting there, ready to go, the light goes green, the first cars start moving, cars behind them see the car in front moving, they start moving and so the initial profile would be, in this case, for cars, would be high density and then zero density. So this would be maximum density sitting at the light and zero density in front of the light and then the light goes green now and what happens? Well, as we all know those first cars pull ahead and the profile is these cars haven't yet started to move. These have started to move and there's a profile there and then zero density. So this represents the first car, the lead car. In front of that car is zero density and that's the fan. That's the x/t expression. So a sample linear expression, which will flatten out as time goes on, these cars will be going faster and faster and then eventually up to full speed. OK, so that's the fan. Maybe now is my chance to speak about this concept of the variation in a function. The idea will be you now, how do we find the difference method to do this automatically?

Well, OK, one way to do it is put in a little viscosity. This rule, the rule that picks out this shock in the one case of colliding characteristics and this fan in the case of separating characteristics, one way to get to those two solution would be to put in a little bit of viscosity and then let epsilon go to zero. And the result would be these two solutions. And in a way that's what finite differences is going to do. Finite difference or finite volume is going to be a numerical viscosity and the epsilon in the Burger's equation would become delta t or a delta x in the numerical method in the finite difference equation. And then what's the concept? So we have to find finite differences. And we want those-- the point is finite differences, sort of takes more care. We can't just create something that's consistent in the smooth case.

Up till now we've just had no hesitation. We replaced derivatives by differences without too many worries, really. And that was fine as long as the solution is smooth. But now if we're going to replace derivatives by differences we have to do it in such a way that this entropy condition, these solutions emerge automatically. And really, in a way that our finite difference or our finite volume methods should somehow conserve mass the same way the true equation does. So those are two different points. First point is this question of total variation. What's the total variation in a function?

Well, let me draw a function. The total variation is-- well, let me see. Total variation in a function u is going to be, by definition, the integral of the absolute value of du/dx. OK now we have to understand this concept. It's beautiful that the analysis has led to a quantity as simple as that, as a guide. So when we create finite difference methods we will ask, are they total variation diminishing? So TVD will be total variation diminishing methods. Methods where this quantity does not increase in time. Let me just understand this quantity first. The absolute value is crucial. If the absolute value was not there then I would just have the integral of the derivative, so it would be u at that point minus u at that. But here I'm taking absolute values, so for awhile the derivative is positive up to there. So the integral gives me this value minus this. Then in this region the derivative is negative, but I'm taking absolute value. Do you see that that integral will split into 4 separate integrals? The total variation will be this difference plus this difference plus this plus this. I'll have 4-- in this example-- have 4 positive contributions to the variation. And that's the key quantity there. And what a shock will do-- can I maybe mention in this figure what a shock does, can do, would be to somehow cut out part of the variation. When a shock occurs it might happen that these things are going into the shock, these are going in from that side, and the shock kind of jumps from this value, u left, to this value, u right following the shock condition when it emerges. And this additional variation is wiped out. You see that the variation for the solution that contains a shock is still-- this part is still there, but this part is all downhill, so it's just this. It's that minus that. And the middle two parts have been cut out. OK, so that's the total variation in a function.

Let me just take the most important example that's not total variation diminishing, which is Lax-Wendroff. The one even in the linear case. The one experiment we've done with that profile, if you remember we had a wall of water and the exact solution brought the wall to the left with velocity c, with speed c, but Lax-Wendroff, which did quite well staying near the discontinuity near the jump had a little oscillation behind. That was essentially our one objection to Lax-Wendroff was that oscillation showing up. And so that's telling us that the variation--see, the variation in the true solution is just that step. The variation in the Lax-Wendroff answer is this bigger step plus this one plus this one plus this, larger than 1. Roughly the match is that monotone methods. So here's the general picture.

We can get monotone or I'll say, all positive coefficients and those will be total variation diminishing. Great. Why don't I draw-- recall a monotome example was Lax-Friedrichs or upwind. Upwind will be our model, so let me draw the upwind one. The upwind one didn't increase the variation, it stayed at 1. We had no oscillations, but it was low accuracy. So first order. We can't do better then first order accuracy, the order of accuracy I'm measuring in the smooth part of the solution because that's where I can take taylor series and match terms and see which is the first one that doesn't match. And Lax-Friedrichs and upwind was first order. Then our problem is that if I allow a negative coefficient then variation increases, ocsillations will come in, but I can get the accuracy up to second order or higher. And Lax-Wendroff is the example, the key example there that I'll continue with. OK, that's the idea of total variation. It just gives us a quantity to work with when we see oscillation. And of course, I'm interested in oscillations, so this is a linear case with a discontinuity. I don't think it's quite fair to call it a shock. At least, the shocks that we meet in conservation laws in the nonlinear problems, they're created-- the ones we saw here are created for nonlinear reasons. Colliding characteristics, whereas the linear case has characteristics all parallel. So here we're in the lecture, we're really into the second lecture into the scientific computing part of our problem. So that's the analysis part, the pd part. Now comes the experience of quite a lot of people working hard to identify a way to get second order accuracy where it smooths and monotonicity or something like it where the shocks are. So that's essentially the goal. Our method is going to turn out to be nonlinear. It's going to have some-- it'll be called a flux limiter. The key quantity will be the flux. And for a smooth function when u right and u left are very close, the flux is small, we don't have a problem. But at a shock where there's a significant jump in u and especially in f of u, the flux limiter is going to come in, into the discrete method. Let me began discrete methods with a new word, finite volumes. And with a new interpretation of capital U.

So the discreet value, capital U jn, which up to now we've been thinking of was that's an approximation to the true solution, j delta xn delta t. Now because we're dealing with discontinuities we better integrate. We integrate over a delta x interval. We take the average over a little delta x interval instead of hitting just the point. So the delta x interval from j minus 1/2 to j plus 1/2, we take the average of u at the time, n and delat t and that's what finite volume-- and let's think of as capital U, the discrete approximation. The appropriateness of that is in the integral form of the conservation law. So the integral form of a conservation law leads us to that natural discretization. d by dt of u, of the average. You're keeping your eye on the integral form of the conservation law. Just come back to it again. The derivative of the mass, the integral of u, so in the discrete case the difference of capital U plus the integrated term, a flux at one end of the integral minus a flux at the other end. And notice that we are getting for the flux, 1/2. So a staggered grid is presenting itself naturally. And so that f j plus 1/2, the fluxes, which is the critical thing is on the j plus 1/2 grid.

All right, now this could be exactly defined as the flux at that grid point. Capital U squared over 2 in the Burger equation at that grid point. That would give us a totally straightforward, first order numerical method. First order, I see first order here in time and it might be better than that in space because in some way that's centered. But if we want higher accuracy, if we want to include Lax-Wendroff for example, that includes several values of u, we go to a more general flux function. The whole job here is create a smart flux function. So a flux function depends not only on u in the center, uj-- at time n, this is all at time n. It can also depend on values to the left of it and values to the right. It may use those values, it may not, but we allow that possibility then that the flux-- so this is the numerical method that you aim for. Is to choose a flux function and of course, this should be consistent with the original problem and that simply says that if all the u's were the same that the numerical flux function would give you the flux function in the given equation. OK, so this is our general pattern. And now the question is really, maybe the first step is I would like to write upwind and Lax-Wendroff. I'd like to find their flux functions. I guess I'm saying that when we look at an equation in this format, we can take methods that we have already used and see, OK, even in the linear case, the one-way wave equation-- what are we picking for the flux function?

OK, now the one distinction is-- the one-way wave equation we had a constant and upwind for example, went in the direction decided by that speed, c. If the speed reversed sign so that the waves were going the other way, the wind is coming the other way, the upwind direction is the other way so our method would change. Let me just write down the flux function that turns out to come for upwind and then for Lax-Wendroff. OK, so for upwind the flux turns out to be f, turns out to be 1/2 of-- because it's at this, this is fj plus 1/2, really. It's at the center point. The flux function will be, I'll put in the a. a delta t-- sorry, a over 2uj plus 1 plus uj. Now plus, or rather minus actually, comes absolute value of a/2 times uj plus 1 minus uj. If we're patient we can track down the two cases, a positive and a negative, and recognize that this form switches direction correctly. The difference method. So it keeps us upwind in other words. If the wind changes direction our flux changes appropriately to follow that direction. So that's an OK flux function. It's a monotone one and it's first order.

Now let me put in Lax-Wendroff. What would Lax-Wendroff be? Would be the same, the same upwind flux, if I call that up for upwind and then a correction, a Lax-Wendroff term that gets the second order accuracy. OK, and that turns out to be, it turns out that I multiply this term by that ratio, r. By the absolute value of r. So it turns out that it's the Lax-Wendroff thing, to get that extra accuracy. You might remember, had an extra delta t over delta x. And so it's this same thing times a over 2. In the nonlinear case it'll be f of uj plus 1. f at uj j plus 1 minus f at uj. Next time I'll come back to this point and pin down this more carefully. All I want to do in completing today is introduce this idea of a flux limiter.

So what people say when they look at Lax-Wendroff is that this term that's been added on is fine in the case of a smooth solution, but in that case of a jump like our picture here, that term is increasing variation. It's an anti-diffusion term. Instead of smoothing things out the way f upwind does, that extra term has the effect of increasing the jump. It's anti-diffusive and therefore when it's a significant term the flux limiter cuts back on it. The result will be a numerical method, which is not linear anymore even if we apply it to a linear problem. The result won't be strictly linear and it will produce a profile that doesn't have that oscillation. So it has what we want. It gives us good answers when the solution is smooth, it stays near the shock the way Lax-Wendroff does, but it doesn't oscillate.

OK, so what's the flux limiter? The flux limiter is, I multiply in here by some function, the limiter function, phi j plus 1 and I multiply in here by phi j. In other words, there's an extra factor, this flux limiting factor. So the question is, how to construct that? What formula do we use for the flux limiter? So the flux limiter, after a lot of experiments is chosen to depend on-- so we want it to be 1. So what do we want? We want it to be 1 or close to 1 and smooth part and we want it to be greater than 1 at a discontinuity. OK, so our question is, how do you identify that discontinuity?

Well, the convenient way has turned out to be look at the ratio-- look at differences. Look at the ratio of-- so what comes in here will be the ratio r say j plus 1/2. That'll be the ratio of u, say j plus 2 minus uj plus 1 to uj plus 1 minus uj. All that time n. That ratio tells us, that ratio is close to 1. So the ratio r is close to 1 for smooth and our function, phi of r then, we want our function phi when the ratio is 1 to give the value 1. So if I try to draw a graph of-- let me draw a graph of peoples flux functions-- of the phi of r. So here's a graph of phi of r.

First on the graph I'll put upwind and Lax-Wendroff without anything. So Lax-Wendroff would take-- phi Lax-Wendroff is identically 1. This is the ratio r and this is phi. So Lax-Wendroff before fixing had no flux limiter function and it just had 1. And actually, if I can keep the same formula, phi upwind would just be zero. Then it'll turn out that-- let me draw these ratios. There is a region here and here where the flux limiter is doing the right thing. So in that region it's acting the right way. It's limiting the anti-diffusive flux. It's controlling the variation, it's controlling the oscillation. So various people have proposed-- let's see, no. I haven't got this right. I think it's- because I know that part is wrong for Lax-Wendroff. So maybe it's there. No, I've forgotten-- sorry that I'm not, I've lost this because I know that one possibility is this as the flux limiter. One possibility would be this one and then there's another possibility that is OK that goes through here. Let me stop for today with this idea of a flux limiter and you see that it's a task to fix the method to do the right thing at smooth parts and also the right thing at discontinuities. OK, I'll pick up on that next time. Thanks.

Free Downloads

Video


Audio


Subtitle

  • English - US (SRT)