Iterative radiation boundary conditions for finite element solutions of the scalar wave equation

Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

www.elsevier.com/locate/cma

Iterative radiation boundary conditions for ®nite element solutions of the scalar wave equation Andrzej J. Safjan Department of Engineering Mechanics, University of Nebraska±Lincoln, W307 Nebraska Hall, PO Box 880526, Lincoln, NE 68588-0526, USA Received 14 April 2000

Abstract New families of radiation boundary operators are proposed for ®nite element simulations of wave propagation problems characterized by the scalar wave equation. These operators have the form Bu…j‡1† ˆ Cu…j† ; j ˆ 1; 2; . . . ; J ; where B is the Dirichlet, Neumann, or the Robin operator, and C is an appropriately chosen local operator. Since B is a standard operator compatible with the wave equation, and C acts on the known iterate u…j† ; the radiation operators can be easily incorporated into a typical ®nite element code. The proposed approaches are shown to produce very good results for the test cases considered. Ó 2001 Elsevier Science B.V. All rights reserved. AMS classi®cation: 35L05; 65M60; 65M99 Keywords: Radiation boundary conditions; Wave equation; Boundary-value problems; Finite elements

1. Introduction Numerical solution of wave propagation problems posed on unbounded domains often requires the restriction of the unbounded domain via the introduction of an arti®cial boundary and the imposition of suitable radiation boundary conditions there. The focus of this work is on developing families of highly accurate radiation boundary conditions for wave propagation problems characterized by the scalar wave equation. In particular, we are interested in developing radiation boundary conditions which: · substantially reduce the (nonphysical) re¯ections from the arti®cial boundaries; · are local; · are well-posed boundary conditions for the governing wave equation; · allow for the possibility of adjusting the aspect ratio of the computational domain; · are computationally inexpensive; and, · can be ®t into the framework of a typical ®nite element code. Various approaches have been proposed for treating computational boundaries (cf., e.g., [1±11,14,15], 1 and the literature cited therein). Usually, a sequence of boundary operators fBn gnˆ1 is derived, which approximate an operator Bexact which provides no re¯ections at the boundary of computational domain. Operators Bn may be nonlocal (see e.g., [10]), or may involve high derivatives of the solution (see e.g., [1,4,5,9]), and, therefore, may be dicult to implement in a standard ®nite element code. A di€erent formulation was proposed in [2,3], which avoids high derivatives of the solution by means of introducing

E-mail address: [email protected] (A.J. Safjan). 0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 1 ) 0 0 2 3 8 - 9

6374

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

auxiliary boundary functions and solving an auxiliary initial boundary-value problem on computational boundary Cc which is coupled with the initial boundary-value problem in computational domain X. Auxiliary boundary functions are also used in [8,11]. Finally, the radiation boundary conditions introduced in [1,6,10] seem to be limited to spherical (or circular) boundaries, and, hence, they do not allow for the possibility of adjusting the aspect ratio of the computational domain. On the other hand, spatially local boundary conditions which apply to half-space or periodic problems can be extended to rectangular computational domains with an aid of appropriate compatibility conditions at the corners, see [2,3]. In this paper, new families of iterative radiation operators are introduced for time-dependent problems and rectangular computational domains. These operators have the form   Bu…j† jCc ˆ Bu…j 1† ‡ a… on A†u…j 1† Cc ; j ˆ 1; 2; . . . ; J ; a 2 R; …1† where Cc is the computational boundary, B is the Dirichlet, Neumann or the Robin operator, on is the normal derivative, … on

A†ujCc ˆ 0

…2†

is an approximate radiation condition on Cc , Au is evaluated by solving an auxiliary initial boundary-value problem on computational boundary Cc ; u…j 1† is a known iterate, and j is the iteration counter. The main advantage of this formulation is that the auxiliary problem on Cc is decoupled from the problem in X which makes the operators easy for implementation in a standard ®nite element code. The iterative boundary operators are characterized by excellent convergence rates and produce very good results for the test cases considered. 2. A model problem Consider the wave equation in Cartesian coordinates u;tt

r2 u ˆ 0

…3†

for 0 < x < a, t > 0, y 2 R. Here u is a real-valued function of x ˆ …x; y† and t, r2 is the Laplace differential operator, r2 u ˆ u;xx ‡ u;yy ;

…4†

and the units have been normalized so that the wave propagation speed is equal to unity. The goal is to construct boundary conditions that cause waves propagating from the interior of X ˆ fx : 0 < x < a; y 2 Rg

…5†

to pass through the boundary Cc ˆ fx : x ˆ 0; y 2 Rg;

…6†

without being re¯ected (see Fig. 1). More speci®cally, we seek a boundary operator Bn , such that if w is the solution of the boundary-value problem (7) in the half-space X1 ˆ fx : x < a; y 2 Rg; 8 w;tt r2 w ˆ 0; x 2 X1 ; t 2 …0; T Š; > > > < wj ˆ g…y; t†; xˆa p …7† lim r…w;r ‡ w;t † ˆ 0; > > r!1 > : w ˆ 0 for t 6 0;

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

6375

y

θ

Γc

0

Ω g(y,t)

x

a

Γu Fig. 1. Two-dimensional Cartesian computational domain.

and u is the solution of problem (8) in the domain truncated to the strip X, 8 u r2 u ˆ 0; x 2 X; t 2 …0; T Š; > > < ;tt ujxˆa ˆ g…y; t†; > > Bn ujxˆ0 ˆ 0; : u ˆ 0 for t 6 0; then the di€erence u

…8†

w; measured in an appropriate norm, k
kX;T , is less than a prescribed tolerance 

ku

wkX;T < kukX;T : …9† p Here r  x2 ‡ y 2 and …
†;r denotes di€erentiation in the radial direction from x ˆ 0: Our notion of the solution is classical; we also assume that R the source term g has compact support in Cu ˆ fx : x ˆ a; y 2 Rg, and satis®es the following conditions: R jgj2 dy < 1; and g ˆ 0 for t 6 0. b …x; x; n† be the Fourier transform of w…x; y; t† with respect to y and t: Let w Z Z 1 b …x; x; n† ˆ p
2 e i…nt‡xy † w…x; y; t† dy dt; …10† w R R 2p where x and n are the Fourier variables and i is the imaginary unit. The inverse Fourier transform is given by Z Z 1 b …x; x; n† dx dn; w…x; y; t† ˆ p
2 ei…nt‡xy † w …11† R R 2p and the following property in conjunction with di€erentiation is essential: oa‡b w a b b: ˆ …ix† …in† w oy a otb

…12†

The transformed wave equation (13) b d2 w ‡ dx2


b ˆ0 n2 ‡ x 2 w

…13†

6376

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

admits solutions of the form b …x; x; n† ˆ w b ‡w b ‡; w

…14†

