Proof of Harten’s Lemma re the Convergence of TVNI Finite-Difference Schemes

An academic paper with a rather unusual history is that of the Israeli mathematician Ami Harten’s “High Resolution Schemes for Hyperbolic Conservation Laws.”  First published in the Journal of Computation Physics in 1983, it was republished in 1997 in the same journal, and is often cited with the later date.

At the time of republishing Peter Lax, who earlier had saved the computer from the hippie radicals, made the following statement about this paper in an introduction:

This paper was a landmark; it introduced a new design principle—total variation diminishing schemes—that led, in Harten’s hands, and subsequently in the hands of others, to an efficient, robust, highly accurate class of schemes for shock capturing free of oscillations. The citation index lists 429 references to it, not only in journals of numerical analysis and computational fluid dynamics, but also in journals devoted to mechanical engineering, astronautics, astrophysics, geophysics, nuclear science and technology, spacecraft and rockets, plasma physics, sound and vibration, aerothermodynamics, hydraulics, turbo and jet engines, and computer vision and imaging.

One point in the paper was a lemma concerning the convergence of TVNI (total variation nonincreasing) finite difference schemes.  Concerning the name of these schemes, Lax points out the following:

Harten originally called his schemes variation diminishing, abbreviated TVD; when Osher pointed out the usual meaning of these initials, the name was switched to total variation nonincreasing (TVNI), but was eventually settled on the more euphonious TVD.

The following is an expansion of Harten’s proof of the lemma.

Schemes which are total variation nonincreasing (TVNI) can be characterized as follows:

We can thus define


and substituting
TV\left(u^{n+1}\right)\leq TV\left(u^{n}\right)

Consider the general expression

C_{{j-1,n}}\geq 0

D_{{j,n}}\geq 0 C_{{j-1,n}}+D_{{j,n}}\leq 1

We should observe that our ultimate goal is to sum these values from negative infinity to positive infinity; thus, we can shift the index at will. The inequalities will still hold but the specific location in space may change. It is also worth noting that the coefficients may themselves change at different points in space.

Let us consider the next spatial step, to wit

Subtracting the previous spatial step from this yields

Some rearranging yields

Taking the absolute value of both sides, we have

At this point we observe that
C_{{j-1,n}}\geq 0
 D_{{j+1,n}}\geq 0
D_{{j,n}}+C_{{j,n}}\leq 1
1-D_{{j,n}}-C_{{j,n}}\geq 0

We can thus limit the absolute values and write the expression as follows:

Taking the summation for both sides,
\sum_{j=-\infty}^{\infty}\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right|\leq \sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|\left(1-D_{{j,n}}-C_{{j,n}}\right)+\sum_{j=-\infty}^{\infty}\left|\left(u_{{j,n}}-u_{{j-1,n}}\right)\right|C_{{j-1,n}}+\sum_{j=-\infty}^{\infty}\left|\left(u_{{j+2,n}}-u_{{j+1,n}}\right)\right|D_{{j+1,n}}

Since, as we observed before, we can shift the indices (as the “centre” of the system is arbitrary with infinite boundaries) we can rewrite the above as follows
\sum_{j=-\infty}^{\infty}\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right| \leq \sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|\left(1-D_{{j,n}}-C_{{j,n}}\right)+\sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|C_{{j,n}}+\sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|D_{{j,n}}

in which case
\sum_{j=-\infty}^{\infty}\left|u_{{j+1,n+1}}-u_{{j,n+1}}\right| \leq \sum_{j=-\infty}^{\infty}\left|\left(u_{{j+1,n}}-u_{{j,n}}\right)\right|

Substituting, we have at last
TV\left(u^{n+1}\right)\leq TV\left(u^{n}\right)

Citation: Harten, A. (1997) “High Resolution Schemes for Hyperbolic Conservation Laws.” Journal of Computational Physics, Vol. 135, pp. 260-278.


Taking the Last Voyage with Newton and Pascal

He’s not widely known outside of the fields he specialised in, but Adhémar Jean Claude Barré de Saint-Venant (1797-1886, usually known in the Anglophone world as simply Saint-Venant) was one of the premier scientists, engineers and mathematicians of the nineteenth century.  His accomplishments were many and include the following:

  • saint-venant-torsionSuccessful derivation of the Navier-Stokes Equations for a viscous flow before Stokes; these equations are the basis of computational fluid dynamics and the analysis of things that fly.
  • Systematisation and development of methods in the theory of elasticity of solids, including his semi-inverse methods for torsion, important in things such as automobile crank shafts.
  • Methods for the analysis of wave mechanics in bars, which we see in many places, from musical instruments to driven foundation piles.

