NASA-C ((-/7'6;055

NASA Contractor Report 178055

lCASE REPORT NO. 86-1

NASA-CR-178055

19860012733

leASE

ON THE BOUNDARY TREATMENT IN SPECTRAL METHODS

FOR HYPERBOLIC SYSTEMS

Claudio Canuto

AUio Quarteroni

Contract Nos. NASI-17070 and NASI-18107

January 1986

LfBRARY Cr~Y

MAR 2 6 l~~o

U'''Gl.lO, IIL~E.I\RCH CENTER

ltSRARY, N"SA

'iAMPTON. VIflGlflil/\

INSTITUTE FOR CO~~UTER APPLICATIONS IN SCIENCE AND ENGINEERING

NASA Langley Research Center, Hampton, Virginia 23665

Operated by the Universities Space Research Association

NI\SI\

National Aeronautics and

Space Administration

Langley Research Center

Hampton, Virginia 23665

ON THE BOUNDARYTREATMENTIN SPECTRALMETHODS

FOR HYPERBOLICSYSTEMS

Claudio Canuto

Istituto di Analisi Numerlca del C.N.R., Pavia, Italy

Alfio Quarteroni

Istituto di Analisi Numerica del C.N.R., Pavia, Italy

ABSTRACT

Spectral methods have been successfully applied to the simulation of slow

transients in gas transportation networks. Implicit time advancing techniques

are naturally suggested by the nature of the problem.

The aim of this paper is to clarify the correct treatment of the boundary

conditions in order to avoid any stability restriction originated by the

boundaries. The Beam and Warming and the Lerat schemes are unconditionally

linearly stable when used with a Chebyshev pseudospectral method. Engineering

accuracy for a gas transportation problem is achieved at Courant numbers up to

I00.

Research was supported by the National Aeronautics and Space Administration

under NASA Contract Nos. NAS1-17070 and NASI-18107 while the authors were in

residence at the Institute for Computer Applications in Science and

Engineering, NASA Langley Research Center, Hampton, VA 23665-5225.

I. INTRODUCTION

Spectral methods have been recently applied to the numerical simulation

of the unsteady flow of a gas in long distance transportation networks (see

[i]). The regularity properties of the solution [I0] make spectral methods

particularly effective for this class of problems. Engineering accuracy is

achieved at a lower computational cost than by using more conventional finiteorder methods.

Transients occurring in the normal operation of pipeline networks are

usually slow. This happens because variations imposed on the physical

variables at the boundaries are very slow. Moreover, their propagation toward

the interior is damped by the presence of a strong friction effect.

These reasons, as well as the need to avoid the severe stability

conditions which arise from the use of explicit methods for spectral methods,

make it highly desirable to use implicit methods in time. If a method is

considered which is unconditionally stable for the pure initial value problem,

the boundary conditions have to be incorporated into the numerical scheme in

such a way that no spurious stability restriction is introduced.

Theoretical and experimental results on the numerical treatment of

boundary conditions for finite difference approximations of hyperbolic systems

are widely available in the literature. In [6], Gottlleb, Gunzburger, and

Turkel address the issue of correctly imposing the boundary conditions in

terms of the physical variables, rather than in terms of the "characteristic"

ones. They show that the boundary conditions can be properly imposed within a

finite difference or finite element method which is explicit in time, by a

procedure of boundary corrections at the end of each time step. There is

computational evidence that their precedure works for spectral methods as well

(see, e.g., [i]).
-2-

It is the present authors" impression that implicit time advancing

methods have not received sufficient attention in [6], hence the reading of

that paper may lead to some misunderstanding about the practical

implementation of the boundary conditions in implicit methods. As a matter of

fact, the boundary corrections proposed there destroy unconditional stability,

as it will be documented below.

The purpose of this paper is to discuss in more detail the correct

implementation of the boundary conditions within an implicit time advancing

scheme when a Chebyshev collocation method is used in space. We stress that

the only unconditionally stable treatment of the boundary conditions consists

of imposing at each endpoint the prescribed physical conditions together with

certain linear combinations of the physical differential equations. The

coefficients of these combinations are those which express the incoming

characteristic variables in terms of the physical ones. Thus the equations at

the boundaries have to be incorporated into the matrix to be inverted at each

time step.

The Beam and Warming scheme and a class of Lerat-type schemes have been

chosen in the following discussion as time-marching methods for the time

discretization of a linear 2 x 2 hyperbolic system. An application to a gas

transient simulation of industrial interest is also presented.

2. THE BOUNDARY TREATMENT ON A LINEAR PROBLEM

The following simple hyperbolic system has been widely considered in the

literature as a model for more complex situations:
-3-

w + Aw = 0 - 1 < x < I, t > 0, (2.1) t x

whe re

w = (w1,w2)T; A = . (2.2)

The system is supplemented by an initial condition

w(x,0) = w0(x), -I < x < i, (2.3)

and the boundary conditions

w1(-l,t) = Wl(1,t) = 0, t > 0. (2.4)

The initial-boundary value problem (2.1) - (2.4) is well-posed (see, e.g.,

[4]); more precisely, for each t > 0 one has

f Ilw(x,t) 12dx < filw0(x) 12dx. (2.5)

-1 -1

The matrix A has two eigenvalues of opposite sign

3 1

_I = r ' _2 = - _ •

(2.6)

It can be diagonalized as

A = TAT-I (2.7)

with

A = diag(_1,_2); T = T-1 =1(:_ _Ii) . (2.8)
-4-

The system (2.1) can be written in diagonal form as

+ Az = 0 (2.9) zt x

where

z = (u,v) T = T -I w (2.10)

are the characteristic variables. The boundary conditions (2.4) become

(u + v)(-l,t) = (u + v)(+l,t) = 0, t > 0. (2.11)

We consider a spatial approximation of (2.1) based on the Chebyshev

collocation method at the points

X. = COS J_

3 --N ' J = O,.-.,N. (2.12)

For each t > 0, the approximate solution, which we still denote by w(x,t),

will be a couple of polynomials of degree N in the x variable. The

degree N will be kept fixed throughout the paper.

We denote by

D = {dij }0_i, j<_N (2.13)

the matrix of the Chebyshev pseudospectral derivative at the points (2.12)

(see, e.g., [7]). We recall that if an N-degree polynomial w is identified

with the vector _ of its values at the points (2.12), then Dw is the

vector of the values of wx at the same points.
-5-

We will now introduce time discretization methods for the previous

system. Hereafter, we will say that a method is stable for a given At if

the computed solutions wn at the times t = nat satisfy an estimate of the n

form

llwnllL2(-I,i) _< C n = 1,2,3,...

with C independent of n. Unconditional stability will mean stability for

all At > 0.

We first give numerical evidence to the fact that if the boundary

correction procedure described in [6] is applied within an implicit time

advancing method, a severe stability restriction occurs.

The differential system (2.1) is collocated at all the points (2.12)

(including the boundary points) and advanced in time by one step of the Beam

and Warming method [2], which here reduces to the classical Cranck-Nicolson

scheme :

_n+l At (~wxn+l) n At (wxn) w. + A = w. ---A j = 0,...,N. (2.14)

j "-_ j j 2 j

(Let us recall for further reference that the Beam and Warming scheme applied

to the conservation law w t + f(w) x = 0 is, before space discretization:

Aw n + --_ At (A(wn)Awn)x = _Atf(wn)x (2.15)

where Aw n = w n+l - w n and A(w) = f'(w)).

Nn+l

The interior values are retained, i.e., w_+I = w. for j = I,...,N-I,

J J

while the boundary values _+i and _n+lwN are corrected by solving the

linear systems
-6-

n+l

"0 = 0

(2.1)

, ,n+l _ _n+l ,~ ,n+l ,~ ,n+l

tWl) 0 + (_2J0 = _Wl) 0 + tw2) 0

I 3n+l

and (Wl. N = 0

(2.1)

(wl_n+l "N - ._W2).n+l N = (W)Nn+l - (_2)_+1

Numerical experiments show that this method is stable only if

At _ .74 At,

where At, - 16 is the stability limit for the Modified Euler method.

3N2

Although the Cranck-Nicolson method is A-stable for a scalar equation, the

stability condition introduced by the boundary treatment is more severe than

that of an explicit method.

The proper way to treat the boundary conditions is naturally derived

using the characteristic form (2.9) of the hyperbolic system. In this case it

is possible to represent the spatially discretlzed system as an O.D.E. system

in the form

dz

d--t+ M_ = 0, (2.18)

where the matrix M already takes into account the boundary conditions. It

is clear that if M is diagonalizable and has the spectrum in the right

complex half plane, then any A-stable time discretization method yields an

unconditionally stable scheme for (2.18). In turn, this gives rise to an

unconditionally stable scheme for the physical system (2.1) - (2.4), via the

transformation (2.10).
-7-

It is well-known that the Chebyshev collocation method at the points

(2.12) for the scalar initial boundary value problem

u t + ux = 0 -i < x < I, t > 0

u(-l,t) = 0 (2.19)

u(x,0)= u0(x)

gives stable numerical results. (Unfortunately, so far no proof of an

estimate of the type (2.5) is known for this method. A stability result

involving global x - t norms can be found in [8].) The method can be

written as an O.D.E. system in the form

du

----+ Du = 0 (2.20)

dt

where _ = (u0,...,UN_l)T and _ is the submatrix of D obtained by

deleting the last row and column. Thus the differential equation is

collocated at the interior points and at the outflow boundary point x = I.

The matrix D has distinct eigenvalues with strictly positive real parts,

hence, any A-stable time discretization method will produce stable solutions

with no stability restriction.

In analogy with the scalar case, the spatial discretization of the system

(2.9) is as follows:
-8-

( + 3 )(x ,t) = 0 j = 0,-..,N-I

ut _ ux j

I ,t)= 0 j = I,...,N

(vt - _ Vx)(Xj

(2.21)

I v(x0,t ) =-u(x0,t )

U(XN,t) = -V(XN,t ).

In order to cast this system into the form (2.18) we eliminate the unknowns

uN and v0 through the boundary conditions. Hence, we set

z = ( ,...... )T (2.22) _ u 0 ,UN_l,V 1, ,vN

and

j = 0,...,N-I

3 3

dij 0 - _ diN i = O,...,N-I

M =

I

I

1⁄2 d%0 [ 0 - _id£m £ = I,...,N

m I m

m = I,...,N

The eigenvlaues of M have been computed for increasing values of N. They

were found to be distinct and with non-negative real parts.
-9-

It follows that the Beam and Warming scheme will produce unconditionally

stable approximate solutions of (2.21). They are defined by the system

n+l 3 At n+l. n 3 At n

u.J + (Ux )j =u.j ----(Uxj 2 2 ) J = 0,...,N-1

n+l i At n+l]

v.3 - ---- 2 2 (Vx _j = v_ + _--2 1 At (Vxn)j j = I,...,N

(2.24)

n+l n+l

uN = -vN

n+l n+l

v0 = -u0 •

We now go back to the physical equations using the transformation (2.10). At

each interior point both the characteristic equations are collocated, hence we

get

wjn+l + --_ At A(wn+l)j x = wjn - -7At A (wxjn_j, j = 1,...,N-I. (2.25)

At the boundary point x0, equation (2.24.1) with j = 0 becomes

At Awn+l At Awn+l TII [wn+l + --_ x ]I + T12 [wn+l + -_ x ]2

(2.26)

n At Aw n] + [wn At n = TII [w - --_ 1 TI2 - --_ AWx]2'

where we have set T-I = {Thk}l<h,k<2. Similarly at xN Eq. (2.24.2) with

j = N gives
-10-

r21 [wn+l +--_At A< +I ]i + T22 [wn+l +-_At Awn+l] x 2

(2.27)

= T21[w • n - --2 At AWxn] i + T22 [wn - --2 At AWxn] 2"

Finally we have the physical boundary conditions

_n+l 3n+l

(Wl. 0 = (Wl. N = 0. (2.28)

We conclude that at each boundary point the equation to be added to the

prescribed boundary condition is obtained by collocating a linear combination

of the physical equations after their discretization by the Beam and Warming

scheme. The coefficients of this linear combination are entries of the

matrix T-I defined by (2.10).

A class of time discretization schemes of implicit type has been recently

proposed by Lerat for finite difference approximations of hyperbolic systems

(see [9], [3]). The interest of such methods lles in the fact that

unconditional stability (in time) is achieved by including a second derivative

term (in space) into the scheme. Unlike the Beam and Warming method, Lerat's

schemes are dissipative in the sense of Kreiss, hence spurious oscillations

are automatically damped.

We will introduce hereafter the analogs of a subclass of Lerat schemes in

the case of spatial discretizations by the Chebyshev collocation method. The

previous considerations on the treatment of the boundary conditions will guide

us in defining a scheme which is unconditionally stable in time.
-11-

Let us recall that for the conservation law

wt + f(w) x = 0 (2.29)

a class of methods of Lerat's type reads as follows ([9])

Awn + B At2 T (Am(wn)" Awn) x = -Atf (wn) x + TAt2 (Am(wn)wn)x (2.30)

where Aw n = w n+l = wn, A(w) = f'(w) and B is a negative parameter. The

method is obtained from the two-term Taylor formula for w at t = tn. Time

derivatives are replaced by space derivatives using (2.29) and w n is

replaced by w n + BAw n in the second-order term.

The same idea can be used in deriving a "Lerat method" for solving a

general O.D.E. system like (2.18). Precisely one gets

n+l At2 M 2 n+l n At2 M 2 n

z + B T z = z - AtMz n + (I + B) T z . (2.31)

The stability properties of this method are easily investigated by a normal

mode analysis. The characteristic root of the method is

_2

I - + (I+ S)

p = (2.32)

_2

i + S _--

where % = Atp and B stands for an eigenvalue of M. If B = -I_ , a

cancellation occurs and the characteristic root is

1+!

2

p - (2.33)

I -- --

2
-12-

hence IPl ! 1 iff _e% _ 0. It follows that for 8 = -1_ the Lerat

discretization (2.31) of (2.18) is unconditionally stable when M is the

matrix associated to the Chebyshev collocation spatial discretization. The

expression (2,33) of the characteristic root shows that the Lerat scheme

coincides with the Beam and Warming scheme on the model equation (2.1), i.e.,

whenever the coefficient matrix A is constant, so that it commutes with

differentiation. The two methods differ for the conservation law (2.29), as

it can be seen by comparing (2.15) and (2.30).

If B € - I_ , p has a singularity at X = ±/27_ hence the scheme is

only conditionally stable. Since the eigenvalues of the pseudospectral

Chebyshev derivative are uniformly bounded away from the origin ([5], Sect.

2), stability is guaranteed if At is chosen sufficiently large.

Finally, let us observe that IPl _ 1 for all 8 d _1_ if _ is

imaginary. Hence, all the Lerat schemes are unconditionally stable when used

with a Fourier method in space.

Remark: The previous analysis shows that unconditional stability is

achieved (at least for 8 = -1_ ) is the stabilizer term in (2.31) is built

up by the square of the matrix of the pseudospectral derivative including the

boundary conditions. Instead, one could think of a stabilizer term which

simply involves the second derivative operator (in analogy to the finite

difference case, see [3]) with no boundary condition or, say, Dirichlet

boundary conditions. The resulting linear system would be of classical

elliptic type, for which efficient solution techniques are available.

Unfortunately, such a method turns out to be unstable, j

Finally, we recall that the scheme (2.31) can be transformed into an

equivalent scheme in terms of the physical unknowns using again the

transformation (2.10).
-13-

3. EXTENSION TO NONLINEAR SYSTEMS

Assume now that the system to be solved is

w + A(w)w + g(w) = 0, (3.1) t x

where g(w) is a vector and A(w) is a 2 x 2 matrix which can be

diagonalized as

A(w) = T(w)A(w)r-l(w) (3.2)

with A(w) = diag{%,_}, %, _ E IN, %_ < 0. The treatment of the boundary

conditions described in the previous section can be applied here, provided the

linear combinations of the physical equations at the boundary points involve

as coefficients the entries of the matrix T-l(wn).

The Beam and Warming method was used to simulate a I hour transient of

methane gas in a pipe of length 150 km. and diameter .5m. The equations of

motion are (see, e.g., [i]).

I t + (PU)x = 0 (3.3) (on) t + (p + pU2)x + fpulu I = 0

where the Moody friction factor f is equal to .01. The boundary conditions

simulate the packing of the pipe (i.e., pu is constant at the outflow and

increasing at the inflow).

Table I contains the values of the computed inflow pressure for different

values of N and At. The exact value is 74.8080... The time step At has
-14-

the form At = _At,, where At, is the stability limit of the Modified Euler

method for the same problem. The CPU times (in hours) on the Honeywell 6040

at the University of Pavia are also reported. The columns on the right

contain the computed values and corresponding CPU times produced by the

Modified Euler method, using At = .9At,.

These results show that the Beam and Warming method considered in this

paper is unconditionally stable and more convenient than an explicit method of

the same order for attaining accuracies of industrial interest. Further

results will appear elsewhere.
-15-

_FE_N_S

I. V. Battarra, C. Canuto, and A. Quarteroni, Comput. Methods Appl. Mech.

Engrg., 48 (1985), 329-352.

2. R.M. Beam and R. F. Warming, J. Comput. Phys., 22 (1976), 87-110.

3. V. Daru, " Contribution fi l'etude d'ume m_thode num_rique implicite pour

r_soudre les fiquations d'Euler," Note Technique 1983-8, Office National

d'Etudes et de R_cherches A_rospatiales, Paris, 1983.

4. K.O. Friedrichs, Comm. Pure Appl. Math., 2 (1958), 333-418.

5. D. Funaro, "Some results about the spectrum of the Chebyshev differencing

operator," I.A.N.-C.N.R. Technical Report n. 452, Pavia, 1985.

6. D. Gottlieb, M. Gunzburger, and E. Turkel, SIAM J. Numer. Anal., 19

(1982), 671-682.

7. D. Gottlieb, M. Y. Hussaini, and S. A. Orszag, in "Spectral Methods for

Partial Differential Equations" (R. G. Voigt, D. Gottlieb, and M. Y.

Hussalni, Eds.), p. 1-54, SLAM, Philadelphia, 1984.

8. D. Gottlieb and E. Turkel, in "Numerical Methods in Fluid Dynamics" (F.

Brezzi, Ed.) Lecture Notes in Mathematics, 1127, p. 115-155, Springer,

New York, 1985.
16 ̧

9. A. Lerat, C. R. Acad. Sci. Paris S_r. A, 288 (1979), p. 1033-1036.

i0. M. Luskin, Math. Comp., 35 (1980), 1093-1112.
-17-

TABLE I

Comparison between Beam and Warming and Modified Euler time discretizatlons,

with the Chebyshev collocation method in space. Exact Pinflow = 74.8080...

Implicit Explicit

N

Pinflow CPU time Pinflow CPU time

8 I. 74.8084 .83 E - 2 74.8076 .58 E - 2

I0. 74.806 .88 E - 3

25. 74.797 .37 E - 3

50. 74.778 .19 E - 3

I00. 74.509 .ii E - 3

200. 74.103 .87 E - 4

300. 73.071 .57 E - 4

16 i. 74.8080 .38 E - I 74.8080 .12 E - 1

10. 74.806 .39 E - 2

25. 74.798 .16 E - 2

50. 74.752 .86 E - 3

I00. 74.651 .51 E - 3

200. 74.274 .28 E - 3

300. 73.458 .21 E - 3

32 I. 74.8080 .79 EO 74.8080 .98 E -I

i0. 74.8080 .83 E - 1

25. 74.807 .33 E - I

50. 74.805 .17 E - 1

i00. 74.798 .89 E - 2

200. 74.752 .45 E - 2

500. 74.460 .19 E - 2

750. 74.360 .13 E - 2
1. Report No. NASA CR-178055 12. Government Accession No. 3. Recipient's Catalog No.

ICASE Report No. 86-1

4. Title and Subtitle 5. Report Date

On the Boundary Treatment in Spectral Methods for January 1986

Hyperbolic Systems 6. Performing Organization Code

7. Author{s) 8. Performing Organization Report No.

Claudio Canuto and AlfioQuarteroni 86-1

9. Performing Organization Name and Address 10. Work Unit No.

Institute for Computer Applicationsin Science & Engrg.

Mail Stop 1320, NASA Langley Research Center 11. Contract or Grant No.

Hampton, VA 23665-5225 NASl-17070, NASl-18l07

13. Type of Report and Period Covered 12. Sponsoring Agency Name and Address

National Aeronautics and Space Administration Contractor Report

Washington, D.C. 20546 14. Sponsoring Agency Code

505-31-83-01

15. Supplementary Notes Submitted to: J. Comput. Phys. Langley Technical Monitor:

J. C. South, Jr.

Final Report

16. Abstract

Spectral methods have been successfully applied to the simulation of slow

transients in gas transportation networks. Implicit time advancing techniques

are naturally suggested by the nature of the problem.

The aim of this paper is to clarify the correct treatment of the boundary

conditions in order to avoid any stability restriction originated by the boundarie!

The Beam and Warming and the Lerat schemes are unconditionally linearly stable

when used with a Chebyshev pseudospectral method. Engineering accuracy for a gas

transportation problem is achieved at Courant numbers up to 100.

17. Key Words (Suggested by Authors{s)) 18. Distribution Statement

Chebyshev methods 64 - Numerical Analysis

hyperbolic systems

boundary conditions

Unclassified-Unlimited

19. Security Classif.{of this report) 120. Security Classif.{of this page) 21. No. of Pages 122. Price

Unclassified Unclassified 19 A02

For sale by the National Technical Information Service, Springfield, Virginia 22161

NASA-Langley, 1986