b ‡ are, respectively, the incoming and the outgoing waves, b and w where w p  2 2 b ˆ C1 …x; n†e i n x x ; w b ‡ ˆ C2 …x; n†ei w

p  2 n

x2 x

;

…15† …16†

b and w b ‡ satisfy the following oneand C1 and C2 are some functions of x and n. It can be veri®ed that w way wave equations: s b dw x2 b ˆ 0; w ‡ in 1 …17† dx n2 and b‡ dw dx

s x2 b ‡ ˆ 0; in 1 w n2

…18†

respectively. To ®nd the exact nonre¯ecting boundary condition on Cc , we restrict (18) to x ˆ 0, and get [4]: s b dw x2 b …0; x; n† ˆ 0; w …19† in 1 dx xˆ0 n2 which, after taking the inverse Fourier transform, reads ow Rw ˆ 0: ox xˆ0 Here, 1 Rw ˆ p
2 2p

Z Z R

R

e

i…nt‡xy †

s x2 b …0; x; n† dx dn: w in 1 n2

…20†

…21†

Boundary condition (20) will absorb waves propagating at any angle in the negative x direction. Unfortunately, it is nonlocal in space and time which may be undesirable from the computational point of view. 3. Iterative boundary operators Iterative boundary operators are constructed as follows. First, the exact boundary condition on Cc ; Eq. (20), is replaced with   o …22† A u ˆ 0; ox xˆ0 where A is an approximant of R which possesses certain computational advantages over R. Thus, the resulting boundary-value problem reads: 8 u;tt r2 u ˆ 0; x 2 X; t 2 …0; T Š; > > < uj ˆ g…y; t†; xˆa
…23† o > A u xˆ0 ˆ 0; > : ox u ˆ 0 for t 6 0:

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

Next, the approximate condition (22) is reformulated so that it is enforced iteratively,     o Bu…j† jxˆ0 ˆ Bu…j 1† ‡ a A u…j 1† ; a 2 R; j ˆ 1; 2; . . . ; J ; ox xˆ0

6377

…24†

where B is some preconditioning operator, a is a free parameter, and j denotes the iteration count. As a result, the original unbounded domain problem (7) is replaced with the following sequence of truncated domain problems: 8 …j† > u;tt r2 u…j† ˆ 0; x 2 X; t 2 …0; T Š; > > > > u…j† xˆa ˆ g…y; t†; > > >  < 0; j ˆ 1; …j†
…j 1†  …25† Bu jxˆ0 ˆ  …j 1† o ; j > 1; > ‡ a A u Bu > ox xˆ0 > > > > u…j† ˆ 0 for t 6 0; > > : j ˆ 1; 2; . . . ; J : The advantage of such a setting is twofold: (i) it allows for selection of B which is compatible with the governing wave equation, and, (ii) it circumvents the diculty arising from nonlocality of R, since the operator A (which may be a local as well as a nonlocal approximant of R) acts on the known iterate u…j 1† . As a consequence, the iterative boundary operators should be easy to incorporate into a standard ®nite element code. The trade-off is that the truncated domain problem (25) has to be solved several times, until

u u…J † < kuk : …26† X;T X;T Here u and u…J † are the solutions of (23) and (25), respectively, k
kX;T is a suitable norm, and  is a prede®ned convergence tolerance. However, it will be shown in Section 6.2 that, assuming that (25) is ®rst discretized in time by using appropriate ®nite difference schemes, iterations j ˆ 2; 3; . . . ; in Eq. (25) can be practically limited to the arti®cial boundary Cc and, consequently, that their cost is insigni®cant as compared to the cost of the ®rst iteration, j ˆ 1; which is the cost of solution of a single bounded domain problem. The choice of the left-hand side operators B is crucial from both the theoretical and the computational point of view. Here, we limit our attention to classical operators associated with the Laplacian (4), that is, to BN , the Neumann operator, def

BN u ˆ u;n jCc ;

…27†

and BD , the Dirichlet operator def

BD u ˆ ujCc ;

…28†

where u;n denotes the normal derivative of u at Cc . In addition, we consider the Robin operator def

BR u ˆ …u;t ‡ u;n †jCc ; which is the exact nonre¯ecting boundary operator for the one-dimensional wave equation 8 u;tt u;xx ˆ 0; x 2 …0; a†; t 2 …0; T Š; > > < ujxˆa ˆ g…t†; u ‡ u;n †jxˆ0 ˆ 0; … > > : ;t u ˆ 0 for t 6 0:

…29†

…30†

Here, (30) can be viewed as the one-dimensional counterpart of (8). We now proceed to discuss the right-hand side operator A. As already alluded to, A can be a local as well as a nonlocal approximant of R, it may be even R itself. However, having an eye on more general

6378

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

computational boundaries, we restrict our attention to operators A which can be localized both in space and in time. Accordingly, we select A in the following form:  ! M X z z d ˆ in 1 ‡ b Au bm u; …31† 1 ‡ cm z 1 cm z mˆ1 x …32† z ; n p which corresponds to approximation of 1 z2 1 by a rational function r…z† having 2M distinct real roots cm in the denominator p 1 z2

1  r…z† ˆ

P2M …z† : Q2M …z†

…33†

Here P2M and Q2M are polynomials of degree 2M, and bm ; cm 2 R, m ˆ 1; 2; . . . ; M; are unknown coef®cients which are evaluated from the following order condition p    4M‡2 r…z† ˆ O jzj : …34† 1 z2 1 In addition, for the Neumann-based iterative operators (that is the ones which have the Neumann operator BN at the left-hand side), coecients bm and cm need to satisfy an extra compatibility condition. More speci®cally, for B ˆ BN we choose P2M and Q2M in the form
P2M …z† ‡ Q2M …z† ˆ z2 1 P2M 2 …z†; …35† which guarantees that there exists an operator f A such that Au ˆ f Aox u:

…36†

Indeed, P2M …z† ‡ Q2M … z† P2M 2 …z† 2 d ˆ in…1 ‡ r…z††b b z Au u ˆ in u ˆ in Q2M …z† Q2M …z†


1 b uˆ

1 P2M 2 …z† d2 b u; in Q2M …z† dx2

…37†

from which (36) follows. Note that the wave equation (13) is used in the last step of (37). The signi®cance of constraint (35) will be discussed in the analysis of the re¯ection coecient in Section 4. Finally, combining (24) and (31), we get c …j† ˆ Cu c …j Bu with c …j Cu

1†

1†

;

…38†

! M X d 2b m c …j 1† ‡ a b u …j 1† a 1 ‡ ˆ Bu inb u …j dx cm mˆ1   M  X 2bm 1 1 1 ‡a ‡ inb u …j 2 1 ‡ c 1 c c z z m m m mˆ1

1†

1†

:

Equivalently, we take the inverse Fourier transform of (38) and (39), and get !! ! M M X X o o 2b 2b o m m 1‡ /…j 1† ; Bu…j† ˆ Bu…j 1† ‡ a u…j 1† ‡ a ox ot cm cm ot m mˆ1 mˆ1