Saint-Venant was born into a royalist, aristocratic, traditionally Roman Catholic family at a time when it was not safe to be any of these: the French Revolution, at that point stumbling from the Reign of Terror to control of France–and most of Europe–by Napoleon Bonaparte.  It was about the latter where Saint-Venant made a statement about himself that got him into trouble with the “new” Europe.  As described in S. Timoshenko’s History of Strength of Materials:

The political events of 1814 had a great effect on Saint-Venant’s career.  In March of that year, the armies of the allies were approaching Paris and the students of the École Polytechnique were mobilized.  On March 30, 1814, they were moving their guns to the Paris fortification when Saint-Venant, who was the first sargeant of the detachment, stepped out from the ranks with the exclamation: “My conscience forbids me to fight for an usurper…” His schoolmates resented that action very much and Saint-Venant was proclaimed a deserter and never allowed to resume his study at the École Polytechnique.

Saint-Venant’s statement of conscience was at once a political and religious statement, and “progressives” of his day didn’t miss either.  The French, then and now innocent of anti-discrimination legislation or sentiments, made his life miserable. The École Polytechnique was and is France’s premier technical institute of higher learning; getting kicked out of it was the equivalent of, say, being expelled from Princeton or MIT.  He worked in the powder industry for nine years, then was admitted to the École des Ponts et Chausées, where his fellow students shunned him.  He graduated first in his class anyway and began his illustrious career in technical things both theoretical and practical.

In spite of his difficulties within France, his reputation outside of her was another matter.  When François Napoleon Moigno wrote his book on statics, he discovered the following:

He (Moigno) wanted the portion on the statics of elastic bodies to be written by an expert in the theory of elasticity, but every time he asked for the collaboration of an English or a German scientist, he was given the same answer: “You have there, close to you, the authority par excellence, M. de Saint-Venant, consult him, listen to him, follow him.” One of them, M. Ettingshausen, added: “Your Academy of Sciences makes a mistake, a great mistake when it does not open its doors to a mathematician who is so highly placed in the opinion of the most competent judges.” In conclusion Moigno observes: “Fatally belittled in France of which he is the purest mathematical glory, M. de Saint-Venant enjoys a reputation in foreign countries which we dare to call grandiose.”

The French finally broke down and admitted Saint-Venant into the Academy of Sciences in 1868.  He continued his work, much of it from his home, up until the time of his death.  When the President of the Academy announced that passing, he made the following statement:

Old age was kind to our great colleague.  He died, advanced in years, without infirmities, occupied up to the last hour with problems which were dear to him and supported in the great passage by the hopes which had supported Pascal and Newton.

Europeans of the time would not have missed the import of the last statement: Pascal and Newton were Christians, and Saint-Venant was being identified with them as one also.  It was also a statement that Saint-Venant, for all of his achievements and interests which have enriched the world, also had an eternal goal as well.

There’s no evidence that Saint-Venant was ostentatious in his faith walk; descriptions of his life show the contrary.  And–shock to today’s atheist–there’s no evidence that it ever impeded the progress of his research or his thought.  As the statistician and eugenicist Karl Pearson, no friend of Christianity, noted:

The more I studied Saint-Venant’s work, the more new directions it seemed to me to open up for original investigation of the most valuable kind. It suggested innumerable unsolved problems in atomic physics, in impact, in plasticity and in a variety of other branches of elasticity, which do not seem beyond solution, and the solution of which if obtained would be of extreme importance. I felt convinced that a study of Saint-Venant’s researches would be a most valuable directive to the several young scientists, whose recent memoirs shew their interest in elasticity as well as their mathematical capacity. Many of the problems raised by Saint-Venant’s suggestive memoirs were quite beyond my powers of analysis, and I recognised that the most useful task I could undertake, was by a careful account of the memoirs themselves to lead the more competent on to their solution.


The biggest impediment he had to face was the blowback from his stand at the École Polytechnique, and that came from his secularist colleagues.  But, when the end came, all of his colleagues knew where he stood, in this life and the next one.

I spend a lot of time on this site and others talking about sea (and sometimes air) voyages.  And I’ve spent most of my career (and all the academic part of same) in the applied sciences.  But when I take my last voyage into eternity, I want to do it in the same hope of Newton and Pascal–and Saint-Venant and Euler for that matter–namely that which comes from following Jesus Christ out of the grave and into eternal life.

