A hybrid computer method for the analysis of time dependent river pollution problems

by R. VICHNEVETSKY

Electronic Associates, Inc. and Princeton University Princeton, New Jersey

and

ALLAN W. TOMALESKY

Electronic Associates, Inc.

Princeton, New Jersey INTRODUCTION

This paper is devoted to the description of work done in the hybrid computer simulation of polluted rivers and estuaries., Our attention in this paper is restricted to the solution of the pollutant concentration equation.

The computer method used to perform the integration is essentially a continuous-space discrete-time method of lines. We have, in a previous paper,! described a con-tinuous-space-discrete-time computer method for the analysis of flows and velocities in a one-dimensional river or estuary. Hence, these two programs, which may be exercised simultaneously, must be viewed as part of the same problem, since the pollutant diffusion parameters in a river (as described in the present paper) may be derived· as explicit functions of the river geometry and water flow.

The kind of problems in partial differential equations to which river flows and pollution studies belong are as a rule computer time-consuming. It is therefore desirable to place emphasis on techniques by which truncation error-correction methods may lead to larger grid sizes in the finite differences processes of approxi-mation. Such a truncation error characterization and correction method is embodied in the present paper, which permits the truncation error induced by larger time steps in the computer simulation to be (in the first approximation) corrected for in a semi-exact fashion.

PROBLEM STATEIVIENT

A simplified analysis of a one-dimensional river in terms of the polluting species is given mathematically

43

as a partial differential equation representing the mass balance on the pollutant in one space dimension and time (see, e.g., Reference 6 for a derivation).

where:

ac = _~ k ac _ a (V • C) _ D ( c)

+ f

(1)

at ax ax ax

C =

c

^{(x, t)} pollutant concentration V = Vex, t) river velocity (ft.jsec.)

D (c) = pollutant degradation or decay function (We shall assume for simplicity of the ensuing discussion that D (c) = D· C where D is a constant.)

k = k (x, t) diffusion constant (ft2/ sec. ) f = f(x, t) pollutant source function x = river length variable (ft.)

t = time variable (sec.)

The boundary conditions associated with this problem are discussed in a later section.

COMPUTER ANALYSIS The C SDT approximation

The hybrid continuous-space-discrete-time method of approximation consists in expressing the solution along equi-distant lines parallel to the x axis in the

(x, t) plane.

Call c^j (x) the approximation of c (x, t^j) where t^j =

.f •

~t;

j = 0, 1, 2, .. '.'

Then equation (1) may be approximated by the

se-44 Spring Joint Computer Conference, 1970

t=o

o 2

4

6 8

10

t=25

o ²

⁴

^{6 8}

¹⁰

t=50

o 2

4

6 8

¹⁰

o ²

⁴

⁶

Distance - x 100C ft,

Figure 1-Typical propagation of pollutant profile

quence of ordinary differential equations (for j 1,2 ... ).

--~ = 0 - k^J+¹ ~- -

-CiH - cJ.' [ d , dci+1 d(Vi+lCi+l)

~t dx dx dx

_ DCi+l

+

P+l]

+

^{(1 _} 0)