…39†

…40†

where /…jm

1†



1 …j p 2 m

1†

‡ qm…j

1†




…41†

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

and pm…j

and qm…j 1† ; m ˆ 1; 2; . . . ; M; satisfy the following equations on Cc
ot ‡ cm oy pm…j 1† ˆ ot u…j 1† ;

6379

1†

…42†

and ot


cm oy qm…j

1†

ˆ ot u…j

1†

:

…43†

Alternatively, the two above equations can be combined into a single one-dimensional wave equation
ott c2m oyy /…jm 1† ˆ ott u…j 1† ; m ˆ 1; 2; . . . ; M: …44† To proceed further, we need to establish a criterion for selecting parameter a (cf. (24) and (39)). To this end, we observe that at z ˆ 0; all the incoming waves b u …j 1† , that is the waves satisfying p db u ˆ in 1 z2 b …45† u ; dx xˆ0 xˆ0 b have to be annihilated by operator C, c …j 1† ˆ 0: Cu

…46†

zˆ0

Introducing (45) into (46) and setting B ˆ BN , we get  i h  p p in 1 z2 a ˆ 0; in 1 z2 ‡ a

…47†

zˆ0

whence 1 a ˆ aN ˆ : 2

…48†

Similarly, for B ˆ BR and B ˆ BD , we get a ˆ aR ˆ 1;

…49†

1 a ˆ aD ˆ ; 2

…50†

and

respectively. We now ®nd coecients bm and cm , Eq. (39),pfor B ˆ BR and B ˆ BD , for M ˆ 1; 2, and 3. Toward this  end, we record the power series expansion of 1 z2 1 in the vicinity of z ˆ 0: p 1 z2

1ˆ

1 2 z 2

1 4 z 8

1 6 z 16

5 8 z 128

7 10 z 256


21 12 z ‡ O z14 ; 1024

and we compare it with the power series expansion of r…z†. Consider ®rst the case M ˆ 1. Comparing (51) with the power series expansion of r…z†,  

z z 2b1 c31 z4 ‡ O z6 ; r…z† ˆ b1 ˆ … 2b1 c1 †z2 ‡ 1 ‡ c1 z 1 c1 z

…51†

…52†

and equating the coecients at like powers of z, we get the following system of algebraic equations: 8 1 > > ; < 2b1 c1 ˆ 2 …53† > 1 > : 2b1 c31 ˆ ; 8

6380

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

which is solved by 1 c 1 ˆ b1 ˆ : 2

…54†

We observe that r…z† de®ned by (52) and (54) reads r…z† ˆ 2

z2 ; 4 ‡ z2

…55†

p which is the 4th order (diagonal) Pade approximant of 1 z2 1. Proceeding similarly, we get the following coecients bm and cm for the case M ˆ 2:     z z z z r…z† ˆ b1 ‡ b2 ; 1 ‡ c1 z 1 c1 z 1 ‡ c2 z 1 c2 z where 8 p > 1 5‡2 5 > > ; < b1 ˆ p 5 5 ‡ 1p > 1 2‡ 5 > > ; : b2 ˆ p p 5 5 1

1 c1 ˆ p ; 5‡1 1 c2 ˆ p ; 5 1

…57†

This time, r…z†; Eqs. (56) and (57), is the 8th order Pade approximant of r…z† ˆ 4z2

16

and

8 17 ‡ 64q 12q2 > >  ; p ˆ b > 1 > > 14 4 ‡ 16q 3q2 > > > < 4q 1 b2 ˆ p ; 14 q > > > > > > 39 68q ‡ 12q2 > > : b3 ˆ p2 ; 14 10 17q ‡ 3q r 28 qˆ2‡ sin 3

1: …58†

1 arctan 3



z 1

c2 z

 ‡ b3

z 1 ‡ c3 z

1 c1 ˆ p ; 2 4 ‡ 16q 3q2 1 c2 ˆ p ; 2 q



z 1

c3 z

;

…59†

…60†

1 c3 ˆ p ; 2 10 17q ‡ 3q2 r! 1 27

! 1 p : 6

This gives the following 12th order Pade approximant of r…z† ˆ 2z2

p 1 z2

2 ‡ z2 : 12z2 ‡ z4

Finally, for M ˆ 3, we have    z z z r…z† ˆ b1 ‡ b2 1 ‡ c1 z 1 c1 z 1 ‡ c2 z with

…56†

16 16z2 ‡ 3z4 : 64 ‡ 80z2 24z4 ‡ z6

…61† p 1 z2

1: …62†

We now proceed to derive coecients bm and cm , Eq. (39), for the Neumann-based operators, B ˆ BN . The procedure is the same as that for B ˆ BR except that the highest possible order condition (that is, the p algebraic equation resulting from comparing the 4M-th term of the power series expansions of 1 z2 1 with that of r…z†) is replaced with constraint equation (35), or, equivalently, with …r…z† ‡ 1†jzˆ1 ˆ 0:

…63†

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

6381

Thus, for M ˆ 1, we have   z z r…z† ˆ b1 ; 1 ‡ c1 z 1 c1 z

…64†

with the constraint …r…z† ‡ 1†jzˆ1 ˆ

2b1 c1 1 ‡ c21 ˆ 0; …1 ‡ c1 †… 1 ‡ c1 †

…65†

and the order equation 2b1 c1 ˆ

1 : 2

…66†

Eqs. (65) and (66) yield 1 b1 ˆ p ; 2 2

1 c1 ˆ p ; 2

…67†

and, together with (64), de®ne the following 2nd order approximation of r…z† ˆ

p 1 z2

1:

z2 : 2 ‡ z2

…68†

It should be noted that (68) is di€erent from Pade approximant (55). Similarly, for M ˆ 2, we have     z z z z r…z† ˆ b1 ‡ b2 ; 1 ‡ c1 z 1 c1 z 1 ‡ c2 z 1 c2 z with

p 8 1 2‡ 2 > > q  ˆ b > 1 p
; > > 8 < 2 2 p
> 2 1 2 > > b ˆ q > p
; > : 2 8 2‡ 2

c1 ˆ

1 2

c2 ˆ

1 2

r   p 2 2 ; r   p 2‡ 2 ;

which de®nes the following 6th order approximation of r…z† ˆ z2

…70†

p 1 z2

1:

4 ‡ 3z2 : 8 8z2 ‡ z4

…71†

Again, we observe that (71) is di€erent from (58). Finally, for M ˆ 3, we have      z z z z z r…z† ˆ b1 ‡ b2 ‡ b3 1 ‡ c1 z 1 c1 z 1 ‡ c2 z 1 c2 z 1 ‡ c3 z with

…69†