And you should too.

Note: I have a more complete account of Saint Venant’s story with a list of references here.

“I Need to Find the Way”

When numerical methods are combined with fluid mechanics, you ultimately get Computational Fluid Dynamics (CFD.)  CFD has been a test of numerical methods for a long time, especially with the solution of the Navier-Stokes equations.  The non-linear nature of these equations has made mating them with CFD a challenging problem.

When you’re first introduced to numerical methods you’re told that there are ways one can predict the stability and consistency of these methods with a given set of physics.  These methods are helpful but things go wrong.  On top of that coding these solutions into a language like C++ or Fortran for parallel processing (or any other computer application) is a challenging business; the endless parade of bugs we see exploited in internet security can be also seen in modelling physics such as fluid flow, just different bugs.  Probably the most important lesson one learns in mastering such an art is that, just because it comes out of a computer, it can’t always be trusted.

That situation confronted a Turkish friend of mine who was working on a particularly sticky CFD coding problem.  One frustrating day he exclaimed, “I need to find the way,” i.e., out of his CFD coding problem.

Finding the way has been a key issue in the Middle East for a long time.  The last night before he was crucified, Jesus Christ had this exchange with his disciples:

“We do not know where you are going, Master,” said Thomas; “so how can we know the way?”

Jesus answered: “I am the Way, and the Truth, and the Life; no one ever comes to the Father except through me. (John 14:5-6 TCNT)

This, then, is a key difference between Jesus Christ and the others.  They knew (or they thought they knew) the way, or they have an idea of what the way is.  But Jesus Christ is the way, something he showed when he rose from the dead shortly thereafter.

My Turkish friend was doubly blessed; not only did he find the solution to his CFD problem, but he also found the way in Jesus Christ, and has been baptised.  You can see a video of another Turk who found the way here.

And for the rest of us…

For more information click here.

Computing Open Channel Flow Using a Pitot Tube

In our last post we discussed the basics of hydraulic jump.  We showed how the flow across the jump could be estimated using the depths of the water before and after the jump.  In this post we will show another method, using a pitot tube, which is also used in a wide variety of other applications.

The calculations in our earlier post assumed the velocity of the water to be uniform.  In a flow channel with parallel straight sides and rectangular cross-section, we assumed the flow to be

Q = uwy

where the variables are as before.  For both open and closed channels, however, the velocity will not be uniform; the velocity we used was an average velocity.  If we actually measured the velocity at various points in the flow, we would get different results for different places, as illustrated below (from Open Channel Flow Measurement):


The following presentation assumes a rectangular cross-section with parallel walls.

Pitot Tube Measurements

We discussed the concept behind the Pitot tube in Wind Tunnel Testing.  The simplest way to show this for open channel flow is to use the presentation of Vennard (1940):


Integrating the Results

Separating the static pressure from the stagnation (dynamic) pressure is fairly simple in open channel flow, although reading the results can be tricky if the flow level varies.  We thus take several measurements at several heights and compute a velocity at each point, getting a velocity profile for the flow.  But how to integrate this relationship for a total flow?  There are several methods.

One would be to take a trend line of the results (as discussed in Least Squares Analysis and Curve Fitting) and integrate it in closed form.  For this to work requires that the trend line have a very good correlation with the original data.

Another would be to take the results and graphically integrate them.  For many engineering applications, graphical methods were popular for many years because they avoided many difficult calculations which strained the limited computational power of the time.  Although the results were approximate, with CAD a higher degree of precision can be obtained than was possible before.  The tricky part of a graphical method such as this is the scaling, which must be understood to properly interpret the results.

Yet another is the use of numerical integration, generally piecewise with methods such as the Trapezoidal Rule and Simpson’s Rule.  With a properly laid out spreadsheet, this can be done with minimal effort, although attention to detail is crucial to success.

An implementation of numerical integration can be used which simplifies the calculations, and is suggested in Open Channel Flow Measurement.  It involves a little pre-planning in that the points where the data is taken need to be pre-determined (they should be in any case.)  Since the more points of data the more accurate the result (all other things equal,) we’ll use six points.  Those points are at the surface, at the bottom of the channel, and at 20% (0.2), 40% (0.4), 60% (0.6) and 80% (0.8) of the total depth.  The mean velocity can be thus computed as follows:


