Applied Numerical Mathematics NorthHolland
I. (1985) 309314
ON BOUNQARY CONDITIONS PROPAGATION *
309
FOR THE NUMERICAL
SIMULATION
OF WAVE
Ludwig WAGATHA Schiinstr. 43, 8000 Miinchen 90, Federal Republic of Germany
Introduction Absorbing boundary conditions are of considerable importance in the applied sciences. Consider for example the simulation of wave propagation on an unbounded surface like the surface of the earth, as is frequently the case in meteorological applications [3]. Either the calculation uses only very few grid points per actual square mile and thus gives only limited accuracy, or artificial boundaries are introduced to limit the computational area. At these boundaries (unphysical) reflections can occur, limiting the accuracy of the results in this case. It is necessary to use boundary conditions that absorb waves incident on the artificial boundary to eliminate these reflections. Metereology is only one example. Situations where it is necessary to eliminate reflections generated by artificial boundaries are rather frequent: examples range from seismology [2] to fluid dynamics [l]. Boundary conditions for these applications have been constructed by several authors (see [12] and the literature cited there). An absorbing boundary conditions asympto?ically valid at, high frequencies was derived by Engquist and Majda [4] for general classes of hyperbolic equations using the theory of reflection of wave front sets for solutions of differential equations (see [8,9,10]9. In this paper we are concerned with constructing a (nonlocal) solution operator for general difference schemes that allows to derive nonreflecting boundary conditions. We will apply this solution operator to the special cases of the 5point star for the approximation of the wave equation and the LaxWendroff scheme and construct the corresponding absorbing boundary conditions. In the third section numerical results are given. For the wave equations we get less reflected energy than with conditions of EngqrristMajda type [4,5,15], while for the shallow water equations good absorption properties can be achieved only at the cost of high storage and computing time requirements.
1. General consideratious In this section we want to motivate tile construction carried out in the next section. The results can be stated in a rigouros form (cf. [12]), but as this requires rather deep functional analytical tools like propagation of singularities, we merely try to give an outline of how absorbing boundary conditions can be characterized. Consider a strictly hyperbolic system u, = Au, + Bu, in y 2 0. Consider a boundary point
* These results are partly contained 01689274/85/$3.30
in the author’s dissertation
at the LudwigMaximilian+UniversitBt,
0 1985, Elsevier Science Publishers B.V. (NorthHolland)
Miinchen
L Wagatha / Numerical simulation of wave propagation
310
(x, 0), where the bicharacteristics intersect the boundary at a positive angle and a boundary condition /IU = g at this boundary point. Assume that a bicharacteristic yi points out of y >, 0 at a given time 2 = t,, yZ points inside the area (Fig. 1). Assume a wave approaches the boundary along yr. An absorbing boundary condition now can be constructed in two ways: (i) Decouple the equations in incoming/outgoing components, determine the incoming reflected waves along yZ using an appropriate projection operator and set them to zero. Then no waves will travel inside from the boundary. This is essentially the EngquistMajda construction (41, Decoupling the equations usually requires a pseudodifferential operator /? in the boundary condition pu = 0 (see [8,9,10]). Numerical implementation requires further efforts. (ii) Determine the solution at some distance from the boundary at times t c t,, apply a solution operator to calculate when the wave will hit the boundary and prescribe the correct solution at this moment as boundary value (i.e. apply a boundary condition u = u0 \ ,,_o, J3= id). Then, as is clear in the case of a difference scheme, no reflections will occur. The drawback in this case is that the boundary condition is nonlocal in time and does not seem suitable for more than two space dimensions in actual computations. We will pursue this second idea in the following section. In the case of difference equations it is possible to predict the value of the solution on the boundary using a discrete solution operator. Thus it is feasible to use this approach in numerical simulations. 2. Absorbing boundary conditions for difference schemes
We first construct a solution operator for general difference schemes. For clarity we restrict the presentation to one space variable and an equation with constant coefficients. Then we use this (nonlocal) solution operator to obtain an absorbing boundary condition. Let u/” denote the discrete solution of the differential equation u, = Au,., discretized by the stable and consistent difference scheme u;+’ = Q,u;*
+ Q,u,n + Q,uJ+,,
(2 .1) (2 .2)
where n denotes the discrete time variable t and j the discrete space variable y. Definition. Consider the difference scheme (2.2) for all j together with initial data uy at n = 0
and
uf= 0
Fig. 1,
for j < 1. Denote the solution of the difference scheme as uy.
311
L. Wagatha / Numerical simulation of wave propagation
A boundary condition uz = & f md n E,,,uy  “’ at j = 0 is called absorbing, if for the solution u of the difference ecluation (2.2) in the halfspace j > 0 satisfying the initial conditions and this boundary condition the following statement holds: ui”  ii; =0
for all j&O, n 2 0.
Theorem. Consider the difference scheme (2.2) in ( j X Ay/j E Z > with initial data uy and u; = 0 for j< 1. Then there is an operator E(n, m) E Rk*k than can be determined recursively, so that the boundary condition
u,“=
c
E(m, O)u;”
(2 .3)
1drncgn
is absorbing. The solution in the halfspace converges there to the unrestricted solution of (2.1) (if the initial data are ‘compatible ‘1. Proof. Let E( n, m) be a (matrix) solution of the difference scheme in (j X Ay/j < 1) with initial and boundary data E(0, m)=am,,Z
(I=identityinIFPkqk)
E(n, ‘I)= 2
(Z=ZeromatrixinIWk*k)
for n>O.
Then uy=
c
E(m,
j)ur,
is a solution uf the difference scheme in ( j X A y/j < 0). This solution and the solution c’: the unrestricted scheme are identical on the line 1 X Ay. As the coefficie.qcs of the schemti are constant in (j x Ay/j < 0) we have
As E is a solution of the scheme itself E(n, j+j’)=
c
E(m,
j)E(nm,
j’+l),
l
For j = 0, j’ =  1 we obtain E(n,
l)=
c
E(m, O)E(n  m, 0).
1 Grn6.z
Inserting this in (2.2) we obtain for j = 0
E(n+l,O)=Q,E(n,l)~Q,E(n,O)+Q_,(
c lfm
E(m,O)E(n

1
m,O) .
Thus we obtain E( n + 1,O) recursively for all n. The boundary condition (2.3) and the unrestricted scheme (2.2) yield identical values on the line y = 0 X Ay. The LaxRichtmyer theory gives convergence. q Remarks
(i) As the boundary condition (2.3) is nonlocal, in real numerical applications the summation in (2.3) has to be truncated to 1~ m < n 0, no < n.
L. Wugarha / Numerical
312
simulation
of wave propagation
(ii) Thisconstruction can be carried out in more general situations. The generalizations are obvious. However, due to the nonlocality of the solution operator, the additional storage required for the ‘boundary layer’ j = 1 increases rapidly in the case of more than two space variables. (iii) Wilson [13] gave a construction of the matrices E for the LaxWendroff scheme using discrete Laplace transformation. The method using the underlying idea of a solution operator seems more straightforward. (iv) Trautenberg, Gebauer, Sachs [ll] constructed a solution operator in the special case of the wave equation. We now apply the construction as an example to the LaxWendroff scheme in two space dimensions. Example Consider the strictly hyperbolic system in two space variables B E Rkgk, and the discrctization by the standard LaxWendroff scheme:
u, = Au,, + Bu,,
A,
X = At/Ax = At/Ay. Because of the second space dimension we obtain as boundary condition uTo=
c c
instead of (2.3)
E(m, 0, /  I’)u;!j”’
1<;nt
l)=E(l,O,
I)k;XA

c I
+
zE(m,O,
Il’)xE(n
 m, 0, I’)}
I’
further terms.
The recursive calculation of E can be performed finite.
by the computer.
The summation
over I’ is
3. Numerical results We use the 2dimensional wave equation u,,  u,,~  u,., = 0 as an example to demonstrate the properties of the solution operator approach. This equati’on is approximated on a 52 X 32 point grid using the usual 7point scheme. Initial conditions are given as a Gaussian wave with support inside the computational area, so that the solution is a circularly expanding wave. Using a larger computational area we obtain the corresponding ‘freespace’ solution and thus are able to calculate the energy reflected at the artificial boundaries. At these artificial boundaries the following boundary conditions are implemented alternatively: (a) The first order boundary condition u,,, 0.9u, I_,+ = 0 (cf [15]).
L. Wagatha / Numerical simulation of wave propagation
313
(b) The second order boundarv _ condition
(cf [15]).
i’yl  4, 0.55u,,I,=,=O
(c) The boundary value obtained using the solution operator for the 7point approximation of the wave equation and no = 10 (compare Section 2). (d) Same as (c) with no = 20. In (c) and (d) we use method (b) at the 8 points next to the corners. In two dimensions the solution operator boundary condition is of the form uzo=
c
c
E(m,
I I’, 0)~:;“’ ..
1
and near corners the summation over 1’ is limited with the result that the corners act as reflectors (Fig. 2). The following results for the energy of the reflected waves are obtained:
Method
w (W (c) (d)
Timesteps 10
20
30
40
50
0.6 0.4 0 0
1.5 0.9 0.2 0
2.2 1.3 0.8 0.6
2.8 1.7 1.2 1.0
2.6 1.8 1.4 1.2
Thus the solution operator approach reflects less energy than the pseudodifferential approach and the absorption properties can be adjusted easily by changing no. In an example using this method for the shallow water equations discretized by a LaxWendroff
Fig. 2. Reflected wave after 65 timesteps (Boundary condition
(d) not modified at the corners).
1, Wagatha / Numerical simulation of wave propagation
314
scheme (details of the example as in [15]) we obtained the following results: _._Method
Timesteps 10
20
30
40
50
60
WI
0.05
0.2
0.3
0.31
0.3.:
0.35
Solution operntor with no = 10
0
0.1
0.25
0.4
0.45
0.55
Modified second order EngquistMajda condition
To get better results than with the modified EngquistMajda conditions it is necessary to choose n, > 19, where the rapidly increasing computing time and storage requirements seem to make this method inefficient: For n, = 20 the computing time is about one order of magnitude higher than for n, = 10.
References 111A. Bayliss and E. Turkel, Radiation boundary conditions for wave like equations, Comm. Pure Appl. Math. 33 (1980) 707725.
PI R. Clayton and B. Engquist, Absorbing boundary conditions for acoustic and elastic wave calculations,
Buff. Seismological Sec. Amer. 67 (1977) 15291540. efficient schemes and boundary conditions for a fine mesh 131 T. Elvius and A. Sundstrom, Computationally barotropic model based on the shallow water equations, Tellus 25 (1973) 132156. 141 B. Engquist and A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comput. 31 (1977) 629651. [51 B. Engquist and A. Majda, Radiation boundary conditions for accoustic and elastic wave calculations, Comm. Pure Appl. Math. 32 (1979) 314357. 161 B, Gustafsson and H.O. Kreiss, Boundary conditions for time dependent problems with an artificial boundary, J. Comput. Physics 30 (1979) 333351. PI H.O. Kreiss, Initial boundary value problems for hyperbolic systems, Comm. Pure Apple Math. 23 (1970) 277298. PI A. Majda and S. Osher, Reflection of singularities at the boundary, Comm. Pure Appl. Math. 28 (1975) 479499. PI L. Nircnberg, Lectures on Linear Partial Differential Equations, Regional Conference Series in Math. 17 (American Mathematical Society, Providence, RI, 1973). PI M.E. Taylor, Reflection of singularities of solutions to systems of differe;ntial equations, Comm. Pure App!. Math. 28 (1975) 457478. WI E. Trautenberg, K. Gebauer and A. Sachs, Numerical simulation of wave propagation with nonreflecting boundary conditions, in: F. Mainardi, Ed., Wave Propagation in Viscoelastic Media (Pitman, London, 1982). fir hyperbolische partielle Differentialgleichungen, Thesis, WI L. Wagatha, Absorbierende Randbedingungen Munchen, 1982. 1131 J.C. Wilson, Derivation of boundary conditions for the artificial boundaries associated with the solution of certain time dependent problems by LaxWendroff schemes, Proc. Edingburgh Math. Sot. 25 (1982) l19. WI F. Breves, Introduction to Pseudodijferential and Fourier Integral Operators II (Plenum Press, New York, 1980). 1151 L. Wagatha, Approximation of Pseudodifferential Operators in Absorbing Boundary Conditions for Hyperbolic Equations, Numer. Math. 42 (1983) 5164.