p 8 p 1  2 > > 2 ‡ 3 p b1 ˆ ; > > > 12 3 1 > > < 1 p 2; b2 ˆ > 12 > p > > > 1 2 > > p p
;
: b3 ˆ 12 3‡1 2‡ 3

 1 pp 2 3 1 ; 4 1 p c2 ˆ 2; 2  1 pp 2 3‡1 ; c3 ˆ 4



z 1

c3 z

;

…72†

c1 ˆ

…73†

6382

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

giving the following 10th order approximation of r…z† ˆ z2

p 1 z2

1:

16 20z2 ‡ 5z4 : 32 ‡ 48z2 18z4 ‡ z6

…74†

4. Re¯ection coecient analysis Since the boundary operators introduced in Section 3 are not exact, then, necessarily, they introduce some nonphysical re¯ections from the arti®cial boundaries. To quantify those re¯ections, we analyze the amplitude of the re¯ected wave for the case of an arbitrary plane wave u‡ u‡ ˆ f …t ‡ xc

ys†;

…75†

travelling in region X ˆ f…x; y† : x 2 R‡ ; y 2 Rg in the negative x direction, and striking the arti®cial boundary x ˆ 0 at the angle h. Here c  cos h;

s  sin h;

…76†

h is the angle between the normal to the front of the wave and the normal to the boundary, and f is some function. Postulating the existence of the re¯ected wave u in the form u ˆ R…j† f …t

xc

ys†;

…77†

we compute R…j† from the relation Bu…j† ˆ Cu…j

1†

;

…78†

which holds on computational boundary x ˆ 0. Here, u u…j† ˆ u‡ ‡ u ˆ f …t ‡ xc

ys† ‡ R…j† f …t

xc

…j†

is the sum of the incident and the re¯ected waves

ys†;

…79†

and R…j† is the amplitude of the re¯ected wave (the re¯ection coecient). We also assume that the previous stage solution, u…j 1† , is already known, and is of the form u…j

1†

ˆ f …t ‡ xc

ys† ‡ R…j

1†

f …t

xc

ys†:

…80†

To illustrate, we derive the re¯ection coecient for the Neumann-based operator de®ned by (71), which reads     d …j† d …j 1† 1 d 4 ‡ 3z2 2 b b ‡ …81† u ˆ u 1‡z in b u …j 1† : dx dx 2 dx 8 8z2 ‡ z4 First, we recall that z  x=n, which implies that 1 ‡ z2

4 ‡ 3z2 2n4 ˆ 4 8 8z2 ‡ z4 8n4

3x2 n2 ‡ x4 : 8x2 n2 ‡ x4

…82†

Next, we premultiply both sides of (81) by the denominator of the right-hand side of (82), and get d 8n4 dx


…j† 8x2 n2 ‡ x4 b u ˆ

1 8n4 2

8x2 n2 ‡ x4


d …j b u dx

1†

which, after taking the inverse Fourier transform, reads 
…j†
o 1 o 8otttt 8ottyy ‡ oyyyy u ˆ 8otttt 8ottyy ‡ oyyyy ox 2 ox

2 2n4

o 2otttt 2 ot


…j 3x2 n2 ‡ x4 inb u

3ottyy ‡ oyyyy




1†

u…j

;

1†

…83†

: …84†

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

Finally, we introduce (79) and (80) into (84), and get 

2 
1 …5† 8 8 1 c2 ‡ 1 c2 f …t ys† 1 R… j † c 2 

2 
1 1 R…j 1† c ‡ 2 2 3 1 ‡ 8 8 1 c2 ‡ 1 c2 2

2

c




‡ 1

c


2 2

6383



1‡R

…j 1†




 ˆ 0: …85†

Since (85) holds for an arbitrary y, the term in the square bracket necessarily vanishes, yielding R…j† ˆ

1 …1 c†4 1 ‡ R…j 2 1 ‡ 6c2 ‡ c4

1†




:

…86†

Note that 1 6 1 ‡ 6c2 ‡ c4 6 8 for 0 6 c 6 1 which means that R…j† is a well-de®ned function of c. Proceeding along the same lines, we get the re¯ection coecient for the Neumann-based operator de®ned by (68),     d …j† d …j 1† 1 d z2 b b u ˆ u ‡ …87† 1‡ in b u …j 1† ; dx dx 2 dx 2 ‡ z2 which reads R…j† ˆ

2

1 …1 c† 1 ‡ R…j 2 1 ‡ c2

1†




:

…88†

Similarly, the re¯ection coecient for the Neumann-based operator de®ned by (74),     d …j† d …j 1† 1 d 16 20z2 ‡ 5z4 b b u ˆ u ‡ 1 ‡ z2 in b u …j 1† ; dx dx 2 dx 32 ‡ 48z2 18z4 ‡ z6

…89†

reads R…j† ˆ

6

1 …1 c† 1 ‡ R… j 2 2 …1 ‡ c †…1 ‡ 14c2 ‡ c4 †

1†




:

…90†

Now, let us investigate the signi®cance of the constraint equation (35) for the Neumann-based operators. To this end, let us consider the following operator:     d …j† d …j 1† 1 d 2 ‡ z2 2 b b ‡ …91† u ˆ u 1 ‡ 4z in b u …j 1† : dx dx 2 dx 16 12z2 ‡ z4 p This operator is constructed in an analogous way to (81) except that 1 z2 1 is approximated by r…z† de®ned by (58) instead of that de®ned by (71). Consequently, operator (91) does not satisfy constraint equation (35). As a result, it generates the re¯ection coef®cient which is unbounded at c ˆ 0, and which takes large values at c close to zero: R…j† ˆ

5

1 …1 c† 1 ‡ R…j 2 c…5 ‡ 10c2 ‡ c4 †

1†




:

…92†

We now list the re¯ection coecients for the Robin-based operators. For the operators de®ned by (55), (58), and (62), which are, respectively,   d z2 ‡ in b u …j† ˆ 2 inb u …j 1† ; …93† dx 4 ‡ z2 

 d ‡ in b u …j† ˆ 4z2 dx 16

2 ‡ z2 inb u …j 12z2 ‡ z4

1†

;

…94†

6384

and

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398



 d 16 16z2 ‡ 3z4 ‡ in b u …j† ˆ 2z2 inb u …j dx 64 ‡ 80z2 24z4 ‡ z6

1†

;

…95†

we have, respectively, 3

R…j† ˆ

2 …1 3 ‡ c2

R…j† ˆ

4… 1 ‡ c 2 † …1 5 ‡ 10c2 ‡ c4

R…j† ˆ

2…c2 ‡ 3†…3c2 ‡ 1† …1 7 ‡ 35c2 ‡ 21c4 ‡ c6

c†R…j

…1 c† ; …1 ‡ c†…3 ‡ c2 †

1†

c†R…j

1†

…96† 5

…1 c† ; …1 ‡ c†…5 ‡ 10c2 ‡ c4 †

…97†