Boundary layer considerations would indicate that the velocity at the bottom of the channel be zero, but if it is possible to take a measurement it would be better.

Once the u_{mean} is calculated, it can be used in the equation at the beginning of the post to estimate the flow in the channel.


Vennard, J.K. (1940) Elementary Fluid Mechanics.  New York: John Wiley & Sons, Inc.


Getting the Jump on Hydraulic Jump

Hydraulic jump is one of the more interesting phenomena in open channel flow.  It is used, for example, in stilling basins, where it is necessary to reduce the fluid velocity coming out from a dam by dissipating the kinetic energy in the flow.  An example of what it looks like is below.

Figure 163.png

This is a very basic treatment of the subject.  It is based on Streeter (1966) and also Dodge and Thompson (1937), where many of the graphics come from.  To keep things simple we will assume that the jump takes place in an open channel with two vertical parallel walls.  But first we need to introduce a couple of concepts.

Froude’s Number

In James Warrington’s monograph about propulsive power, he mentions “the investigations of the elder Froude.”  The Froude Number, which is named after him, is basically the ratio of the inertia force to the gravity force, and is defined in Wind Tunnel Testing as

Fr = \frac{u^2}{g_cD}

Although it’s confusing, the Froude Number is frequently the square root of this value, and for this study will be designated by the equation

Fr = \frac{u}{\sqrt{g_cy}}

In this case u is the average velocity of the fluid in a one-dimensional channel; the fact that it’s average is very important.  The constant g_c is the gravitational constant and y is the depth of the fluid (usually water) in the channel.  This is the way we will use the Froude Number in this study.

Since we’re restricting ourselves to vertically-walled channels with rectangular cross-sections, we can make another modification to the Froude Number.  If we define the flow as

Q = u w y

where Q is the flow in volume per unit time, w is the width of the channel and y is defined as before, the Froude Number can be stated as

Fr = \frac{Q}{w{g_c}^{\frac{1}{2}}y^{\frac{3}{2}}}

Keep in mind this is only for channels with vertical walls and rectangular flow cross-sections.

Specific Energy

Let us consider a vertical cross section of a one-dimensional horizontal channel flow.  Bernoulli’s Equation states that the fluid energy/head at that cross section should be

E = y + \frac{u^2}{2 g_c}

What we have done here is to eliminate the pressure term in Bernoulli’s Equation, which does not apply for a surface.  We are left with the potential term (the depth of the fluid in the channel) and the kinetic term (due to the movement of the fluid.)  The energy E is thus in units of length, or is a head term.

Again assuming vertical, parallel walls and substituting the flow term, this equation can be stated as

E = y + \frac{Q^2}{2w^2 y^2 g_c}

The resulting curve looks like this:

Figure 165

Obviously this figure allows for non-parallel sides, but the basic theory is the same.  The depth which we call y is called d in this diagram.  At a depth of zero the specific energy is infinite, decreasing to a minimum at the point where the curve in (c) comes closest to the ordinate, and then decreases to the point where the second term on the right hand side vanishes and the specific energy equals to the depth.  We can find the location of this minimum by differentiating with respect to depth as follows:


Solving this equation for y yields a cubic equation with one real and two complex roots.  Using the real root only, the depth at which the energy is minimised is given by the equation

y_c = \frac{Q\frac{2}{3}}{w^\frac{2}{3}g_c^\frac{1}{3}}

This is referred to as the critical depth, and flow at this point is referred to as critical flow.  For depths shallower than this, the flow is referred to as rapid or supercritical.  For depths above this, the flow is referred to as subcritical or tranquil.

Some algebra will also reveal the following:

  1. At critical flow, the Froude Number is unity.
  2. At critical flow, the specific energy is \frac{3}{2}y_c .
  3. Since the Froude Number is inversely proportional to y^\frac{3}{2} , for values of y greater than critical (subcritical/tranquil flow) the Froude Number is less than unity.
  4. Conversely, for values of y less than critical (supercritical/rapid flow) the Froude Number is greater than unity.

Solving for Hydraulic Jump

Now that we’ve laid the groundwork for some basics of flow in channels, it’s time to consider hydraulic jump itself, which is illustrated below.

The first thing worth noting is that the incoming flow (from the left) must be supercritical for hydraulic jump to take place.  From that the flow on the right must be subcritical.  In the process of the jump energy is dissipated, which is a big purpose of inducing hydraulic jump in, say, a stilling basin.  But this also means that conservation of energy does not apply here.  What does apply is a) conservation of momentum and b) conservation of mass flow.

