Algorithms for transport equations: 1-D flow

This section begins with some very simple transport equations __ examples about which much is known, and therefore it will be quite obvious if or when and sometimes how the numerical algorithms are not getting it right.

The wave equation

The differential equation is

Ftt - c2Fxx = (t2 - c2x2)F = 0,(1)

where c is a constant, F is F(x,t), and the subscripts, 'x' and 't', denote partial derivatives. Less formally, Fx is the rate at which F varies (with x) at a fixed time t, and Ft is the rate at which F varies (with t) at a fixed point x. The use of tF = Ft(x,t), etc, is just a different way of talking about the same thing that can sometimes be quite useful. For example, some of the basic ground rules are: tf(x) = 0 , xg(t) = 0 and (usually) tFx = xFt. (The condition is that the mixed partial derivatives have to be a continuous function of x and t, and proving that for either one will do.)

From the time of Fourier (ca. 1825) much work on the wave equation has employed Fourier Analysis. Here the emphasis is on the older (ca. 1746) d'Alembert's method. To introduce it, consider first a function F(q(x,t)). One of the very most useful pieces of calculus-lore, often called the chain rule, states that tF(q(x,t)) = Fqqt, where Fq is the rate at which F varies (with q). Likewise, xF(q(x,t)) = Fqqx. Now let q = x ± ct, and assume that the second derivative of F with respect to its argument (qFq) is continuous; then (incidentally) tFx = xFt and (importantly) F is a solution of eq(1). So d'Alembert's general solution of eq(1) is

F(x,t) = F(x - ct) + G(x + ct) ,(2)

where F(q) and G(q) are arbitrary functions, except for the requirement that they be smooth enough to qualify as 'waves' (piecewise continuous will do quite nicely). By the same kind of argument (Gt = cGq, Ft = -cFq),

Ft + cFx = 0   and   Gt - cGx = 0 ,(3)

and those are examples of the very simplest of transport equations. The excuse for including one of them here is the hope that workings of the semi-Lagrangian algorithm will be clearer in a case where the answer is known.

A particularly simple case is

F(x,t) = F(x - ct,0)   with   F(x,0) = îìí 1 for x < 00 for x > 0 .(4)

For x > 0, nothing happens until t = x/c, then the front arrives, and thereafter nothing new happens. This one is a good test case for algorithms that are meant to treat more general cases.

More generally
Ft + cFx = 0   implies   F(x,t) = F(x - ct,0),(5)

no matter what the shape of F(x,0) may be, and another way to say that is to introduce a 'path' x = X(t) along which dX/dt = c. Then X(t) = X(0) + ct and

 ddt----F(X(t),t) = Ft + cFx = 0   implies   F(X(t),t) = F(X(0),0).(6)

Eq(6) is the formulation that has easy generalizations to transport equations in which the constant velocity 'c' is replaced by u(x,t) and Ft + uFx ¹ 0.

This case (u = c) isn't complicated enough to require a semi-Lagrangian algorithm, but some of the basic ideas can introduced here. The algorithms for 1-d flow will be made a bit more complicated than they need to be in order to introduce ideas that will be needed for 2-d flows and for selective grid refinement. For now, these are just notes: Note: 't' never appears explicitly in the V's, so neither V0 nor V1 knows what time it is. (That will go away presently.) In this case the displacements, x(k) - X(k) are all equal to cdt, and it is only F(k) that varies with time in the V's. The steps needed to follow the progress of Fk(t + ndt) for n = 0, 1, . . .   are: That generates Fk(t + dt) and Fk(t + 2dt), and it can be repeated indefinitely. More about the evaluations of F(X(k))0 and F(X(k))1 will follow in the next example.

Advection

A relatively simple example of Eulerian and Lagrangian formulations of a transport equation is

qt + uqx = s(x,t)   and   dX/dt = u(X,t) , dq/dt = s(X,t) ,(7)

where u(x,t) (dimensions L/T) is the velocity in a one-dimensional, unsteady flow of a liquid or gas, q(x,t) is a quantity (e.g. thermal energy) that is transported by the motion, and s(x,t) is a source (e.g. thermal radiation) that also contributes to dq/dt. At the far right dq/dt, often called the material derivative, is (d/dt)q(X,t). According to the chain rule dq/dt = qt + (dX/dt)qx, and the two formulations agree. The process is called advection, and as long as u(x,t) and s(x,t) are prescribed functions, known in advance, algorithms for the evolution of q(X(t),t) (given q(X(0),0), say) are fairly easy to formulate.