and c†R… j

1†

…1 c†7 : …1 ‡ c†…7 ‡ 35c2 ‡ 21c4 ‡ c6 †

…98†

Likewise, for the Dirichlet-based operators de®ned by (55), (58), and (62), which are, respectively,     1 1 d 1 z2 …j† …j 1† b b u ˆ ‡1 b u 2 u …j 1† ; 2 in dx 2 4 ‡ z2 b u …j† ˆ

1 2

b u …j† ˆ

1 2

and





 1 d ‡1 b u …j in dx

1†

 1 4z2 2 16

 1 d ‡1 b u …j in dx

1†

  1 16 16z2 ‡ 3z4 b 2z2 u…j 2 64 ‡ 80z2 24z4 ‡ z6

 2 ‡ z2 b u …j 12z2 ‡ z4

1†

;

…99†

…100†

1†

;

…101†

we have, respectively, 3

R…j† ˆ

1 c2 ‡ 2c ‡ 5 …1 2 3 ‡ c2

R…j† ˆ

1 c4 ‡ 4c3 ‡ 14c2 ‡ 4c ‡ 9 …1 2 5 ‡ 10c2 ‡ c4

R…j† ˆ

1 c6 ‡ 6c5 ‡ 27c4 ‡ 20c3 ‡ 55c2 ‡ 6c ‡ 13 …1 2 7 ‡ 35c2 ‡ 21c4 ‡ c6

c†R…j

1†

1 …1 c† ; 2 3 ‡ c2 c†R…j

1†

…102† 5

1 …1 c† ; 2 5 ‡ 10c2 ‡ c4

…103†

and c†R…j

1†

7

1 …1 c† : 2 7 ‡ 35c2 ‡ 21c4 ‡ c6

…104†

The following comments are in order. For the Robin- and the Dirichlet based operators, the re¯ection coecient R…j† is of the form R…j† ˆ g…c†R…j

1†

‡ h…c†;

where g…c† and h…c† are rational functions of c. Since g…c† is of ®rst order in the powers of …1 g…c† ˆ O…j1

cj†;

and h…c† is of order M ˆ 3; 5; 7, speci®c for a particular operator,
h…c† ˆ O j1 cjM ;

…105† c† …106†

…107†

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

6385

the resulting jth stage solution u…j† is better than u…j 1† in the sense that R…j† , the re¯ection coef®cient associated with u…j† , is of higher order than that associated with u…j 1† (that is, R…j 1† ). The improvement continues until O …jg…c†R…j 1† j† > O…jh…c†j†, or, equivalently, until   M‡1 R…j† ˆ h…c† ‡ O j1 cj : …108† On the other hand, the Neumann based operators generate the re¯ection coecient of the form
R…j† ˆ h…c† R…j 1† ‡ 1 ;

…109†

which is superior to (105) in the sense that condition (108) is achieved in two iterations only. (Here, h…c† is of order M ˆ 2; 4; 6.) Finally, in order to assess various boundary operators, it is important to consider not only the order of R…j† , (i.e., the power of …1 c†), but also K, the associated constant:
M K ˆ lim …1 c† R…j† ; …110† c!1

…j†

where jR j ˆ O j1 Kˆ

cj

M


. Assuming that (108) is satis®ed, all the operators reach the optimal value of K,

1 : 2M

…111†

5. Iteration convergence We now investigate the convergence of iterative operators (112)   o Bu…j† ˆ Bu…j 1† ‡ a A u…j 1† ; j ˆ 1; 2; . . . ; J : ox

…112†

The analysis is limited to the plane wave case. Let us ®rst consider the Neumann-based operators:    o …j† o 1 o u ˆ ‡ A u…j 1† ; j ˆ 1; 2; . . . ; J ; ox ox 2 ox xˆ0 xˆ0

…113†

and, for focus, let us analyze the operator de®ned by (71). (The other Neumann-based operators can be analyzed along the same lines.) For the plane wave case, we have u…j

1†

ˆ f …t ‡ xc

ys† ‡ R…j

1†

f …t

xc

ys†

…114†

and u

…j†

ˆ f …t ‡ xc

ys† ‡

1 …1 c†4 1 ‡ R…j 2 1 ‡ 6c2 ‡ c4

1†




! f …t

xc

ys†;

…115† …j†

where we use the notation introduced in the previous section. Let us decompose u…j† into u0 -part which satis®es the homogeneous boundary condition …j† …116† ox u0 ˆ 0; xˆ0

and the remaining part, which is the contribution from u…j …j†

u…j† ˆ u0 ‡ bf …t ‡ xc

ys† ‡ cR…j

1†

f …t

xc

ys†;

1†

. Such a decomposition is of the form …117†

where …j†

u0 ˆ a‰ f …t ‡ xc

ys† ‡ f …t

xc

ys†Š;

…118†

6386

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

and a; b; and c are some functions of c. To ®nd them, we compare (115) and (117), and get 4

aˆcˆ

1 …c 1† ; 2 1 ‡ 6c2 ‡ c4

…119†

and bˆ

1 …c ‡ 1†4 : 2 1 ‡ 6c2 ‡ c4

…120†

Thus, on the boundary x ˆ 0, we have …j 1†

vN

ox u…j





1†

xˆ0

ˆ

cf 0 …t ‡ ys† ‡ R…j

1†

cf 0 …t ‡ ys†

…121†

and …j†

vN 

ox u…j† xˆ0 ˆ

bcf 0 …t ‡ ys† ‡ cR…j

1†

cf 0 …t ‡ ys†:

…122†

We notice that 06c6

1 2

for 0 6 c 6 1;

…123†

1 6 b 6 1 for 0 6 c 6 1; 2

…124†

and

(cf. Fig. 2), which means that the amplitudes of both the outgoing and the incoming waves are decreasing in the absolute value ( cf 0 …t ‡ ys† is mapped onto bcf 0 …t ‡ ys†, and, likewise, R…j 1† cf 0 …t ‡ ys† is mapped onto cR…j 1† cf 0 …t ‡ ys†). Now, we analyze the convergence of re¯ection coecient R…j 1† . Toward this end, let us introduce the following de®nitions [12]. De®nition 1. Let F : X ! X be a mapping. A point x 2 X is called a fixed point of F if x ˆ F …x†:

…125†

Fig. 2. Functions b and c for the Neumann-based operator.

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

6387

De®nition 2. Let …X ; d† be a metric space and F a mapping of X into itself. (Here, d is a metric.) The function F is said to be a contraction mapping if there exists a real number k; k 2 ‰0; 1†; such that d … F …x†; F …y†† < kd…x; y† 8x; y 2 X :

…126†

We now record the following theorem known as the principle of contraction mappings. Theorem 3. Let …X ; d† be a complete metric space and F : X ! X be a contraction mapping. Then F has a unique fixed point x 2 X . Moreover, if x0 is any point in X and xn ˆ F …xn 1 †, then xn ! x as n ! 1. Proof. See for example, [12].

 …j 1†