Starting with conservation of momentum, after a little algebra and using our y notation for depth (as opposed to d in the figures) we have

g_cy_1^2 + u_1^2y_1 = g_cy_2^2+u^2_2y_2

After a great deal more algebra and application of the Froude Number,

\frac{1+2Fr_{1}^2}{Fr_1^\frac{4}{3}} = \frac{1+2Fr_{2}^2}{Fr_2^\frac{4}{3}}

Conservation of mass flow requires that

Q_1 = Q_2



There is more than one way to combine the conservation of momentum and mass flow equations.  One way results in determining the relationship between the two Froude Numbers, thus

Fr_2 = \frac{2\sqrt{2}Fr_1}{\left( \sqrt{1+8Fr_1^2}-1 \right)^\frac{3}{2}}


Plot of the Froude Number before the jump (Fr1) and after (Fr2)

We can also reapply the conservation of mass flow equation and eliminate Fr_2 in this way:

\frac{y_1}{y_2} = \frac{2}{ \sqrt{1+8Fr_1^2}-1 } 

We now have a definite relationship between a) the two Froude Numbers and b) the ratio of the water depths on both sides of the jump to the first Froude Number (and thus the second.)

Determining the Flow

At this point we have the basics of hydraulic jump.  There are many ways we can apply this theory.  In this monograph we’ll look at two: using the results to estimate the flow, and using the results to estimate the energy loss.

We’ll start with presenting this graphic overview of hydraulic jump and the effects of varying the upstream Froude Number (and thus the downstream one, as we saw earlier) from Research Studies on Stilling Basins, Energy Dissapators and Associated Appertuances.


Let’s consider a worked example, taken from Streeter (1966).  A hydraulic jump takes place where the upstream depth is 5′ and the downstream depth is 31′.  The sluice is 50′ wide.  Determine the Froude Numbers on either side of the jump and the flow through the jump in cubic feet per second.

We start by rewriting the equation that relates the ratio of depths to the upstream Froude Number (1) as follows:

\frac{y_1}{y_2} - \frac{2}{ \sqrt{1+8Fr_1^2}-1 } = 0

We then note that \frac{y_1}{y_2} = \frac{5}{31} = 0.161.   We could go through more algebra and solve for Fr_1 , but it’s simpler to use goal seek and compute it by finding the zero of the above expression; it comes out to Fr_1 = 4.72.

We then substitute this result into this equation

Fr_2 = \frac{2\sqrt{2}Fr_1}{\left( \sqrt{1+8Fr_1^2}-1 \right)^\frac{3}{2}}

and compute Fr_2 = 0.306 .

Since we now know the Froude Numbers for both sides, we can compute the flow.  We defined the Froude Number earlier in terms of flow; rearranging that equation yields

Q = Fr w{g_c}^{\frac{1}{2}}y^{\frac{3}{2}}

Taking either the upstream Froude Number and depth or the downstream Froude Number and depth and substituting these, the sluice width (50′) and the gravitational constant, the result Q  = 14,986.4 \frac{ft^3}{sec} .

Computing the Energy Loss

One of the main purposes for inducing hydraulic jump is to dissipate kinetic energy in a fluid flow.  Earlier we saw that at any point in the flow there is a specific energy, or head, in the flow.  If we subtract the specific energy at the upstream (Point 1) from the downstream (Point 2) the difference between the two heads is

\Delta E = \frac{\left( y_2 - y_2 \right)^3}{4y_1y_2}

or strictly a function of the two water depths.  If we substitute the two water depths from the last example, this yields an energy/head change of 28.35′.  Although the units are of length, it is more accurate to say that this is an energy per unit weight of water, or ft-lb/lb.  The power dissipated during the hydraulic jump is thus

P = Q\gamma\Delta E

Substituting the head change we just computed with the flow we computed earlier and the unit weight of water yields P = 26,510,196 ft-lb/sec = 48,200 hp.

You are urged to set these equations up in a spreadsheet and substitute the values given to confirm the validity of your spreadsheet before using them for another case.

References (other than ones on this site)

  • Dodge, R.A., and Thompson, M.J. (1937) Fluid Mechanics. First Edition.  New York: McGraw-Hill Book Company
  • Streeter, V.L.  (1966) Fluid Mechanics.  Fourth Edition.  New York: McGraw-Hill Book Company.