This time the databases (for vertices) are Again the fixed increments of 'x' and 't' are dx and dt, and this time the line Xk(t) that passes through x(k) (at time t) is at a position Xk(t - dt) that is defined for each 'k' by


X(t) = x = X(t - dt) +  ó õ   tt-dtu(X(t'),t')dt'.(8)

This is not a very easy equation for X(t-dt), especially if values of u(x,t) are known only at the vertices and at time steps separated by the increment dt. Several layers of approximation will be introduced for numerical algorithms that almost get eq(8) right.

First, let be assumed that u(x,t) is known everywhere _ between the vertices and at all times. Then the first layer of approximation is replacement of the integral by a quadrature formula, and the three of those that will be considered are:

x = X(t-dt) + u(X,t-dt/2)dt

x = X(t-dt) + (u(X,t-dt) + u(x,t))dt/2

x = X(t-dt) + (u(X,t-dt) + 4u(X,t-dt/2) + u(x,t))dt/6		 the midpoint rule

the trapezoid rule

Simpson's rule


Note, u(X,t') is shorthand for u(X(t'),t') in the formulas, and the only one of them that is a functional equation (one function, one variable, any number of fixed parameters) is the second one, the trapezoid rule.


A digression: The midpoint rule and Simpson's rule can be rewritten with a further approximation that converts them both to functional equations, as follows: Simpson's rule is by far the better choice, but it entails quite a lot of extra computation, and it is unlikely that it will be implemented for this project. But just in case, the basic formulas are

U(t - dt/2) = (U(t - dt) + u)/2 + (dU/dt - du/dt)dt/8
and
x = X(t-dt) + (U(t - dt) + u)dt/2 + (dU/dt - du/dt)dt2/12 ,

where dU/dt is evaluated at t-dt . Both results are exact if U(t') is a cubic function, but there is still a functional equation to be dealt with.


Approximations of X(t-dt) are generally found by a sequence of approximations (successive approximations), and two commonly used methods are successive substitution and Newton's method. For the trapezoid rule they are

 
X <--  x - (U + u)dt/2 successive substitution
X <-- X + dX , dX = (x - X - (U + u)dt/2)/(1 + Uxdt/2) Newton's method

In these formulas X is X(t-dt), u is u(x,t), U is u(X,t-dt), Ux is ux(X,t-dt), and the symbol '<--' signifies: evaluate the expression to the right using the current value of X, and assign to X the updated result. Either method can be initiated with X = x, if there is no better estimate to be found. Newton's method, though somewhat more complicated, generally converges much more rapidly than successive substitutions, and the first pass (X = x) just might be close enough for all practical purposes. From X = x (and U = u) the first updated estimate (Newton's method for the trapezoid rule) is

X = x - udt/(1 + uxdt/2) ,(9)

and this should be used to define at least one alarm: if 1 + uxdt/2 is near zero (less than 1/2, say), dt is too large. (There will be other alarms.)

For the quantity q(x,t) the situation is different because values of q are known only at the vertices and at the time intervals 'dt'. The evaluation of Q(t-dt) = q(X(t-dt),t-dt) generally requires an interpolation between vertices of q(xk,t-dt), and the same applies to s(x,t) if it is known only at vertices. An approximate evaluation of q(x,t) at x(k), from the trapezoid rule, is

q(x,t) = q(X(t-dt),t-dt) + (s(X,t-dt) + s(x,t))dt/2 .(10)

To evaluate Q(t -dt) = q(X(t-dt),t-dt) the interpolator (see Interpolation on edges) needs to know the interval, the values of q(x,t-dt) at vertices 'ww', 'w', 'e' and 'ee' that surround the interval, and one or both of the fractional displacements (e+w=1). ('ww' and 'ee' are not needed for linear interpolation.) In one dimension, with a fixed separation of vertices, the index of 'w' is a convenient identifier of the interval, and the vertices of the kth interval are { k-1 , k , k+1 , k+2 }. Though it is not needed yet, the second kind of database, of which only one is needed, is (I for intervals, will be followed later by P for patches, E for edges, F for faces, G for the global grid, and a few others.) So the scheme is that the algorithm that updates X(t) and q(t) in V0 or V1 calls for something like

ntrp1( V , kw , e , q , s , V' , k , q' , s' ),

where kw and I locate data in V = V0 or V1 for interpolations defined by 'e' that are then assigned as data in the kth element of V' = V0 or V1. If values of 's' are known everywhere, interpolations of 's' are not included unless the evaluation of s(X,t-dt) is more trouble than the interpolation.

Burgers' equation - I.

The simplest model of the development of steep fronts (shock-waves) is Burgers' equation,

ut + uux = nuxx .(11)

The function u(x,t) (dimensions L/T) is to be thought of as the velocity in a one-dimensional, unsteady fluid flow, and the constant n (dimensions L2/T) is the kinematic viscosity. (Note n is the greek letter nu, not the latin v.) Of particular interest is the case where n is zero, and ut + uux = 0 is called the inviscid Burgers' equation. The three formulations of it are

(a): ut + uux = 0   or   (b): dX/dt = u --> du/dt = 0   or   (c): ut + (u2/2)x = 0 .(12)

The Eulerian form (a) will not be used to define numerical algorithms, so consider now the Lagrangian form (b). A semi-Lagrangian algorithm (look backward in time) will follow presently, but first note that the initial value problem has a really easy exact solution that works for a while, but not forever. According to eqs(12b), the value U(0) of u(x,0) at x = X(0) does not vary with time on the (straight) line X(t) = X(0) + U(0)t. In principle, all that need be done is to construct all the lines through values of X(0) and look forward (or backward) in time to find u(x,t) (cf eq(5)).

A simple example (cf eq(4)) and a good test case is


u(x,t) = u(x - ut,0)   with   u(x,0) = îìí 1 for x < 0 1 - x for 0 < x < 1 0 for x > 1 .(13)


For this one the lines for 0 < X(0) < 1 all intersect at X(1) = 1 where U(t) has every value in [0,1] __ in other words, the solution quickly evolves to the front of eq(5), displaced one unit to the east. That will suffice for t < 1, and for t > 1 we look for a steady (ut = 0) shock wave that travels with a velocity U that is neither 1 nor 0. In the moving coordinate system the conservation law (eq(12c)) is (u - U)2x = 0, with u = 1 to the west and u = 0 to the east, and U = 1/2 does the trick. More generally, the shock velocity is the average of the the velocities at the discontinuity. This treatment is consistent with similar treatments of shockwaves in gas dynmics. The more important trick will be to find an algorithm that gets both aspects of these flows qualitatively right.

With dU/dt = 0 the functional equation for X(t-dt) is

find   X   such that   X + u(X,t-dt)dt = x,(14)

and no quadrature formula is required. In general, X and U(X) will be found by interpolations of data at t-dt, and if Newton's method is to be employed, estimates of ux(xk,t-dt) will also be needed for interpolated values of UX(X). Then the iteration is

X <-- X + dX where (1 + UXdt)dX = x - X - Udt.(15)

This needs an initial estimate of u for a first displacement X = x - udt/(1 + uxdt), and the sensible one is the discrete form of the conservation law,

u(k,t+dt) <-- u(k,t-dt) - (dt/dx)(u2(k+1,t) - u2(k-1,t))/2.(16)

This is called a leap-frog algorithm; it is suitable for capturing shockwaves, but because of its numerical dissipation it needs a correction everywhere else. Eq(15), with corrections for the numerical dissipation that comes with linear interpolation, will be implemented and results will be compared.

Burgers' equation - II.

Another model of the development of steep fronts is the nonhomogeneous Burgers' equation,

ut + uux = nuxx + s(x,t).(17)

Again, the 'viscosity' will be neglected (n = 0) in order to concentrate on the fact that the X- and u-equations are coupled now. The use of the rectangle rule gives the functional equations,

u = U + (S + s)dt/2 + O(dt2) and x = X + (U + u)dt/2 + O(dt2) ,
(18)

where X = X(t-dt), U = u(X,t-dt), S = s(X,t-dt) and O(dt2) indicates the power of dt in the error estimates.

In this particularly easy case, it makes some sense to compare eqs(18) with results obtained from the use of Simpson's rule. The Hermite cubic approximation of Q(t-dt/2) is (Q + q)/2 + (dQ/dt - dq/dt)dt/8, and the improved quadrature formulas are

u = U + (S + s)dt/2 + (S' - s')dt2/12 + O(dt4)

x = X + (U + u)dt/2 + (U' - u')dt2/12 + O(dt4)

= X + Udt + (2S + s)dt2/6 + O(dt4) ,


(19)



where Q' = dQ/dt and q' = dq/dt. Note: the neglected term in the third of eqs(19) is (S' - s')dt3/24 = O(dt4), and the improvement of the error estimates for all three of eqs(19) is a welcome gift _ Simpson's rule is one of the best bargains in all of numerical analysis. The reason for including these is that tests will be implemented to see if algorithms based upon eqs(19) are qualitatively different from the ones that use eqs(18).

Shallow water equations

The physical situation is a layer of water, flowing with a velocity that depends upon (x,z,t) for 0 < z < h(x,t). The transport equations are

ht + uhx + hux = 0 and ut + uux + ghx + px = 0,(20)

where u(x,t) is the mean velocity (averaged over z) and 'g' is the acceleration of gravity (dimensions L/T2). The prescribed source is p(x,t) = pA/r, where pA is the atmospheric pressure (at z = h) and r is the (constant) density. Eqs (20) govern the simplest of several nonlinear models of shallow water waves. It contains shockwaves (called bores), but not undular jumps (to be included later because they are somewhat like mountain waves).

This time (at last) the transport equations are almost like gasdynamics; they are coupled equations in which the material derivatives, dh/dt and du/dt, and the gradients, ux and hx, are both present. The conservation laws that follow from them are

ht + (hu)x = 0 and (hu)t + (hu2 + gh2/2)x + hpx = 0,(21)

and the simplest (leapfrog) predictors for h and u are

h(x,t+dt) = h(x,t-dt) - 2dt(hu)x(x,t) and (hu)(x,t+dt) = (hu)(x t-dt) - 2dt((hu2 + gh2/2)x + hpx)(x,t). (22)

Estimates of qx employ values of q at points w _ c _ e that are separated by the fixed distance dx, and corrections that use ww _ w _ c _ e _ ee can also be used. The estimates are

qxcdx <-- (qe - qw)/2 and qxcdx <-- (1 + a)(qe - qw)/2 - a(qee - qww)/4 , (23)

where a varies from zero to one according to which estimates are used for values of qxe and qxw. The corrected formula comes from a quartic Hermite interpolation (qw, qc, qe, qxw and qxe), and if central differences (the first expression above) are used to estimate the slopes at 'e' and 'w', a = 1/2.

Before the introduction of Lagrangian correctors, let it be noted that there is also a leapfrog predictor for X. If it is agreed that no further corrections of the trapezoid rule will be used, the approximatate solution of dX/dt = u(X,t) is a functional equation for a variable X(x,t) __ fairly near X(t-dt), and defined by

x = X(x,t) + (u(X(x,t),t-dt) + u(x,t))dt/2 , (24)

from which follows (1 + Uxdt/2)Xt = -(Ut + ut)dt/2, where U, Ux and Ut (at X and t - dt) are all to be defined by interpolation. The predictor that follows from the second of eqs(20) is

X(x,t+dt) <-- X(x,t-dt) + (UUx + gHx + Px + uux + ghx + px)dt2/(1 + Uxdt/2) . (25)

Note: the result suggests that X(t+dt) = X(t) would usually have been a better first guess than X(t+dt) = x.

The databases are

V0 = V1 = { k , t , x , X , h , u , hx , ux , px , (p) },

and the transition from time 't' to time 't+dt' goes as follows:
The crudest corrector uses successive substitutions in the equations for X, h and u, possibly in that order but not necessarily. That gives

X <-- x - (U + u)dt/2 , h <-- H - (HUx + hux)dt/2 and u <-- U - (g(Hx + hx) + Px + px)dt/2 ,(26)

and that's it, except for a few notes.

A perfect gas

This section and the next won't have much to say about new algorithms. They are more about introducing things that happen in models of atmospheric phenomena, where the numbers of transport equations become ever larger. The variables are density r(x,t), velocity u(x,t) and temperature q(x,t), and the transport equations are

rt + urx + rux = 0 , ut + uux + px/r = 0 , Cv(qt + uqx) + Rqux = r(x,t) , (27)

where R is the gas constant in p = rRq  and Cv (not a derivative) is the constant specific heat at constant volume (usually somewhere between 3R/2 and 3R). The source r(x,t) generally includes radiant heating and cooling among other things, and it will be taken to be zero in this section. (Note, q has been used instead of T for temperature, and Q will be q(X,t-dt), as in other sections.)

The cconservation laws that follow from eqs(27) are

rt + (ru)x = 0 , (ru)t + (r(u2 + Rq))x = 0 , (r(u2/2 + Cvq))t + (r(u2/2 + Cpq))x = rr . (28)

The substitution p = rRq has been made in eqs(28), and Cp ( = Cv + R ) is the specific heat at constant pressure.

(More to follow.)

Dispersive waves in shallow water

This will not be a very good model of an atmosphere, but it serves to introduce several of the ideas that will appear in better models. We begin with a relatively general (compressible) fluid: The special case where density is constant comes later in this section, and a perfect gas is reintroduced in the next section.

This is about one-dimensional models of two-dimensional flows that are defined by any number of quantities q(x,z,t) and a density r(x,z,t). Among the q's are velocity components u(x,z,t) and w(x,z,t), and mass is conserved, according to

rt + (ru)x + (rw)z = 0 or dr/dt + r(ux + wz) = 0 , (29)

where dq/dt = qt + uqx + wqz. Conservation laws and transport equations for the q's are all of the forms

(rq)t + (ruq)x + (rwq)z + f = rs and dq/dt + f/r = s , (30)

and each 'q' has its own sources f(x,z,t) and s(x,z,t). The only q's that will appear in this section are u, v and z, and their transport equations are

dz/dt = w , du/dt + px/r = 0 and dw/dt + pz/r + g = 0 ,(31)

with a pressure p(x,z,t) to be discussed later.

To build a shallow layer model for flow taking place between z = b(x,t) and z = h(x,t) introduce mass-weighted mean values and deviations of the q's as follows:


<f> = óõb(x,t)ôh(x,t) f(x,z,t)dz , m(x,t) = <r> , Q(x,t) = <rq>/m , q'(x,z,t) = q - Q. (32)


(Q will be replaced by qm later, so Qm can be used for qm(X,t-dt).) The fluid boundary conditions at z = b and z = h are

w(x,b,t) = bt + u(x,b,t)bx and w(x,h,t) = ht + u(x,h,t)hx,

and it follows with relative ease that

(mQ)t + (mUQ)x + <ru'q'>x + <f> = mS .(33)

The special case q = 1 (q' = f = s = 0) implies mt + (mU)x = 0 , and that can be used to convert eqs(33) to transport equations. All told, the conservation laws for mean values are

mt + (mU)x = 0 (mU)t + (mUU)x + <ru'u'>x + <px> = 0

(mZ)t + (mUZ)x + <ru'z'>x = mW (mW)t + (mUW)x + <ru'w'>x + <pz> + mg = 0.(34a)

and the transport equations are

dm/dt + mUx = 0 mdU/dt + (<ru'u'> + <p>)x = hxp(h) - bxp(b)

dZ/dt + <ru'z'>x/m = W m(dW/dt + g) + <ru'w'>x + p(h) = p(b),(34b)

where p(b) and p(h) are p(x,b(x,t),t) and p(x,h(x,t),t), and dQ/dt is Qt + UQx.

Eqs(34) are all exact, but they are not complete, and one-layer approximations of the 2-d flow require approximations of u', w' and p -- all functions of x, z and t. (z' = z - Z is the only easy one.) Estimates of u' and p are entirely different from one another, so q's will be put off until the next section (i.e. set u' = 0.), and some estimates of p will be provided in this section. Note:


That leaves p(b) and <p>, both functions of x and t, to be dealt with.

(More to follow.)

A one-layer model of an atmosphere