…j†

Now, we consider the map vN ! vN ; which reads

cf 0 …t ‡ ys† R…j 1† 1 ! cf 0 …t ‡ ys† cR…j 1† b ;

…127†

and we notice that it induces another map hN : ‰ 1; 0Š ! ‰ 1; 0Š; hN …U† ˆ cU ‡ c b;

…128†

U  R…j

…129†

where 1†

1:

But (128) is a contraction. Indeed, for any U1 ; U2 2 ‰ 1; 0Š; we have jh…U1 †

h…U2 †j ˆ jc…U1

U2 †j ˆ jcjj…U1

U2 †j 6

1 j …U 1 2

U2 †j:

…130†

Finally, we ®nd U , the ®xed point of (128). To this end, we solve the equation hN …U † ˆ U ;

…131†

and get U ˆ

c c

b : 1

…132†

This, together with the de®nitions of U; b and c; gives the following re¯ection coecient R : R ˆ U  ‡ 1 ˆ

…1

4

c†

4

…1 ‡ c†

:

…133†

Now, let us consider the Robin-based operators, … ox ‡ ot †u…j† xˆ0 ˆ ‰… ox ‡ ot † ‡ …ox A†Šu…j



1†

xˆ0

;

j ˆ 1; 2; . . . ; J ;

…134†

and, for focus, let us analyze the operator de®ned by (58). For the plane wave case, we have u…j

1†

ˆ f …t ‡ xc

ys† ‡ R…j

1†

f …t

xc

ys†

…135†

and u…j† ˆ f …t ‡ xc

ys† ‡

4… 1 ‡ c 2 † …1 5 ‡ 10c2 ‡ c4

c†R…j

1†

! 5 …1 c† f …t …1 ‡ c†…5 ‡ 10c2 ‡ c4 †

xc

ys†;

…136†

6388

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

which can be decomposed into …j†

u…j† ˆ u0 ‡ bf …t ‡ xc Here,

…j† u0

…

ys† ‡ cR…j

1†

f …t

xc

ys†:

…137†

satis®es the homogeneous boundary condition …j† ox ‡ ot †u0 ˆ 0;

…138†

xˆ0

…j†

and b and c are functions of c. A general form of u0 reads     1 c …j† u0 ˆ a f …t ‡ xc ‡ ys† f …t xc ‡ ys† ; 1‡c

…139†

where a ˆ a…c†. Comparing (136) with (137), we get, 4

…1 c† ; 5 ‡ 10c2 ‡ c4 4…1 ‡ c2 † bˆ …1 ‡ c†; 5 ‡ 10c2 ‡ c4 4… 1 ‡ c 2 † …1 c†: cˆ 5 ‡ 10c2 ‡ c4

aˆ

…140† …141† …142†

We also note that 4 5

…143†

4 6b61 5

…144†

06c6 and

for 0 6 c 6 1 (cf. Fig. 3). Thus, on the boundary x ˆ 0, we have …j 1† vR  … ox ‡ ot †u…j 1† xˆ0 ˆ …1 c† ‡ …1 ‡ c†R…j and …j†

vR  …

ox ‡ ot †u…j† xˆ0 ˆ b…1

c† ‡ c…1 ‡ c†R…j

1†

1†







f 0 …t ‡ ys†

f 0 …t ‡ ys†;

Fig. 3. Functions b and c for the Robin-based operator.

…145†

…146†

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

6389

which means that the amplitudes of both the outgoing and the incoming waves are decreasing. Next, we consider the map …j 1†

vR

…j†

! vR ;

…147†

and we observe that it induces another map hR : ‰0; 1Š ! ‰0; 1Š; hR …U† ˆ cU ‡ …b c†…1

…148†

c†;

where c† ‡ …1 ‡ c†R…j

U  …1

1†

:

…149†

We now show that (148) is a contraction. Indeed, jhR …U1 †

U2 j 6

hR …U2 †j ˆ jcjjU1

4 jU 1 5

U2 j;

…150†

where we have used inequality (143). Finally, we ®nd the ®xed point U of (148): U ˆ 8c 1 ‡ c2


…5 ‡

10c2

1‡c ‡ c4 †…c

1†

;

…151†

and the corresponding re¯ection coecient R : R ˆ

U

1‡c ˆ 1‡c

c†

5

…1 ‡ c †

5

…1

:

…152†

Now, let us consider the Dirichlet-based operators:   1 …j† …j 1† ot u xˆ0 ˆ ot ‡ …ox A† u ; j ˆ 1; 2; . . . ; J ; 2 xˆ0

…153†

and, for focus, let us analyze the operator de®ned by (58). For the plane wave case, we have ot u…j

1†

ˆ f 0 …t ‡ xc

ys† ‡ R…j

1† 0

f …t

xc

ys†

…154†

and ot u

…j†

0

ˆ f …t ‡ xc

ys† ‡

1 c4 ‡ 4c3 ‡ 14c2 ‡ 4c ‡ 9 …1 2 5 ‡ 10c2 ‡ c4

c†R

…j 1†

5

1 …1 c† 2 5 ‡ 10c2 ‡ c4

! f 0 …t

xc

ys†; …155†

which reduces on the boundary x ˆ 0 to …j 1†

vD

 ot u…j



1†

xˆ0

ˆ f 0 …t

ys† ‡ R…j

1† 0

f …t

ys†;

…156†

and …j† vD  ot u…j† xˆ0 ˆ bf 0 …t

ys† ‡ cR…j

1† 0

f …t

ys†:

…157†

Here, b and c are functions of c and are given below: 5

bˆ1

1 …1 c† ; 2 5 ‡ 10c2 ‡ c4

…158†

6390

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

and cˆ

1 c4 ‡ 4c3 ‡ 14c2 ‡ 4c ‡ 9 …1 2 5 ‡ 10c2 ‡ c4

c†:

…159†

We note that 06c6

9 10

for 0 6 c 6 1;

…160†

and 9 6 b 6 1 for 0 6 c 6 1; 10 (cf. Fig. 4) which means that the amplitudes of both the outgoing and the incoming waves are decreasing. Next, we consider the map …j 1†

vD

…j†

! vD ;

…161†

and we observe that it induces another map hD : ‰0; 1Š ! ‰0; 1Š; hD …U† ˆ cU ‡ …b

…162†

c†;

where U  1 ‡ R…j

1†

:

…163†

But (162) is a contraction: jhD …U1 †

hD …U2 †j ˆ jcjjU1

U2 j 6

9 jU 1 10

U2 j:

…164†

Finally, we ®nd the ®xed point U of (162): U ˆ

b c

c ; 1

…165†

Fig. 4. Functions b and c for the Dirichlet-based operator.

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

6391

and the corresponding re¯ection coecient R : R ˆ U 