[.!i

^{kJ dc·i' _ d(Vici)}

dx dx dx

0,

- Dei

+ p]

⁽²⁾

where 0 is a constant which must be chosen in the interval

72 <

0 ~ 1 to ensure stability in the time marching process.²

To produce a recurrence relation for the time

march-ing solution of (2), we solve that equation for (Ci+l) : d . d , . dci+¹

(1

^{dVi+1 )}

_kJ+I-C1H - Vl+1-- - -

+ - - +

^{D Ci+l}

dx dx dx O~t dx

cⁱ

(1 - ^{0) [d d}

- O~t - fi+l - - 0 - dx ki dx ci

d(ViC^{i )} ]

- -Dci+p

. dx (3)

Let the right-hand side equal:

where:

iY

=

ci

+

^{(1 - 0) •}

~t [.!!:.-

k idci _ d(Vicⁱ)

dx dx dx

- Dcⁱ +p] (4) I t is easily shown that E·i satisfies the recurrence relation :3,4

(5) For convenience, we call Ri the entire right-hand side of equation (3) :

We can see that equation (3) can be rewritten as:

d dCi+1 dCi+1

(1

^{dV )}

_ki+l_- - Vi+1- - - - -

+ - +

D Ci+l

dx dx dx 0 • ~t dx

= Ri (6) Ri is a known function of x and this equation can now be solved at each time step, together with the algebraic calculation of Ei+l as expressed by (5).

Stability problem and application of the method of decomposition

Equation (6) is of the second order in x. For constant V and k, its characteristic equation is:

k"Y2 -

V

^{"y -}

(~

_O~t

⁺ D) ⁼

^{0 (7)}

or

V

~V2

^{( 1 )}

- y = - ± - + k - - + D

2 4 () • t::.t (8)

These two values of -yare real and of opposite sign.

Thus, direct integration of equation (6) as an initial value problem of the second order will have unstable error propagation properties whi~h may impair the validity of the computer results.

The Method of Decomposition (Vichnevetsky,2.3) consists in avoiding the difficulty by transforming this second order differential equation into two first order ordinary differential equations, for which directions of stable x-integration may be chosen independently.

This is obtained as follows:

The second order operator

appearing in equation (6) is (arbitrarily) decom-posed into the product of two first order operators, LB and L F, which are intended to yield stable integrations in the backward and forward direc-tions, respectively:*

Conditions for the stable integration in these respec-tive directions are AB ~ 0 and AF ~ O.

By identification of (10) with (9), we find:

d d d

(1

^{dV )}

L = - k - . - V - - - - + - + D dx dx dx () • t::.t dx

d d dAF d d

= - k - - - - AF - - kAB -

+

AFAB

dx dx dx dx dx

or:

(11)

* The operator LF( . ) is said to be forward-stable if all solutions of the equation LF(V) = 0 are stable in the classical sense. The operator LB( . ) i~ said to be backward-stable if all solutions of the equation LB(V) = 0 are stable in the classical sense when the integration variable (-dx) is used instead of dx An operator L( . ) is said to be unstable when it is neither forward-stable, nor backward-stable.

A Hybrid Computer IVlethod 45

BACKWARD INTEGRATION

~F(JC .... )ao

BACKWARD STABLE INTEGRATION

v (x..,.

0

FORWARD. STABLE INTEGRATION

~=-tl+ ~+ (g+g+ D)

I~

Figure 2-Computing sequence block diagram and:

"1

...

_t,

•

...

dAF 1 dV

- - - A F A B = - - + - + D (12) dx () • t::.t dx

AF may be obtained by the integration of the Ricatti equation:

dAF _ A (V - AF) = _1_ dV D dx F k () • t::.t

+

^{dx + (13)}

and AB subsequently obtained by the application of equation (11).

Now, a particular solution of equation (6) is obtained by the following sequence of computer integrations:

(a): LB(y(x»)

== d~

Y - AB(Y) = Ri(x) (14)

(b): (15)

Indeed, that Ci+l satisfies (6) is easily shown:

L(Ci+l) = LB . LF(Ci+l) = L B(LF(ci+¹») ⁼ LB(y) = Ri q.e.d.

In summary, the sequence of equations solved at each time step in this problem is that of Figure 2.

46 Spring Joint Computer Conference, 1970

BOUNDARY CONDITIONS

The equation (1) is of the first order in time and second order in space. Hence, solutions are specified by one initial condition function (i.e., the initial pollutant concentration profile ^c(x, 0)) and by two spatial boundary conditions. One of these boundary conditions, c(O, t), is well defined (the pollutant concentration at the inlet of the river section under analysis), while the second boundary condition, c (Xmax, t) is not easily defined in terms of the problem formulation. However, that end-boundary condition has, mathematically, a very small influence on the solution c (x, t) for x E [0, Xmax] except for a small region which is close to xmax. Hence, one may choose this end-boundary condi-tion to take any convenient form. The one chosen in the computer implementation described hereafter is that:*

-

acl

-0 ax Xmax

(16) This condition can be seen to be automatically satis-fied by choosing as boundary conditions for the AF and yequations (equations (13) and (14), respectively):

AF(Xmax) = 0; } simple fluid transport equation:

au = _ vau

at ax ⁽¹⁸⁾

introduces a truncation error which has the effect to

"disperse" the solution u(x, t) by a diffusion-like phe-nomenon. Hence, one may look at the CSDT

approxi-* This assumption is not a limitation of the method described in this paper. Any other boundary condition C (XmaxJ t) could be chosen. This then would require the independent calculation of solutions of the homogeneous equation L(w) =

o.

as shown for instance in Reference 4.

mation of equation (1) as' an approximation process in which the diffusion coefficient k(x, t) results from that which is introduced explicitly by the computing process described in an earlier section of this paper, plus a truncation-induced diffusion effect of the CSDT ap-proximation of equation (18) followed by an experi-mental computer verification of the applicability of these theoretical results to the more general CSDT approximation of the transport-diffusion equation (1), as described earlier. The partial differential equation

(18) describes a pure fluid transport phenomenon.

The CSDT approximation of equation (18) is

The solution of (19) approximates that of a transport diffusion equation of the form:

(20) where the diffusion constant k* is a spurious diffusion coefficient, introduced strictly by the approximation process, and which depends on the parameters appear-ing in (19).

An equivalent value of k* may be estimated analyti-cally. To that effect, we express the different terms of

For the exact solution, we have the relation:

which yields:

and

au

- =

at _V

au

ax

Thus, after using these relations, (23) becomes:

(24)

au au a²u

- = -

v - +

^{(8 - ~) • V2 • I1t • -}

+ ...

(26)

at ax ax²

By identification of (26) with (20), we find the equiva-lent diffusion constant:

k* = (8 - ~) V2 • I1t (27) For 8

<

~, k* becomes negative: It is of interest to note that this corresponds exactly to the values of 8 for which the CSDT approximation is unstable.² Computer verification of the analysis

Computer verification of the preceding analysis has confirmed its first-order validity within a range of parameters which applies to river pollution problems.

In order to perform this verification, the homogeneous equation:

ac = _ V ac k a²c

at ax

+

^{c ax2 (28)}

was integrated in the manner described in an earlier section on the computer, with initial conditions c(x, 0) corresponding to a Gaussian distribution; i.e.,

A (x - Xp)2

c(x,O)

=

,..(0) • exp

-v 2U(0)2 (29)

where A is a positive constant, ^Xp the point where the peak of the distribution occurs at ^{t =} 0, and u(O) the initial standard deviation of the c(x,O) distribution.

The exact solution of (28) With (29) as initial condi-tion is (at least if the boundaries are assumed to be far enough not to have any effect upon the solution) :

( ) A [x - (xp

+

V • t)

J2

C x, t = - ( ) • exp - ( ) 2

ut 2'ut (30)

A Hybrid Computer lVlethod 47

where u(t) is the solution of:

du k

-dt u

Hence,

U(t)2 - u(0)2 = 2k • t (31) and, generally, between two instants of time tl and t2

>

tl:

Equation (30) expresses the fact that the "peak'"

moves with the flow at the velocity V, that the Gaus-sian distribution property remains preserved in time, and that the standard deviation u(t) grows as Vt.

Experimental measurement of u (t) is easily achieved, either by measuring the "peak" of the solution:

A

Cmax (t) = u(t)

.

⁽³³⁾

or by measuring the "2u" of c(x, t) at 1/

Ve

^{of the}

peak: (for c = Cmax • e-¹ⁱ² • x = Xpeak ± u). A typical input is shown in Figure 1.

For the purpose of this study, the computer program described in an earlier section was utilized for equation (28), where kc .was chosen "small" (specifically equal to .01 ft/sec^2), and the results shown on Figure 3 are,

K

500"---..---r--~~-/~'

·9

Figure 3-K* vs. (J for V = 4. and At

=

.50

-e 1·

4:8 Spring Joint Computer Conference, 1970

.5 rJ*

.4

. 3

.8 .9 1.0

e

Figure 4-'l1 * vs. ()

in effect:

k* = kmeasured ~ kc ~ kmeasured (kc small)

(Note that the range of k of interest for rIver studies is of 500 and over.)

Experimental results are shown in Figure 3. Con-sideration of this figure shows that there is a reasonably good agreement between the predicted and experimental

v~lues. of the truncation-induced diffusion constant k*.

I ndueed diffusion in the transport-diffusion equation The preceding analysis and computer verification of k* is concerned with the simplified transport equation

ae

^{= _} V

ae

at ax'

or at least with the transport diffusion equation with

"small" values of the diffusion constant k.

We are in practice interested in the equations of the

rJ

o

1 .

.8

.2

0

form:

.5

Figure 5--rJc as a function of 'l1

ae

at

1.0

where k is not "small" in the sense previously described.

The question thus arises as to what happens to k*

as a function of 0, for non-small values of k. An answer to this question was searched experimentally, by per-forming experiments' similar to that described in the preceding section, but with k as an additional free parameter. In this analysis, it was first recognized that an analysis of the dimensionless "induced" diffu-sion constant 1/* = k* /V2 • ilt could be obtained as a function of the dimensionless "explicit" diffusion con-stant 1/c = kc/V2 • ilt, thereby providing a relationship where, for a fixed value of 0, 1/* would be a function of

'J1c alone.

The argument here is that there is no reason why Buckingham's ^'If' principle cannot be applied to the behavior of computer program solutions as it applied to the field of mechanics and thermodynamics. Hence, relationships between dimensionless parameters must be absolute, save for round-off phenomena, if and where they occur.)

Experience confirmed this theory, and Figure 4 shows a dimensionless chart of the induced 1/* = k* / V2 • Ilt as a function of

°

^{and 1/} = k/V2 • ilt.

1.5

f)

Figure 6-Diffusion-corrected computer method

Weare reminded at this point that any practical computer program will entail a fixed value for (), and that 'YJ* will thereby become a function of 'YJ alone.

Figure 5 shows this more useful relationship for various values of ().

DIFFUSION-CORRECTED COMPUTER METHOD

The "diffusion.,..corrected" method simply consists in deriving the function 'YJc ('YJ) for the particular value

A Hybrid Computer l\1ethod 49

of () being used from Figure 5 and introducing the cor-rected diffusion constant kc = 'YJc V2Llt into the procedure discussed earlier.

This correcting method is applied continuously as a function of time and space, as shown in Figure 6.

CONCLUSIONS

The method of s~mulation of river pollution described in this paper has been implemented as a computer program. Experimental results have confirmed the use-fulness of the diffusion correction of Section 6 as a means of allowing larger time steps t~ be utilized. It has also been found that the general "hybrid approach"

which consists in approximating the problem in the form of ordinary differential equations offers a con-venient way to implement pure-numerical simulations.

REFERENCES

1 R VICHNEVETSKY

Computer integration of hyperbolic partial differential equations by a method of lines

Proc Fourth Australian Computer Conference Adelaide South Australia August 1969

2 R VICHNEVETSKY

Hybrid computer methods for partial differential equations Course Notes for the In-Service Seminar on Engineering Applications of Hybrid Computation Pennsylvania State University June 19-21 1969

3 R VICHNEVETSKY

A pplication of hybrid computers to the integration of partial differential equations of the first and second order

Proceedings of the IFIP Congress 68 Edinburgh Scotland August 5-10 1968

4 R VICHNEVETSKY

A new stable computing method for the serial hybrid computer integration of partial differential equations

Proceedings SJCC Conference Vol 32 Thompson Book Company Washington DC 1968

5 LANGHAAR

Dimensional analysis and theory of models J Wiley and Sons New York New York 1951 6 R BIRD W STEWART E LIGHTFOOT

Transport phenomena

J Wiley and Sons New York New York 1960

Im Dokument SPRI NG JOI NT COMPUTER CONFERENCE (Seite 51-59)