…1

1ˆ

5

c†

…1 ‡ c†5

:

…166†

6. Iterative boundary operators: generalizations 6.1. Rectangular computational domains Consider the following initial-value problem for the wave equation: 8 > w r2 w ˆ 0; x 2 R2 ; t 2 …0; T Š; > > ;tt < wjtˆ0 ˆ u0 ; w pˆ u1 ; > > ;t jtˆ0 > : lim r…w;r ‡ w;t † ˆ 0; r!1

where u0 and u1 are the Cauchy data, r  from x ˆ 0: Let Xˆ…

a; a†  …

…167†

p x2 ‡ y 2 and …
†;r denotes di€erentiation in the radial direction

b; b†

S4

…168†

be the computational domain, and let Cc ˆ oX ˆ iˆ1 Cc;i be the arti®cial boundary. Here, Cc;1 ˆ fx : x ˆ a; jyj < bg, Cc;2 ˆ fx : y ˆ b; jxj < ag, Cc;3 ˆ fx : x ˆ a; jyj < bg, and Cc;4 ˆ fx : y ˆ b; jxj < ag. Also, it is assumed that the initial data u0 and u1 have compact support in X. The truncated domain problem consists of the series of initial boundary-value problems 8 …j† > u;tt r2 u…j† ˆ 0; x 2 X; t 2 …0; T Š; > > > < u…j† ˆ u0 ; tˆ0 …169† …j† ˆ u1 ; u > ;t > > tˆ0 > : …j† …j 1† Bu ˆ ri on Cc;i ; i ˆ 1; 2; 3; 4; j ˆ 1; 2; . . . ; J ; where …j 1† ri

ˆ Bu

…j 1†

‡a

M X o 2bm 1‡ ot cm mˆ1

o on

!! u

…j 1†

! M X 2bm o …j 1† / ‡a ; cm ot m;i mˆ1

i ˆ 1; 2; 3; 4; j ˆ 2; 3; . . . ; J ; …0† ri

…j 1† ˆ 0; and /m;i ; m ˆ 1; 2; . . . ; M; satisfy the following equation on Cc;i :
ott c2m oss /…jm;i 1† ˆ ott u…j 1† :

…170†

…171†

Here, on Cc;1 and Cc;3 ; o=on ˆ o=ox, o=os ˆ o=oy, respectively. Likewise, on Cc;2 and Cc;4 , o=on ˆ o=oy, o=os ˆ o=oy, respectively. Eq. (171) is appended by homogeneous initial conditions …j 1† /…jm;i 1† ˆ 0; ot /m;i …172† ˆ 0; tˆ0

tˆ0

and by appropriate corner conditions. Here, we use the corner conditions developed by Collino in [2,3]. For a typical corner, say, …x ˆ a; y ˆ b†; these conditions have the following form:   M X o o o …j 1† …j 1† cm;n /n;2 … x ˆ a; t†; m ˆ 1; 2; . . . ; M; ‡ am /m;1 … y ˆ b; t† ‡ bm u…j 1† … x ˆ a; y ˆ b; t† ˆ oy ot ot nˆ1 …173†

6392

and

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398



 o o o …j 1† ‡ am /m;2 … x ˆ a; t† ‡ bm u…j ox ot ot

1†

… x ˆ a; y ˆ b; t† ˆ

M X nˆ1

…j 1†

cm;n /n;1 … y ˆ b; t†;

m ˆ 1; 2; . . . ; M; …174†

where am ; bm ; cm;n ; are coecients which depend on bm and cm . Let w be the solution of the wave equation w;tt

r2 w ˆ 0;

x 2 X; t 2 …0; T Š;

…175†

with homogeneous Dirichlet, Neumann or Robin boundary conditions, Bw ˆ 0

on Cc;i ; i ˆ 1; 2; 3; 4;

…176†

and non-zero initial conditions wjtˆ0 ˆ u0 ;

w;t jtˆ0 ˆ u1 ;

…177†

and let E denotes the energy of w Z def 1 2 2 Eˆ jrwj ‡ jw;t j dx: 2 X Then,



Z

E;t ˆ

Cc

w;t w;n ds ˆ

0 R Cc

for B ˆ BN ; BD ; 2 …w;t † ds for B ˆ BR

…178†

…179†

and, consequently, E…t2 † 6 E…t1 †;

t2 > t1 :

…180†

Here, (179) is the standard energy identity for the wave equation, which can be obtained by taking the L2 inner product of (175) with w;t Z w;tt w;t r2 ww;t dx ˆ 0; …181† X

and integrating by parts. Inequality (180) can be viewed as a well-posedness result for the initial boundaryvalue problem (175)±(177). 6.2. Temporal approximations Any approximation of a transient problem must involve discretization both in time and space variables. Although a simultaneous discretization in all variables is possible (e.g., space-time ®nite elements), we will adopt the assumption here that the ®nal approximation is obtained by using ®nite di€erences in time and ®nite elements in space variables. Two di€erent approaches are possible. In the classical method of lines, an approximation in space variables converts the original initial boundary-value problem into a system of ordinary differential equations (ODEs) which next is discretized in time using one of many time integration schemes for ODEs. An alternative procedure known as the method of discretization in time, consists of the same two steps but done in the reverse order. By discretizing in time ®rst, the initial boundary-value problem is converted into a sequence of boundary-value problems which, in turn, give a basis for a spatial approximation and, consequently, a fully discretized scheme. Problem (169) can be solved numerically by using any of the two approaches. Here, we outline the method of discretization in time. Let us introduce a partition of ‰0; T Š into K subintervals of size Dt Dt ˆ

T ; K

tk ˆ kDt;

k ˆ 0; 1; . . . ; K;

…182†

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

and let us consider a typical step tk ! tk ‡ Dt of the following form: 8 …j† < u…j† ˆ Tutk …x†; x 2 X; tk ‡Dt …x†

…j 1† : Bu…j† ˆ B ‡ a oxo A utk ‡Dt ; j ˆ 1; 2; . . . ; J : tk ‡Dt Cc

6393

…183†

Cc

Here T is a transient operator which corresponds to a particular temporal discretization scheme applied to …j† …j† (169), utk is the (known) solution at time tk , and utk ‡Dt is the solution at tk ‡ Dt which is being sought. …J † …j† Furthermore, since the converged solution utk may provide a better starting point for iteration than utk ; …J † j ˆ 1; 2; . . . ; J 1, we replace all the iterates at tk with utk ; and arrive at the following problem: 8 …j† …J † utk ‡Dt …x† ˆ Tu > > 8 tk …x†; x 2 X;

> > > † > < < B ‡ a oxo A u…J j ˆ 1; tk ; …j† C c …184† Bu ˆ

> …j 1† o > tk ‡Dt Cc > : B ‡ a A u ; j > 1; > tk ‡Dt ox > > Cc : j ˆ 1; 2; . . . ; J : The advantage of (184) over (183) is that the iteration in (184) is limited to the layer of thickness Dt adjacent to computational boundary Cc . This follows from the ®nite speed of wave propagation, which, for the normalized wave equation, is equal to unity. 7. Numerical examples The test problem is de®ned as follows (see Fig. 5): · governing wave equation is to be solved in X1  …0; T Š; where X1 ˆ R2 n f…r; /† : r 6 r0 g;

…185†

and r0 ˆ 0:5; · the boundary conditions are de®ned below 5

u;n jrˆr0 ˆ A‰ sin…xt†Š cos…k/†‰H…t†

H…t

1†Š;

…186†

where A ˆ 0:005, x ˆ p, k ˆ 0 or k ˆ 4, and H…
† is the Heaviside step function;

y

Γc = ARTIFICIAL BOUNDARY

r φ

2a 2r0

x

VIBRATING CYLINDER ACOUSTIC FLUID

Γu = VIBRATING BOUNDARY

Fig. 5. Vibrating cylinder problem.

6394

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

computational boundary

vibrating boundary

Fig. 6. Finite element mesh. Element spectral order p ˆ 5.

Fig. 7. Solution at t ˆ 0:5. /-independent driving term, 6-term Neuman-based radiation operator.

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

6395

· the initial conditions are assumed to be zero; · truncated domain X is introduced X ˆ f…x; y† ˆ …

a; a†  …

a; a†g \ X1 ;

…187†

with a ˆ 0:55; · the radiation boundary operators are applied on the resulting computational boundary, and the results are compared with the exact solution. (Note that the exact solution for the above problem can be found by means of the Fourier±Bessel series.) The resulting bounded domain problem was solved by using Taylor±Galerkin schemes [13], with a time-step Dt ˆ 1=240, and a ®xed mesh of elements of order p ˆ 5, see Fig. 6. Note that, since the error introduced by the boundary truncation can be isolated at the continuous level only (that is, when no temporal and spatial approximations are introduced), it was necessary to select a rich enough ®nite element space, and a highly accurate approximation in time. Figs. 7 and 8 show a 3D rendering of the solution u at t ˆ 0:5 and t ˆ 1, respectively, in the segment X \ f…r; /† : p=4 < / < p=4g for /-independent boundary condition (186) (that is, with k ˆ 0 in (186)). Likewise, Figs. 9 and 10 depict a 3D rendering of the solution u at t ˆ 0:5 and t ˆ 1, respectively, in the segment X \ f…r; /† : p=4 < / < p=4g for /-dependent boundary condition (186) with k ˆ 4. These solutions are obtained by using the 6-term Neumann-based operator de®ned by (74). Figs. 11 and 12 show the graphs of the relative error ku uh k eˆ …188† kuk as a function of the number of partial fraction terms in (39). Here, k
k is the energy norm Z kuk ˆ jruj2 ‡ jot uj2 dx; X

Fig. 8. Solution at t ˆ 1. /-independent driving term, 6-term Neuman-based radiation operator.

…189†

6396

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

Fig. 9. Solution at t ˆ 0:5. /-dependent driving term, k ˆ 4, 6-term Neuman-based radiation operator.

Fig. 10. Solution at t ˆ 1. /-dependent driving term, k ˆ 4, 6-term Neuman-based radiation operator.

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

Fig. 11. Convergence of Neumann- and Robin-based operators, k ˆ 0, time t ˆ 1.

Fig. 12. Convergence of Neumann- and Robin-based operators, k ˆ 4, time t ˆ 1.

Fig. 13. Convergence Robin-based operators, k ˆ 0, time t ˆ 1.

6397

6398

A.J. Safjan / Comput. Methods Appl. Mech. Engrg. 190 (2001) 6373±6398

Fig. 14. Convergence of Neumann-based operators, k ˆ 0, time t ˆ 1.

u is the exact solution of the unbounded domain problem, uh is the solution obtained by using iterative boundary operators, and the graph is plotted in the log±log scale. As can be seen, the boundary operators produce results convergent toward the exact solution. Finally, depicted in Figs. 13 and 14 are the graphs of the relative error as a function of the number of iterations, and a very fast convergence can be observed. In conclusion, the results of these (and other) numerical experiments fully support the theoretical results, and suggest a high potential of the new families of iterative radiation boundary operators. Acknowledgements Support of this work by the National Science Foundation (NSF) under the grant #9713893 and by the Oce of Naval Research (ONR) under the grant #N00014-99-1-0778 is gratefully acknowledged. References [1] A. Bayliss, E. Turkel, Radiation boundary conditions for wave-like equations, Commun. Pure Appl. Math. 33 (1980) 707.  e pour des Modeles de Propagation d'Onde dans des Domains Rectangulaires, [2] F. Collino, Conditions Absorbantes d'Ordre Elev INRIA Rapports de Recherche, No. 1790, 1992.  e pour l'Equation des Ondes 3D, INRIA Rapports de Recherche, [3] F. Collino, Conditions aux Limites Absorbantes d'Ordre Elev No. 2932, 1996. [4] B. Engquist, A. Majda, Radiation boundary conditions for the numerical simulation of waves, Math. Comput. 31 (1977) 629. [5] B. Engquist, A. Majda, Radiation boundary conditions for acoustic and elastic wave calculations, Commun. Pure Appl. Math. 32 (1979) 313. [6] M.J. Grote, J.B. Keller, On nonre¯ecting boundary conditions, J. Comput. Phys. 122 (1995) 231. [7] M.N. Guddati, J.L. Tassoulas, Continued fraction absorbing boundary conditions for the wave equation, J. Comput. Acoust. 8 (2000) 139. [8] T. Hagstrom, On high-order radiation boundary conditions, NASA Contractor Report 198404, ICOMP-95-15. [9] R.L. Higdon, Absorbing boundary conditions for di€erence approximations to the multi-dimensional wave equation, Math. Comput. 47 (1986) 437. [10] J.B. Keller, D. Givoli, Exact nonre¯ecting boundary conditions, J. Comput. Phys. 82 (1989) 172. [11] E.L. Lindman, Free-space boundary conditions for the time dependent wave equation, J. Comput. Phys. 18 (1975) 66. [12] J.T. Oden, L.F. Demkowicz, Applied Functional Analysis, CRC Press, Boca Raton, 1996. [13] A. Safjan, J.T. Oden, High-order Taylor±Galerkin methods for ®rst-order linear hyperbolic systems, J. Comput. Phys. 120 (1995) 206. [14] A. Safjan, Highly accurate non-re¯ecting boundary conditions for ®nite element simulations of transient acoustics problems, Comput. Methods Appl. Mech. Engrg. 152 (1998) 175±193. [15] L.N. Trefethen, L. Halpern, Well-posedness of one-way wave equations and absorbing boundary conditions, Math. Comput. 47 (1986) 421.