- Email: [email protected]
Non-reflecting boundary conditions for elastic waves
WAVE M O T I O N 12 (1990) 261-279 NORTH-HOLLAND
261
N O N - R E F L E C T I N G B O U N D A R Y C O N D I T I O N S FOR ELASTIC WAVES Dan GIVOLI Department of Aeronautical Engineering, Technion, Haifa 32000, Israel
Joseph B. K E L L E R Departments of Mathematics and Mechanical Engineering, Stanford University, Stanford, CA 94305, USA Received 2 March 1989, Revised 16 October 1989
An exact non-reflecting boundary condition is devised for time-harmonic two-dimensional elastodynamics in infinite domains. The domain is m a d e finite by the introduction of a circular artificial boundary on which this exact condition is imposed. In the finite computational domain a finite element method is employed. Numerical examples are presented in which the accuracy a n d efficiency of the method using the exact non-local boundary condition are compared with those of methods based on approximate local boundary conditions. The method is also used to solve problems in large finite domains by reducing them to smaller domains. In addition, local boundary conditions are derived which are exact for waves with a limited n u m b e r of angular Fourier components.
1. Introduction To solve a wave problem numerically in an unbounded domain one usually introduces an artificial boundary 90, thereby defining a finite computational domain O. On 90 one must impose some boundary condition which should allow waves f r o m / 2 to pass through 90 without generating spurious reflections. We shall find an exact boundary condition having this property, for problems involving two-dimensional time-harmonic elastic waves. In O we shall employ the finite element method. For problems of wave propagation there has been a considerable amount of previous work to devise appropriate artificial boundary conditions. The time dependent scalar wave equation and the reduced wave equation were considered by Engquist and Majda [1], Bayliss and Turkel [2], Feng [3], Ting and Miksis [4], Higdon [5], MacCamy and Marin [6], Goldstein [7] and Keller and Givoli [8]. In references [2] and [5] approximate local boundary conditions are derived. In [1] and [3], exact nonlocal conditions are derived and approximated by local ones, whereas in [4], [6], [7] and [8] exact non-local boundary conditions are employed. To get an exact boundary condition for elastic waves we use the procedure presented in [8] for the simpler case of the reduced wave equation. Artificial boundary conditions for elastic waves have been considered mainly in the context of geophysics. In two dimensions the artificial boundary 90 is traditionally composed of straight line segments and the computational domain O is usually a rectangle. Often the original domain is the whole space, possibly with some holes in it, which corresponds to wave propagation deep inside an elastic body. We now describe some of the previous methods. We denote the coordinates normal and tangential to 90 in two dimensions by Xl and x2 respectively. We also let Ul and u2 denote the displacements and T1 and T2 the tractions in these directions. 0165-8641/90/$03.50 @ 1990--Elsevier Science Publishers B.V. (North-Holland)
262
D. Givoli, ZB. Keller / Non-reflecting boundaries
Lysmer and Kuhlemeyer [9] suggested the following boundary condition on ~ : Oul =
apcL O----t T,,
bpCT~t= T 2 .
(1)
Here p is the mass density, t is time, CL is the longitudinal wave speed, CT is the transverse wave speed, and a and b are dimensionless parameters. The parameters a and b were chosen to minimize the reflected energy for an incident plane wave which hits the boundary ~ at a given angle of incidence. The choice of a and b was considered separately for an incident longitudinal wave, for an incident transverse wave, and for a surface wave. It was suggested that a = b = 1 is a good choice in general. In the computational domain ~ a finite element scheme was used. Smith [10] showed how the reflections from a piece-wise planar boundary ~ can be eliminated by adding together the numerical solutions of several problems, in each of which a certain combination of Dirichlet and Neumann boundary conditions is used on ~. If reflections are to be eliminated exactly on n plane surfaces, 2 n such solutions must be added together. This procedure is applied separately to dilatational waves and to surface waves. Cerjan et al. [11] and Sochacki et al. [12] each presented a method in which the amplitudes of ul and u2 are gradually reduced in a strip of nodes adjacent to ~. Thus, the solution is artificially damped in the vicinity of ~. On ~ , the usual Dirichlet or Neumann boundary condition is used. Engquist and Majda [13] found an exact non-local boundary condition on ~ in terms of pseudodifferential operators, and then approximated it by local boundary conditions. The simplest of these approximate conditions is -
CL
Ul = 0,
--
CT
U2 =
0.
(2)
This condition is perfectly absorbing for plane waves at normal incidence. Their next boundary condition involves a linear combination of the operators 02/0x, at, 02/0t 2, 02/aX2 0t and a2/ox 2. These boundary conditions are used in a difference scheme employed in a rectangular domain. A special procedure is suggested to avoid instabilities near the corners. Some numerical experiments with these boundary conditions are presented in Clayton and Engquist [14]. Scandrett et al. [15] obtain time-harmonic solutions of the equations of elastodynamics by solving numerically the time-dependent problem and using the limiting amplitude principle. They derive an approximate time-dependent boundary condition on ~ , which is equivalent to the first-order Engquist and Majda condition. An analogous boundary condition for the three-dimensional case is also obtained. A finite difference scheme in time and space is employed. Higdon [16] suggested using a difference scheme in ~ with a boundary condition on ~ of the form
Here the fl~ are parameters which can he adjusted so that the boundary condition is perfectly absorbing for an incident plane wave which hits ~ at a given angle. In the case m = 1 with normal incidence, (3) reduces to the Engquist-Majda condition (2). If m is large enough, the resulting high-order boundary condition can be made perfectly absorbing for a number of chosen incidence angles and for both longitudinal and transverse waves. Surface waves were not treated. The boundary conditions (1), (2) and (3) can be used in solving either the time-dependent elasticity equations, or those for time harmonic waves with frequency ~o. In the latter case the equations involve
D. Givoli,J.B. Keller / Non-reflectingboundaries
263
the space variables alone. Then (1)-(3) can still be used with auj/~t replaced by itouj. Since they are approximate, they generate spurious reflection from ~. The magnitude of the reflected wave depends upon the specific problem. It is hard to estimate because of the presence of various types of waves. We shall derive an exact non-local boundary condition on ~ which will eliminate all reflections for any type of wave. It has the form - Ti(x) = Mijuj(x)
m .=0 ~"
f~ m~(x, x') us(x') dx'.
(4)
The summation convention for repeated indices is used here and elsewhere. This boundary condition is analogous to the condition introduced in [8] for the reduced wave equation. In Section 2 we derive the boundary condition (4) for time-harmonic waves in two dimensions, when the artificial boundary ~ is a circle. We consider problems of plane stain and problems of plane stress. In Section 3 we describe the finite element scheme employed in/2. A similar formulation is discussed in more detail in [8] in the context of the reduced wave equation. We also show how the finite element method with the boundary condition (4) can be used to solve efficiently problems in large finite domains by eliminating an interior circular region. Numerical experiments are presented in Section4. The results obtained by using the exact boundary condition are compared with those obtained by using the approximate Lysmer-Kuhlemeyer boundary condition (1) and by using some of the approximate boundary conditions obtained in the appendix. We compare both the accuracy and the efficiency of the various schemes. In the Appendix we localize the exact boundary condition and obtain a sequence of local approximate boundary conditions. The M t h boundary condition in that sequence is exact for waves that consist of the first M Fourier components in the angle variable 0. We also obtain analogous local boundary conditions for the reduced wave equation.
2. An exact non-reflecting boundary condition We consider the equations of linear elastodynamics in an infinite two dimensional region. To make the region finite, we introduce an artificial boundary ~ which is a circle of radius R. Our goal is to derive an exact boundary condition of the form (4) on ~. We shall do so by finding an exact solution in the region exterior to ~. In that region the medium is assumed to be linear, homogeneous, isotropic and free of body forces. Then outside ~ the elastodynamic equations in Cartesian coordinates are
p (3, =~ u,.jj+(,+~) u.,.
(5)
Here U~ is the displacement in the i direction, i = 1 or i = 2, p is the mass density, and a comma indicates partial differentiation. The constants /~ and A are related to Young's modulus E and to Poisson's ratio v via
E ft = 2(1+ v) ;
f rE~(1+ v)(1 A = /[ v E / ( l _ v 2 ) '
-
2v),
plane strain, plane stress.
(6)
To treat time-harmonic solutions we set
Uj(x, t) = uj(x) e i~'',
(7)
substitute (7) into (5) and obtain p. u~,~+ (A +/z) ujd~+~o2put =0.
(8)
D. Givoli, J.B. Keller / Non-reflecting boundaries
264
We now use the Helmholtz decomposition theorem to write ui in terms of two potentials 4) and 0: u l = 4<1 + 0,2,
(9)
u2 = 4<2- 0,1.
Then (8) will be satisfied if 4) and 0 satisfy V2~b+ k2~b = 0,
V20 + k20 = 0.
(10)
Here k L and kT are the longitudinal and transverse wave numbers: k2=to 2
P A +2/.~'
k2=to 2 p . /~
(11)
At infinity the Sommerfeld radiation conditions must hold:
!ira rl/2(C~,,--ikLCk)=0,
!im rl/2(0.,--ikT0) =0.
(12)
The general solutions of (10) and (12) are, with a. and b. arbitrary constants,
cb(r, O)= ~ a.H~l)(kL r) e i.o ,
O(r, O) = ~ b.Ht.')(kTr)
n=--oo
e i"°.
Here H(,]) is the Hankel function of the first kind. Since the displacements u~ and polar coordinate form of (9),
U,.=qb.r+lo.o,
(13)
n=--oo
uo are
Uo=--O.,.+I~,o,
given by the
(14)
equations (13) and (14) yield
Ur= ~ [a.kLH(1),(kLr)+b inH(,)(kTr)] e i.o, n = --oo
r
(15)
Now we express the constants a. and b. in terms of the values u,(R, 0) and uo(R, boundary r -- R. To this end we first expand u,(R, 0) and uo(R, O) in Fourier series:
[1I?
=
e -i~°' u.(R, 0 ' ) d O '
]
O) on
the artificial
ei'S,
n =-oo
(16) u o ( R , O) = . . ._~.
Then
e -i"°'
us(R, 0') dO'
e i'°.
we evaluate (15) at r = R and compare it with (16) to obtain
a. -- 2~rA. ~ H.TUo(R) 1 I o " e-i"°' [ - k T H n T' U r ( R ) O ' ) - -in
0') ] dO',
(17) 2erA.
Here
, Uo(R. 0') - - ~ e -i.o' kLHnL
HnL u,.(R, 0')
dO'.
265
D. Givoli,J.B. Keller / Non-reflectingboundaries An
~
n2 I f --kLk'rHnLHnT -I-_R_ 2 HnLHnT
(18)
and
H~L = H(.'(kLR),
H.T
dH("l)(x)
HrnL
=
H'-r
dx
H(.'(kTR),
(19)
dH("l)(x)x=k~a'
dx
X=kL R '
Substituting (17) into (15) gives
-uAr, O)= ~
[sTdr, O,O')u,(R,O')+sT~(r,O,O')uo(R,O')]dO',
n=--oo
-uo(r,O)=
I
(20)
~
[s~l(r,O,O')u~(R,O')+s~2(r,O,O')uo(R,O')]dO',
n=-oo
where
r
]
in<°-°')L kLkrH,TH~l),(kLr)_._g_ s'il(r, O, 0 ' ) - e2~a-------~ ~r H~LI-I~l)(krr) , n
_
_
1?!
t
(1)
kLHnTHO)'(kLr)----r kLHnLHn (krr) ,
$12(r, 0, 0') -- 2~rA.
(21) [ ~ kTH.TH~ , ~)(krr)-'-~ in kTH.LH~),(kTr)], - 2,rra,, sEt(r, O, O') - ein(#-°')
kLkTHnLH" (kTr) "
s~2(r, 0, 0 ' ) = 2,rrA.
Our goal is to express the tractions T,(R, O) and To(R, O) in terms of u,(R, O) and uo(R, 0). We first consider plane stress, for which
(
v
v )
T~ - (1 - v 2) u,,, +rUO.O+rU, ,
E (lu,,o+uo._luo)
To = 2(1+ v)
.
(22)
We use (22) to calculate the tractions from (20), set r = R, and after some algebra we get
- T,(R, O) = E'
[ M d o , 09 ur(R, o')+pl%(o, 09 uo(R, o')]R dO',
rl=O
(23)
-To(R,O)= E'
[p~l(O,O')u,(R,O')+p2~(O,O')uo(R,O')]RdO'.
n=O
The prime on the sum indicates that a factor of 1/2 multiplies the term with n = 0. We eliminate H~L and H~-r from the kernels p~ by using the identity
H~l)"(x)=(~2-1) H~l)(x)-lH~l)'(x)'x
(24)
266
D. Givoli, J.B. Keller I Non-reflecting boundaries
Then we finally obtain p T , ( o - o ' ) =( l-+-Eu-) cos ,'itn(O RA,,
,
pT~(o - o') = (1 + v)
o,)r|
2 , t. 3 " k L k ' r H ' L H ' r +
~g~,~.
k L k v H,LH,T--~' ' R
n2 HnLHnT],
k~k,~H'~H',;
]. (25)
p~,( O - 0') = -p'~2( O - 0 %
E c o s n ( 0 - 0 ' ) [ l l " 1.2u, u ±1 n2 ] (1 + u-----) ~-~-~ L~'~L'~V'"nLt'nTT'RkLkTH~nLH~nT----~ H n L H n T _ ,
p~2(O - 0')
3' = 1/(1 -- l')
(plane stress).
(26)
Equations (23) and (25) constitute the exact boundary condition on ~ in polar coordinates. To obtain the boundary condition in rectangular coordinates, as in (4), we transform T,, To, u,, uo in (23) into T~, T2, ul, u2. This leads to the following relation between the polar kernels p~ and the rectangular kernels rots: (27)
m~(O - 0') =p~a(O - 0') V~g(O) Vj,(O').
Here Vl,(X) = V22(x) = cos x,
V~2(x)= - s i n x,
V21(x) = sin x.
(28)
It is easy to see that the kernels m ~ ( O - 0') possess the symmetry (29)
m ~ ( O - 0') = mj~( O ' - 0).
This property is essential to preserve symmetry in the finite element scheme to be employed in the computational domain/2, as we shall see in Section 3. Equations (25) and (26) were obtained for plane stress. Now we replace (22) by the corresponding relations for plane strain. Then we find that the formulae (25) for p~ are the same, with the definition (26) of 3' replaced by 3' =
l-z,
(plane strain).
1-2v
(30)
3. Finite element formulation
We shall now combine the boundary condition (4) with the finite element method in the computational domain O. Since this formulation is analogous to that presented in [8] for the reduced wave equation, it will be described briefly. The computational domain ~2 is bounded in part by the circular artificial boundary ~, and in part by some interior boundary F = F s w Fh. T h e problem in 12 is: ( CoklUk.r). j + tO2pUi = 0 ui = gi
in O,
(32)
on Fs,
Ti = CijklUk.tnj = hi - Ti = Mijuj
on ~.
(31)
on Fh,
(33) (34)
D. Givoli, ZB. Keller / Non-reflecting boundaries
267
Equation (31) is the generalization of equation (8) to a non-isotropic and non-homogeneous medium, which may be present in O. In (33), nj is the component of the unit outward normal to Fh. Equation (34) is the boundary condition (4). To obtain a weak form of the problem (31)-(34), we introduce the two sets
S¢= { u l u e n l ( 1 2 ) x n l ( o )
& u,= g, on Fg},
6eo={wlweHl(O)xH'(O)
& w , = 0 on rs}.
(35) Here H~(O) is the Sobolev space containing functions which are square-integrable and which have square-integrable first derivatives. We also introduce the notation
v(i,j) I- 3( vi,j + vj,i).
(36)
Then the weak form of the problem is: Find u ~ 5e such that for all w ~ Seo
(37)
a(w, u) + b(w, u) = (w, h)r where
f.
w,h,
,..
,39.40,
The finite element method is obtained by approximating the weak form (37). The domain D is discretized into a finite number of elements, and each element is associated with a finite number of nodes. Then u and w are approximated by the sums u h and wh: uh(x) = ~ dA,NA(X)+ E gA,NA(X)',
gA,=g,(XA),
(41)
AE~g
A~71
w~(x)= ~ CAiNA(X).
(42)
A~rl
Here rig is the set of nodes on the boundary Fg, 7/is the set of all other nodes, NA is the shape function associated with node A, :cA is the location of node A, and CAi and dai are constants. We now substitute (41) and (42) into (37). This results in the finite element matrix form of the problem for the vector d:
Kd=F,
K=K°+K
K a = [K~4iBj],
b,
(43,44)
K b = [K~,aj],
K~AiBj = a( NAei, NBej),
d = {dsj},
K bA~Bi= b( NAe~, N ~ ),
Far = (Naei, h)r -~, Y. gBia(Nae~, Na~).
F = {FA,},
(45) (46) (47)
J BE~ B
Here e~ is the unit vector in the x~-direction. On the left side o f (46), Ai and Bj are the positions in the matrix corresponding to degree of freedom i at node A and degree o f freedom j at node B. A similar meaning applies to the Ai on the left side of (47). Equation (44) shows that the effect of the boundary condition on the standard finite element scheme is just the inclusion o f K b in the matrix K in (43). The entries of the matrix K b can be written explicitly
D. GitJoli, J.B. Keller / Non-reflecting boundaries
268 as
b
K A~aj= b( Nae,, NBe1) = = ,=o ~
L NAM~jNBd~3
f~ f~ NA(X)m~(x,x')Ns(x')dx' dx;
To evaluate the double integral in (48) we note that 4
(48)
m~(O-0') given by (27) can be written in the form
4
m~(O- 0') = E E a~k' F~,(O) F'~(O') k=l
(49)
I =1
where
F~(O)=cosnOcosO, F~(O) = sin nO cos O,
F~(O)=cosnOsinO, F~(O)= sin nO sin O,
(50)
and the a~ kl are constants. For example, all,,1 = a~,33 is the coefficient multiplying cos equation of (25). Therefore (48) gives .=o k=, ,:, a,j ~ j ~
NA(x) Fg(x) dx
N.(x) FT(x) dx .
n(O- 0') in the first
(51)
b Thus, only one-dimensional integrals have to be computed to form K b. Of course, KA~Sj is nonzero only b if both nodes A and B lie on the boundary ~ . Also, to compute KA~nj in practice, the infinite sum in (51) has to be truncated after a finite number of terms. Other computational aspects of the method are discussed in [8]. The matrix K b is symmetric. This follows from (48) and from the symmetry of the kernels m~ in (29). Thus the exact boundary condition preserves the symmetry of the standard finite element formulation. After the finite element solution in g2 has been found, the solution outside g2 can be determined from equations (20). The finite element solution Uhr(R,0') and uho(R, 0') can be used to replace u,(R, 0') and uo(R, 0') in those equations. The present finite element method can also be used to solve problems in large finite domains. To see this, we consider a finite two-dimensional domain ~ with an irregularly shaped outer boundary 0~. We suppose that the medium is homogeneous and isotropic, except perhaps in a region near a~. To obtain a small computational domain g~ we cut out of ~ a circular domain D within which the medium is homogeneous and isotropic. The boundary of D is denoted ~. Thus the computational domain ~ = ~ - D is b o u n d e d internally by ~ and externally by aR. By finding an analytic solution u in D, we derive an exact boundary condition on ~ of the form (4). Then we combine this boundary condition with the finite element method in £2. The finite element formulation is identical to the one given before. To get the exact boundary condition we go through a procedure analogous to the one outlined in Section 2. Then we find that equations (8)-(30) hold with the Hankel function H~1) replaced everywhere by the Bessel function .In. Since the boundary condition is real in this case, the whole problem is real.
4. Numerical experiments We now compare the results obtained by using the finite element scheme employing the exact boundary condition with results obtained by using the scheme employing local approximate boundary conditions
D. Givoli, J.B. Keller / Non-reflecting boundaries
269
and with the exact solution. We consider first the case where the domain is infinite. The local boundary conditions that we use are the one by Lysmer and Kuhlemeyer (1) with the parameters a,= b = 1, and two local conditions obtained in the appendix, namely the condition (A11) and its asymptotic counterpart (A14). The local boundary conditions (2) and (3) are not compatible with the finite element method, because they involve the N e u m a n n operator Ou~/On rather than the traction T~. We consider an infinite plate with a circular hole of radius a. Now from (15) it is easy to see that the following are two solutions of the elastodynamic equations:
Solution 1 (longitudinal mode): Ur(r, O) = k L H ~ ( kLr) c o s nO,
uo(r, O) =
n_ Hn(kLr) sin nO. r
(52)
Solution 2 (transverse mode): n
u,(r, O) = - - H , ( k T r ) sin nO, r
uo(r, O) = - k T H ' ( k T r ) cos nO.
(53)
In order to have (52) or (53) be the exact solution for the plate with a hole, we use as boundary condition on the hole boundary that obtained by setting r = a in (52) or (53). We choose a = 0.5 and we set the artificial b o u n d a r y ~ at r = R = 1. The other parameters are to = 1, p = 1, E = 1 and ~, = 0.3. The mesh is composed of 120 bilinear quadrilateral elements. It is shown in Figure 1. In the exact boundary condition (23) we use the first five terms. In fact, only the term corresponding to n in (52) or (53) contributes to the solution.
Fig. 1. Finite element mesh for the infinite plate with the circular hole.
We first solve the axisymmetric problem with the longitudinal mode. This case corresponds to n = 0 in equations (52). Figure 2 shows the real part o f the displacement u~ on ~. Five solutions are compared: the exact solution, the one obtained with the exact boundary condition, and the three solutions obtained by using the local boundary conditions (1), ( A l l ) and (A14). We see that the exact boundary condition and its local version ( A l l ) give results which coincide with the exact solution. The local condition (A11) is indeed expected to give identical results to those obtained with the non-local exact b o u n d a r y condition because ( A l l ) is exact in the axisymmetric case. The solution with the asymptotic condition (A14) is about 8% off the exact solution, whereas the L y s m e r - K u h l e m e y e r solution is almost 50% off at the peaks. Our next problem corresponds to n = 1 in equations (52), still in the longitudinal mode. Figures 3 and 4 present the real parts of u I and u2 respectively along ~ . In both figures the results obtained with the
D. Oivoli, J.B. Keller / Non-reflecting boundaries
270 Real
. ., . , ,. . ,.
ul
0.6"
...... y.' .....................
0.5" 0.4" 0.3 ~ 0.2 ~ 0.1"
...////'//'"
""''".....,,%.
0.0 : -0.1 -0.2 -0.3 -0.4 -0.5 -0.6
~
.....
-0.7 -0.8 -0.9 -1.0 '
.
.
.
.
.
.
.
'
120
I
'
'
'
'
.
.
.
.
'
I
130 node
LEGEND
'
'
'
'
'
'
'
'
'
140
.
.
.
.
.
.
.
.
.
160
number
Lys.-Kuh. exact b.c. local N=O
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
I
150
- ...............
asymp, b.c exact sol.
Fig. 2. The infinite plate problem: comparison of results obtained with various boundary conditions for longitudinal axisymmetric waves. The real part of the displacement u I is shown along the artificial boundary r = 1.
Reol
ul 2.2'
1.8 1.6
! 1.4
i
1.2
/
1.0
0.6
i
0.4
x ",
/
t/
,
\
/
\
/,
I
i
i
i
///"
//.
/
,, \
',.
"',. \
•
".
\\
.
.
.
.
.
.
.
.
.
.
///"
/ .,,"
., x
",,\ \
.
\
-...x.--.;i I
.
.
.
.
.
.
.
.
.
130
I
. . . . . . . . .
140 node
LEGEND
i
/ , ///
",,
"-.'2~- .>.." 120
/// "'\'k /" '\ .....'""...,.
/..,.'"'"'"',,,,, ;/ "...
2.0
.................................. L y s . - K u h . ........... exact, b.c. . . . . . . . . local 14=0
I
.
.
.
.
.
.
.
.
.
150
160
number - ...............
asymp, b.c exact sol.
Fig. 3. The infinite plate problem: comparison of results obtained with various boundary conditions for longitudinal 'cos 0 waves'. The real part of the displacement ul is shown along the artificial boundary r = 1.
D. Givoli, J.B. Keller / Non-reflecting boundaries Real
271
u2 1.0
.:..-'""'.,..
0.9 0.8
.......''"".
/."
/ /
)\_
0.7
0.6
/-
\
0.5
0.4 0.3 0.2 0.1 O.O -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7
\
/
/
i :
!"
x
~.,~
-~. -
"x
/ ~.~ ":/"
"~
.i/
"~:X
-\
./
'..
\~,~
/-/
/
\-\
/
.,, .~ /
x\
y/:
'"",.,,,,,
/." t.\.
". /.":: "~........."
".........f
-0,9 ~1.0 . . . .
120
'
'
'
'
'
1
'
'
'
.
.
.
.
.
.
130
.................................. L y s . - K u h . . . . . . . . . . . . e x a c t b.c.
........
I
,
'
'
,
. . . .
140 node
LEGEND
,,
/
/
\
-0.8
.
j':"
local
'
1
'
'
,
,
. . . .
150
,
160
number - ...............
asymp,
exact
b.c sol.
M=O
Fig. 4. The infinite plate problem: comparison of results obtained with various boundary conditions for longitudinal 'cos 0 waves', The real part of the displacement u 2 is shown along the artificial boundary r = 1.
exact b o u n d a r y condition deviate from the exact solution only slightly. On the other hand, the results using the three local approximate boundary conditions deviate considerably. The condition (1) gives unacceptable results. It is a coincidence that the local condition ( A l l ) is even worse than its asymptotic counterpart (A14); it so happens that the error due to the asymptotic approximation cancels part of the finite element error. N o w we solve the problem which corresponds to n = 1 in equations (53), in the transverse mode. In Figure 5 the real part of u2 on ~ is shown. The exact boundary condition again gives excellent results whereas the approximate conditions do not. We note that the results obtained from the conditions (A11) and (A14) tend to err more symmetrically with respect to the exact solution than condition (1). This is apparent both in Figure 3 and Figure 5. Next we move to a problem in a finite domain. We consider time-harmonic waves in a thin circular plate of radius a. The displacements on r = a are prescribed. We delete a circular domain of radius R. Then we impose the exact boundary condition on r = R. The computational domain is the ring R < r < a. Equations (52) and (53) are still solutions of the ¢lastodynamic equations if the Hankel functions H are replaced by the Bess¢l functions J. On r = a we impose the boundary condition obtained by setting r = a in (52) or (53) (with H replaced by J). We choose a = 1 and R = 0.5. The other parameters are the same as for the infinite domain problem. The mesh shown in Figure 1 is used again, but this time the non-reflecting boundary condition is imposed on the inner circle. We compare the exact solution with that obtained by using the exact boundary condition (A11) (with H replaced by J) and the local asymptotic condition (A25). Figure 6 depicts ul along the artificial b o u n d a r y r = 0.5 in the axisymmetric longitudinal mode. As in the infinite case, solutions with the exact boundary condition and the local condition (A11) coincide with the exact solution. The asymptotic condition gives bad results. In fact, the argument kLR of the Bessel
D. Givoli, J.B. Keller / Non-reflecting boundaries
272
Ree] u2 1.2 1.0
......'"""'",,...,,
0.8
.\
0.6
\
/" .;" ..."
/
/,..'"'"'",,,,.., ~\\ "..\
\
i
0.4
\
/
0.2
i / ":..... " • ",
"--'-:~" \ '"",.
0.0 -0.2
-"-
-0.4
,,.
~ " - ~ x " .x.- ~
i .,/'
x
,'"
~,.
i.~ , / , :
i ,.',"
~ . \ . ~ . . .. //'/z
-0.6 .....f,-
-0.8
".......,."
-1.0 .
.
.
.
.
'
'
'
'
120
I
'
'
'
'
'
.
.
.
.
I
130
.
.
.
.
.
.
.
.
.
I
140
.
.
.
.
.
.
'
1.50
'
'
160
node number ..................................Lys.-Kuh. ........... exact b.c. ........ l o c a l N-O
LEGENO
- ...............
asymp, b.c exact sol,
Fig. 5. The infinite plate problem: comparison of results obtained with various boundary conditions for transverse 'cos O waves'. The real part of the displacement u2 is shown along the artificial boundary r = 1.
ul 0.7 0.6
............................... .....",..
0.5 0.4
...""...... "",., "
0.3 0.1 0.0 -0.1 -0.2 -0.3
~
//"
~
/
"",.
.:/
""',.,",,.
-0.4
../""
...........................................
-0.6-0"P -0.7 -0.8
-0.92 -1.0 2
. . . . . . . . .
120
I
. . . . . . . . .
130
I
. . . . . . . . .
140
I
.
.
.
.
.
.
.
.
.
150
160
node number LEGEND
.................................. esymp, b.c exact s o l .
exact b.c. l o c a l H=O
Fig.&Thecircu~arp~r~b~m:c~m~ari~n~£resu~ts~btainedwithva~usb~undaryc~nditi~ns~r~n~tudina~axisymm~t~c waves. T h e d i s p l a c e m e n t u t i s s h o w n a i o n g t h e a ~ i f i d a l b o u n d a r y r = 0 . 5 .
D. Givoli, J.B. Keller I Non-reflecting boundaries
273
functions has the value 0.48 here, and for this value the asymptotic approximation invoked in (A25) is very bad. On the other hand, in the previous problem the argument was twice as large and, more important, the asymptotic approximation of the Hankel functions is quite good even for such small arguments. In Figure 7 the displacement u2 is shown in the case n = 1 of the longitudinal mode. The exact boundary condition is again superior by far to the approximate ones. We have measured the CPU times needed to obtain the results above with the various boundary conditions. After averaging the times used in all the cases considered, we estimate the running time in Table 1. The time corresponding to the Lysmer-Kuhlemeyer condition (1) is normalized to 1. We see that the computation using 5 terms of the exact boundary condition takes about 10% more time than that using the local conditions. This is a small price to pay for the tremendous gain in accuracy obtained in the preceding examples. Finally, to demonstrate a possible application of the method when an explicit analytic solution cannot be obtained, we consider the following problem. An infinite plate with a circular hole is given. Young's modulus of the plate material is Ep = 10. Into the circular hole an inclusion with a square hole is fitted, which is of a material ten times 'softer' than that of the plate, i.e. with Young's modulus of E I : 1. Periodic expansion and contraction of the square hole produces elastic waves in the inclusion and in the plate. u2 0.26" 0.21
/
0.16 0.11 .,i
"\ \
0.06
\ ,, / . . . .
.\
i,/
',.\
i
\",,
,/"
0.01
.,,'/
"\
/
"..'\
/ :"
-0.09
/
,, \
]'
/'~ , , / /
""'k
-0.14 -0.19 -0.24
//
'-'.,
/] -0.04
•
.........
120
k'.,X,,, / :/" ',.,,.,..,,./ I .........
I .........
130
LEGEND
..................................
asymp,
execL
"
'k,,,., // ",.,.,...,../' I .........
140 node
/ ' " ~--
150
160
number
b.e sol.
................ "..........
exact local
b.c. N=O
Fig. 7. The circular plate problem:comparison ofresults obtained with various boundary conditionsforlongitudinal'cos 0 waves'. The displacement u2isshown aiongthe a~ificialboundaryr=0.5. Table 1 Normalized CPU times for the finite element schemes employing various non-reflecting boundary conditions B.C.
Lys.-Kuh. B.C. (1)
Asymptotic B.C. (A14)
Local M = 0 B.C. (All)
Exact B.C. 2 terms
Exact B.C. 5 terms
Normalized time
1
1
1.02
1.07
1.11
D. Givoli, £B. Keller / Non-reflecting boundaries
274
The configuration is shown in Figure 8. The parameters are: v = 0.3 and p = 1 for both materials, to = 10, the radius of the inclusion is 2.5 and the edge of the square hole is 2. On the square boundary the normal displacement is prescribed to have the value 1, whereas the tangential displacment is unspecified. The artificial boundary ~ is set at r = R = 3. Thus of the plate material only a ring of thickness 0.5 is included in the computational domain. We solve this problem using the exact boundary condition on with 5 terms in the expansion (23)• The mesh includes 128 bilinear quadrilateral elements: the inclusion comprises three 'rings' of 32 elements each, and the domain of the plate is composed of an additional ring of 32 elements. Figure 9 shows the result for the real part of u~ along the radial line y = 0, starting from r = 1 at the middle of the side of the square hole and ending at r = 3 on the artificial boundary. The point r = 2.5 is •
.
°
.
.
.
.
•
.......
• " " : lZli kl'O:
: i i ! ! i i i i-:~'~": . . . . . . .
~: i ! i i i i
f..
.
. % . . . . . . .
::::::/: . . . . .
:,~:::::: J..
.
:::::': :::::r:
i iii i~:
.~. . . . . . .
:i::::: :l:::::
i,'iiiii
"::::::.~:..
.:/::::::"
Fig. 8. The set up of the plate-inclusion problem.
Real
ul 1
-1
-2 '
'
'
'
'
'
'
'
'
I
'
'
'
'
'
'
'
'
'
2 r
Fig. 9. The plate-inclusion problem. The real part of the displacement u~ is shown along the radial line y = 0.
D. Givoli, J.B. Keller / Non-reflecting boundaries
275
at the plate-inclusion interface. The displacement u2 is zero on that line due to symmetry. In Figure 10 the imaginary part of ul is shown along the artificial boundary in the first quadrant. Due to symmetry the distribution is similar in the other quadrants. In fact, the symmetry can be exploited by discretizing only a quarter of the plate and by using symmetry boundary conditions. If that is done the exact non-reflecting boundary condition must be modified according to the guidelines given in [17].
Appendix
Approximate local boundary conditions Sometimes it may be desirable to approximate the nonlocal boundary condition (4) by a local one. We shall now obtain local boundary conditions which are exact for waves that consist of the first M harmonics. For such waves, u,(R, 0) and u0(R, 0) can be expanded in the finite series M
M
u,(R,O)= ~.' (AmcosmO+BTsinmO),
uo(R,O)= Y/ (A'~cosmO+B'~sinmO).
m=O
(A1)
m~O
The kernels p~ in (25) can be written in the form
Z~icosn(O-O'),
i=j;
(~)
p~(0-0')= L~--R Z~j n sin n(O - 0'),
i Pj.
In (A2) the Z~ are functions of R. We now substitute (A1) and (A2) into the boundary condition (23). Imag u l 0.014 0.013 0.012 0.011 0.010 0.009 0.008 0.007
0.006 0.005 0.0042 0.003 0.002 0.001 0.000 -0.001 -0.002 -0.003 -0:004 -0.005
-0.006 I . . . . . .
129
' ' ' 1 ' '
130
. . . . . .
' 1 ' ' ~ ' ' ' ' ' ' 1
131
......
132
'''1
~'
.....
133
' ' 1 ' ' ' ' ' ' ' ' ' 1
134
. . . . . . .
135
¢ ' 1 ' ' '
136
. . . . . .
137
node number
Fig. 10. The plate-inclusion problem. The imaginary p a ~ o f the displacement u~ is shown along the a~iflcial boundary r = 3 in the fi~t quadrant.
276
D. Oivoli, J.B. Keller / Non-reflecting boundaries
The orthogonality of the trigonometric functions yields
M - Tr(R, 0) = ~,' [Z~'I(cos nOAh+sin nOB"r)+nZ~2(sin n O A h - c o s nOB~)],
n=0 M
-
(A3)
To(R, 0) = ~' [Z~2(cos nO A~ + sin nO B~) + nZ~'l(sin nO A~ - cos nO BT)].
n=0
Next we write Z~ as a finite sum of even powers of n, namely
M
n
Zij =
Z
m
m=0
Oto (n)2m;
n =0, 1 , . . . , M.
(A4)
For a fixed i and j, (A4) is a system of M + 1 linear algebraic equations for the coefficients ot~. More explicitly, this system has the form 1o
12
14
...
o M E M 4 ...
M2
ol b
Zb
La,"MJ
/z',T]
.
(A5)
Here 0° is defined to be 1. We now substitute the expansion (A4) into (A3). In doing so we make use of the following identities: (-1) m
(cos nO) = n 2m cos nO,
( - 1 ) m+l d02m+l d 2 ~ + ' (COS nO) = n 2m+' sin nO,
2~ (-1)"
(A6)
d2m+l (sinnO)=n2msinnO,
(--1)"~TX~--d~+l(sinnO)=n2"+~ cosnO. dlt
Then we obtain - T , ( R , 0) = Y.'
n=0 m = O
o~~ ( - 1)"
(A'~cosnO+B'~sinnO)
dEm+l
]
+ t~2(--1) m+l d02m+l (A~ cos nO+B~ sin nO) ,
M M[ - T o ( R , O ) = ~' Y. a ~ ( - 1 ) m n=O
(A7)
2~ (A~cosnO+B;sinnO)
m=O
d2m+l
]
+ a ~ ( - 1 ) m+~ d02m+~ (Am cos nO+B'~ sin nO) . But from (A1) this reduces to
d2mur
md2m÷%]
- T r = m=o ~ (--1)m a ~ 2dO ---~-a12 d--~J F ~ d2"Uo - T o = ,=0 ~ (-1)"L~22~-~21
'
.,
(A8) ~
J"
Equations (A8) are the desired local boundary condition. Since (A8) follows from the exact boundary condition, it is exact for waves which consist of the first M harmonics.
D. Givoli,J.B. Keller / Non-reflectingboundaries
277
We shall now derive explicit expressions for the first two boundary conditions in (A8). With M = 0, the system (A5) reduces to 0
O
(A9)
a o - Z u. From (A2), (25) and (18) we deduce E
Z°'-(l+~~)~Y
/' k H o L + l ~ L'~LOL "~'],
z%
=
z °, =
0,
(A10) E l Z°2 = (i ~-~v) (~kT ,--~,r+ I ~ Therefore, for M = 0 the boundary condition (A8) becomes E
(
HOL I \
E
- T , . = ( l + v ) TkL-'~OL+-~)Ur,
(! k HOT+I__~
- T o = ( l + v ) \ 2 THe)T RJUO.
(All)
This boundary condition is perfectly absorbing for axially symmetric cylindrical waves. For M = I, the system (AS) yields
o tou_- Z oo,
o ~l _u - Z u1 - Z ou.
(A12)
From (A8), the local boundary condition is d2ur - Tr = Z °, Ur -- ( Z l , - Z ° , ) - ~ +
t (Z,2
-
d3ur Z°2) dO3 , (A13)
-To=Z°2uo-(Z212-
z O ~ cl u e + : Z ~
d3uo
22, d02 t 2,-Z°,) dO3-
This condition is perfectly absorbing for axially symmetric waves as well as for those involving sin 0 and cos 0. When kLR and kTR are large, we can replace the Hankel functions by their asymptotic forms. For example, by using HoL/H~L~--i and HoT/H~T~--i in (All) wc get -T,=
(l+v)
u,,
TkLi-
-To=-(l+v-----j ~kTl----~ uo.
(A14)
To express the preceding boundary conditions in terms of Cartesian quantities we use vector transformations, as in Section 2. For example (All) and (A14) have the forms
-Tr(O)=fllur(O),
-To(O)=fl2uo(O).
(A15)
After transformation they become
- T~(O) = Aij(O) uj(O),
(A16)
where 2
Aq(O) = • J~k Vik(O) Vjk(O). k=l
The functions V0(0 ) are defined in (28).
(A17)
278
D. Givoli, J.B. Keller / Non-reflecting boundaries
For the sake o f completeness, we now use a similar procedure to localize the exact boundary condition derived in Keller and Givoli [8] for the reduced wave equation. This exact boundary condition is m " ( O - 0') u ( R , 0') R dO',
u~(R, O) = - Z'
(A18)
n=0
k H~I)'(kR) ~rg H ~ ) ( k g ) cos n ( O - 0').
m " ( O - 0') =
(A19)
Here u~ is the outward normal derivative of the unknown scattered field u, and k is the wave number. We write m " ( O - 0') =1_1_ Z "
~R
cos n(O - 0'),
H~)'(kR) Z " = - k H~1)(kR) .
(A20, A21)
Proceeding as above we obtain the local boundary condition M .
-m
~2 m U
--U~= ~ (--1) ama--~,
(A22)
m=0
where am is found from the solution of the system (A5), ignoring the indices ij. For M = 0 we have the boundary condition H~ol)'(kR) - u~ = - k H
(A23)
When k R is large we can use the asymptotic approximation H ' o / H o - i to obtain (A24)
- u~ = - iku.
This boundary condition has exactly the same form as the Sommerfeld radiation condition. In the case considered at the end of Section 3 where the original domain is finite, equations (A1)-(A13) hold with the Hankel function H~ ~) replaced everywhere by the Bessel function Jn. Thus, in the boundary condition ( A l l ) all the functions H are replaced by J. To obtain an asymptotic boundary condition analogous to (A14), we use the asymptotic approximations JoL/J~L~--cot(kLR--~/4 ) and JOT/J~T~ - c o t ( k T R -~r/4), which are valid for large values of kLR and kTR. Then we get the boundary condition T~-
-To
(l+v)
(, ootb
).
(1+,,)
We note that if a local boundary condition of the form (A16) is used instead of the exact boundary condition (4), the finite element formulation is exactly the same, except that the definition of b ( . , • ) becomes b(w, u) = J~ wiAijuj d ~ .
(A26)
Although (A26) has the same form as (39), the operator Aq in (A26) is local whereas M~j in (39) is not.
Acknowledgments This work was supported in part by ONR, AFOSR and NSF. The first author was also supported by the Technion's postdoctoral grant no. 27000827 and by the L. Kraus Research Fund.
D. Givoli, J.B. Keller / Non-reflecting boundaries
279
References [1] B. Engquist and A. Majda, "Absorbing boundary conditions for the numerical simulation of waves," Math. Comp. 31, 629-652 (1977). [2] A. Bayliss and E. Turkel, "Radiation boundary conditions for wave-like equations," Comm. Pure Appl. Math. 33, 707-725 (1980). [3] K. Feng, "Finite element method and natural boundary reduction," Proceedings of the International Congress of Mathematicians, Warsaw, 1439-1453 (1983). [4] L. Ting and M.J. Miksis, "Exact boundary conditions for scattering problems," J. Acoust. Soc. Am. 80, 1825-1827 (1986). [5] R.L. Higdon, "Numerical absorbing boundary conditions for the wave equation," Math. Comp. 49, 65-90 (1987). [6] R.C. MacCamy and S.P. Matin, "A finite element method for exterior interface problems," Int. J. Math. Math. Sci. 3, 311-350 (1980). [7] C.I. Goldstein, "A finite element method for solving Helmholtz type equations in wave guides and other unbounded domains," Math. Comp. 39, 309-324 (1982). [8] J.B. Keller and D. Givoli, "Exact non-reflecting boundary conditions," J. Comp. Phys. 82, 172-192 (1989). [9] J. Lysmer and R.L. Kuhlemeyer, "Finite dynamic model for infinite media," J. Eng. Mech. Div. ASCE 95, 859-877 (1969). [10] W.D. Smith, "A nonreflecting plane boundary for wave propagation problems," J. Comp. Phys. 15, 492-503 (1974). [11] C. Cerjan, D. Kosloff, R. Kosloff and M. Reshef, "A nonreflecting boundary condition for discrete acoustic and elastic wave equations," Geophysics 50, 705-708 (1985). [12] J. Sochacki, R. Kubicheck, J. George, W.R. Fletcher and S. Smithson, "Absorbing boundary conditions and surface waves," Geophysics 52, 60-71 (1987). [13] B. Engquist and A. Majda, "Radiation boundary conditions for acoustic and elastic wave calculations," Comm. Pure App. Math. 32, 313-357 (1979). [14] R.W. Clayton and B. Engquist, "Absorbing boundary conditions for acoustic and elastic wave equations," Bull. Seismol. Soc. Amer. 67, 1529-1540 (1977). [15] C.L. Scandrett, G. A. Kriegsmann and J.D. Achenbach, "Application of the limiting amplitude principle to elastodynamic scattering problems," S I A M J. Sci. Star. Comp. 7, 571-590 (1986). [16] R.L. Higdon, "Absorbing boundary conditions for elastic waves," SIAM J. Num. Anal., to appear. [17] D. Givoli and J.B. Keller, "A finite element method for large domains," Comp. Math. Appl. Mech. Eng. 76, 41-66 (1989).
261
N O N - R E F L E C T I N G B O U N D A R Y C O N D I T I O N S FOR ELASTIC WAVES Dan GIVOLI Department of Aeronautical Engineering, Technion, Haifa 32000, Israel
Joseph B. K E L L E R Departments of Mathematics and Mechanical Engineering, Stanford University, Stanford, CA 94305, USA Received 2 March 1989, Revised 16 October 1989
An exact non-reflecting boundary condition is devised for time-harmonic two-dimensional elastodynamics in infinite domains. The domain is m a d e finite by the introduction of a circular artificial boundary on which this exact condition is imposed. In the finite computational domain a finite element method is employed. Numerical examples are presented in which the accuracy a n d efficiency of the method using the exact non-local boundary condition are compared with those of methods based on approximate local boundary conditions. The method is also used to solve problems in large finite domains by reducing them to smaller domains. In addition, local boundary conditions are derived which are exact for waves with a limited n u m b e r of angular Fourier components.
1. Introduction To solve a wave problem numerically in an unbounded domain one usually introduces an artificial boundary 90, thereby defining a finite computational domain O. On 90 one must impose some boundary condition which should allow waves f r o m / 2 to pass through 90 without generating spurious reflections. We shall find an exact boundary condition having this property, for problems involving two-dimensional time-harmonic elastic waves. In O we shall employ the finite element method. For problems of wave propagation there has been a considerable amount of previous work to devise appropriate artificial boundary conditions. The time dependent scalar wave equation and the reduced wave equation were considered by Engquist and Majda [1], Bayliss and Turkel [2], Feng [3], Ting and Miksis [4], Higdon [5], MacCamy and Marin [6], Goldstein [7] and Keller and Givoli [8]. In references [2] and [5] approximate local boundary conditions are derived. In [1] and [3], exact nonlocal conditions are derived and approximated by local ones, whereas in [4], [6], [7] and [8] exact non-local boundary conditions are employed. To get an exact boundary condition for elastic waves we use the procedure presented in [8] for the simpler case of the reduced wave equation. Artificial boundary conditions for elastic waves have been considered mainly in the context of geophysics. In two dimensions the artificial boundary 90 is traditionally composed of straight line segments and the computational domain O is usually a rectangle. Often the original domain is the whole space, possibly with some holes in it, which corresponds to wave propagation deep inside an elastic body. We now describe some of the previous methods. We denote the coordinates normal and tangential to 90 in two dimensions by Xl and x2 respectively. We also let Ul and u2 denote the displacements and T1 and T2 the tractions in these directions. 0165-8641/90/$03.50 @ 1990--Elsevier Science Publishers B.V. (North-Holland)
262
D. Givoli, ZB. Keller / Non-reflecting boundaries
Lysmer and Kuhlemeyer [9] suggested the following boundary condition on ~ : Oul =
apcL O----t T,,
bpCT~t= T 2 .
(1)
Here p is the mass density, t is time, CL is the longitudinal wave speed, CT is the transverse wave speed, and a and b are dimensionless parameters. The parameters a and b were chosen to minimize the reflected energy for an incident plane wave which hits the boundary ~ at a given angle of incidence. The choice of a and b was considered separately for an incident longitudinal wave, for an incident transverse wave, and for a surface wave. It was suggested that a = b = 1 is a good choice in general. In the computational domain ~ a finite element scheme was used. Smith [10] showed how the reflections from a piece-wise planar boundary ~ can be eliminated by adding together the numerical solutions of several problems, in each of which a certain combination of Dirichlet and Neumann boundary conditions is used on ~. If reflections are to be eliminated exactly on n plane surfaces, 2 n such solutions must be added together. This procedure is applied separately to dilatational waves and to surface waves. Cerjan et al. [11] and Sochacki et al. [12] each presented a method in which the amplitudes of ul and u2 are gradually reduced in a strip of nodes adjacent to ~. Thus, the solution is artificially damped in the vicinity of ~. On ~ , the usual Dirichlet or Neumann boundary condition is used. Engquist and Majda [13] found an exact non-local boundary condition on ~ in terms of pseudodifferential operators, and then approximated it by local boundary conditions. The simplest of these approximate conditions is -
CL
Ul = 0,
--
CT
U2 =
0.
(2)
This condition is perfectly absorbing for plane waves at normal incidence. Their next boundary condition involves a linear combination of the operators 02/0x, at, 02/0t 2, 02/aX2 0t and a2/ox 2. These boundary conditions are used in a difference scheme employed in a rectangular domain. A special procedure is suggested to avoid instabilities near the corners. Some numerical experiments with these boundary conditions are presented in Clayton and Engquist [14]. Scandrett et al. [15] obtain time-harmonic solutions of the equations of elastodynamics by solving numerically the time-dependent problem and using the limiting amplitude principle. They derive an approximate time-dependent boundary condition on ~ , which is equivalent to the first-order Engquist and Majda condition. An analogous boundary condition for the three-dimensional case is also obtained. A finite difference scheme in time and space is employed. Higdon [16] suggested using a difference scheme in ~ with a boundary condition on ~ of the form
Here the fl~ are parameters which can he adjusted so that the boundary condition is perfectly absorbing for an incident plane wave which hits ~ at a given angle. In the case m = 1 with normal incidence, (3) reduces to the Engquist-Majda condition (2). If m is large enough, the resulting high-order boundary condition can be made perfectly absorbing for a number of chosen incidence angles and for both longitudinal and transverse waves. Surface waves were not treated. The boundary conditions (1), (2) and (3) can be used in solving either the time-dependent elasticity equations, or those for time harmonic waves with frequency ~o. In the latter case the equations involve
D. Givoli,J.B. Keller / Non-reflectingboundaries
263
the space variables alone. Then (1)-(3) can still be used with auj/~t replaced by itouj. Since they are approximate, they generate spurious reflection from ~. The magnitude of the reflected wave depends upon the specific problem. It is hard to estimate because of the presence of various types of waves. We shall derive an exact non-local boundary condition on ~ which will eliminate all reflections for any type of wave. It has the form - Ti(x) = Mijuj(x)
m .=0 ~"
f~ m~(x, x') us(x') dx'.
(4)
The summation convention for repeated indices is used here and elsewhere. This boundary condition is analogous to the condition introduced in [8] for the reduced wave equation. In Section 2 we derive the boundary condition (4) for time-harmonic waves in two dimensions, when the artificial boundary ~ is a circle. We consider problems of plane stain and problems of plane stress. In Section 3 we describe the finite element scheme employed in/2. A similar formulation is discussed in more detail in [8] in the context of the reduced wave equation. We also show how the finite element method with the boundary condition (4) can be used to solve efficiently problems in large finite domains by eliminating an interior circular region. Numerical experiments are presented in Section4. The results obtained by using the exact boundary condition are compared with those obtained by using the approximate Lysmer-Kuhlemeyer boundary condition (1) and by using some of the approximate boundary conditions obtained in the appendix. We compare both the accuracy and the efficiency of the various schemes. In the Appendix we localize the exact boundary condition and obtain a sequence of local approximate boundary conditions. The M t h boundary condition in that sequence is exact for waves that consist of the first M Fourier components in the angle variable 0. We also obtain analogous local boundary conditions for the reduced wave equation.
2. An exact non-reflecting boundary condition We consider the equations of linear elastodynamics in an infinite two dimensional region. To make the region finite, we introduce an artificial boundary ~ which is a circle of radius R. Our goal is to derive an exact boundary condition of the form (4) on ~. We shall do so by finding an exact solution in the region exterior to ~. In that region the medium is assumed to be linear, homogeneous, isotropic and free of body forces. Then outside ~ the elastodynamic equations in Cartesian coordinates are
p (3, =~ u,.jj+(,+~) u.,.
(5)
Here U~ is the displacement in the i direction, i = 1 or i = 2, p is the mass density, and a comma indicates partial differentiation. The constants /~ and A are related to Young's modulus E and to Poisson's ratio v via
E ft = 2(1+ v) ;
f rE~(1+ v)(1 A = /[ v E / ( l _ v 2 ) '
-
2v),
plane strain, plane stress.
(6)
To treat time-harmonic solutions we set
Uj(x, t) = uj(x) e i~'',
(7)
substitute (7) into (5) and obtain p. u~,~+ (A +/z) ujd~+~o2put =0.
(8)
D. Givoli, J.B. Keller / Non-reflecting boundaries
264
We now use the Helmholtz decomposition theorem to write ui in terms of two potentials 4) and 0: u l = 4<1 + 0,2,
(9)
u2 = 4<2- 0,1.
Then (8) will be satisfied if 4) and 0 satisfy V2~b+ k2~b = 0,
V20 + k20 = 0.
(10)
Here k L and kT are the longitudinal and transverse wave numbers: k2=to 2
P A +2/.~'
k2=to 2 p . /~
(11)
At infinity the Sommerfeld radiation conditions must hold:
!ira rl/2(C~,,--ikLCk)=0,
!im rl/2(0.,--ikT0) =0.
(12)
The general solutions of (10) and (12) are, with a. and b. arbitrary constants,
cb(r, O)= ~ a.H~l)(kL r) e i.o ,
O(r, O) = ~ b.Ht.')(kTr)
n=--oo
e i"°.
Here H(,]) is the Hankel function of the first kind. Since the displacements u~ and polar coordinate form of (9),
U,.=qb.r+lo.o,
(13)
n=--oo
uo are
Uo=--O.,.+I~,o,
given by the
(14)
equations (13) and (14) yield
Ur= ~ [a.kLH(1),(kLr)+b inH(,)(kTr)] e i.o, n = --oo
r
(15)
Now we express the constants a. and b. in terms of the values u,(R, 0) and uo(R, boundary r -- R. To this end we first expand u,(R, 0) and uo(R, O) in Fourier series:
[1I?
=
e -i~°' u.(R, 0 ' ) d O '
]
O) on
the artificial
ei'S,
n =-oo
(16) u o ( R , O) = . . ._~.
Then
e -i"°'
us(R, 0') dO'
e i'°.
we evaluate (15) at r = R and compare it with (16) to obtain
a. -- 2~rA. ~ H.TUo(R) 1 I o " e-i"°' [ - k T H n T' U r ( R ) O ' ) - -in
0') ] dO',
(17) 2erA.
Here
, Uo(R. 0') - - ~ e -i.o' kLHnL
HnL u,.(R, 0')
dO'.
265
D. Givoli,J.B. Keller / Non-reflectingboundaries An
~
n2 I f --kLk'rHnLHnT -I-_R_ 2 HnLHnT
(18)
and
H~L = H(.'(kLR),
H.T
dH("l)(x)
HrnL
=
H'-r
dx
H(.'(kTR),
(19)
dH("l)(x)x=k~a'
dx
X=kL R '
Substituting (17) into (15) gives
-uAr, O)= ~
[sTdr, O,O')u,(R,O')+sT~(r,O,O')uo(R,O')]dO',
n=--oo
-uo(r,O)=
I
(20)
~
[s~l(r,O,O')u~(R,O')+s~2(r,O,O')uo(R,O')]dO',
n=-oo
where
r
]
in<°-°')L kLkrH,TH~l),(kLr)_._g_ s'il(r, O, 0 ' ) - e2~a-------~ ~r H~LI-I~l)(krr) , n
_
_
1?!
t
(1)
kLHnTHO)'(kLr)----r kLHnLHn (krr) ,
$12(r, 0, 0') -- 2~rA.
(21) [ ~ kTH.TH~ , ~)(krr)-'-~ in kTH.LH~),(kTr)], - 2,rra,, sEt(r, O, O') - ein(#-°')
kLkTHnLH" (kTr) "
s~2(r, 0, 0 ' ) = 2,rrA.
Our goal is to express the tractions T,(R, O) and To(R, O) in terms of u,(R, O) and uo(R, 0). We first consider plane stress, for which
(
v
v )
T~ - (1 - v 2) u,,, +rUO.O+rU, ,
E (lu,,o+uo._luo)
To = 2(1+ v)
.
(22)
We use (22) to calculate the tractions from (20), set r = R, and after some algebra we get
- T,(R, O) = E'
[ M d o , 09 ur(R, o')+pl%(o, 09 uo(R, o')]R dO',
rl=O
(23)
-To(R,O)= E'
[p~l(O,O')u,(R,O')+p2~(O,O')uo(R,O')]RdO'.
n=O
The prime on the sum indicates that a factor of 1/2 multiplies the term with n = 0. We eliminate H~L and H~-r from the kernels p~ by using the identity
H~l)"(x)=(~2-1) H~l)(x)-lH~l)'(x)'x
(24)
266
D. Givoli, J.B. Keller I Non-reflecting boundaries
Then we finally obtain p T , ( o - o ' ) =( l-+-Eu-) cos ,'itn(O RA,,
,
pT~(o - o') = (1 + v)
o,)r|
2 , t. 3 " k L k ' r H ' L H ' r +
~g~,~.
k L k v H,LH,T--~' ' R
n2 HnLHnT],
k~k,~H'~H',;
]. (25)
p~,( O - 0') = -p'~2( O - 0 %
E c o s n ( 0 - 0 ' ) [ l l " 1.2u, u ±1 n2 ] (1 + u-----) ~-~-~ L~'~L'~V'"nLt'nTT'RkLkTH~nLH~nT----~ H n L H n T _ ,
p~2(O - 0')
3' = 1/(1 -- l')
(plane stress).
(26)
Equations (23) and (25) constitute the exact boundary condition on ~ in polar coordinates. To obtain the boundary condition in rectangular coordinates, as in (4), we transform T,, To, u,, uo in (23) into T~, T2, ul, u2. This leads to the following relation between the polar kernels p~ and the rectangular kernels rots: (27)
m~(O - 0') =p~a(O - 0') V~g(O) Vj,(O').
Here Vl,(X) = V22(x) = cos x,
V~2(x)= - s i n x,
V21(x) = sin x.
(28)
It is easy to see that the kernels m ~ ( O - 0') possess the symmetry (29)
m ~ ( O - 0') = mj~( O ' - 0).
This property is essential to preserve symmetry in the finite element scheme to be employed in the computational domain/2, as we shall see in Section 3. Equations (25) and (26) were obtained for plane stress. Now we replace (22) by the corresponding relations for plane strain. Then we find that the formulae (25) for p~ are the same, with the definition (26) of 3' replaced by 3' =
l-z,
(plane strain).
1-2v
(30)
3. Finite element formulation
We shall now combine the boundary condition (4) with the finite element method in the computational domain O. Since this formulation is analogous to that presented in [8] for the reduced wave equation, it will be described briefly. The computational domain ~2 is bounded in part by the circular artificial boundary ~, and in part by some interior boundary F = F s w Fh. T h e problem in 12 is: ( CoklUk.r). j + tO2pUi = 0 ui = gi
in O,
(32)
on Fs,
Ti = CijklUk.tnj = hi - Ti = Mijuj
on ~.
(31)
on Fh,
(33) (34)
D. Givoli, ZB. Keller / Non-reflecting boundaries
267
Equation (31) is the generalization of equation (8) to a non-isotropic and non-homogeneous medium, which may be present in O. In (33), nj is the component of the unit outward normal to Fh. Equation (34) is the boundary condition (4). To obtain a weak form of the problem (31)-(34), we introduce the two sets
S¢= { u l u e n l ( 1 2 ) x n l ( o )
& u,= g, on Fg},
6eo={wlweHl(O)xH'(O)
& w , = 0 on rs}.
(35) Here H~(O) is the Sobolev space containing functions which are square-integrable and which have square-integrable first derivatives. We also introduce the notation
v(i,j) I- 3( vi,j + vj,i).
(36)
Then the weak form of the problem is: Find u ~ 5e such that for all w ~ Seo
(37)
a(w, u) + b(w, u) = (w, h)r where
f.
w,h,
,..
,39.40,
The finite element method is obtained by approximating the weak form (37). The domain D is discretized into a finite number of elements, and each element is associated with a finite number of nodes. Then u and w are approximated by the sums u h and wh: uh(x) = ~ dA,NA(X)+ E gA,NA(X)',
gA,=g,(XA),
(41)
AE~g
A~71
w~(x)= ~ CAiNA(X).
(42)
A~rl
Here rig is the set of nodes on the boundary Fg, 7/is the set of all other nodes, NA is the shape function associated with node A, :cA is the location of node A, and CAi and dai are constants. We now substitute (41) and (42) into (37). This results in the finite element matrix form of the problem for the vector d:
Kd=F,
K=K°+K
K a = [K~4iBj],
b,
(43,44)
K b = [K~,aj],
K~AiBj = a( NAei, NBej),
d = {dsj},
K bA~Bi= b( NAe~, N ~ ),
Far = (Naei, h)r -~, Y. gBia(Nae~, Na~).
F = {FA,},
(45) (46) (47)
J BE~ B
Here e~ is the unit vector in the x~-direction. On the left side o f (46), Ai and Bj are the positions in the matrix corresponding to degree of freedom i at node A and degree o f freedom j at node B. A similar meaning applies to the Ai on the left side of (47). Equation (44) shows that the effect of the boundary condition on the standard finite element scheme is just the inclusion o f K b in the matrix K in (43). The entries of the matrix K b can be written explicitly
D. GitJoli, J.B. Keller / Non-reflecting boundaries
268 as
b
K A~aj= b( Nae,, NBe1) = = ,=o ~
L NAM~jNBd~3
f~ f~ NA(X)m~(x,x')Ns(x')dx' dx;
To evaluate the double integral in (48) we note that 4
(48)
m~(O-0') given by (27) can be written in the form
4
m~(O- 0') = E E a~k' F~,(O) F'~(O') k=l
(49)
I =1
where
F~(O)=cosnOcosO, F~(O) = sin nO cos O,
F~(O)=cosnOsinO, F~(O)= sin nO sin O,
(50)
and the a~ kl are constants. For example, all,,1 = a~,33 is the coefficient multiplying cos equation of (25). Therefore (48) gives .=o k=, ,:, a,j ~ j ~
NA(x) Fg(x) dx
N.(x) FT(x) dx .
n(O- 0') in the first
(51)
b Thus, only one-dimensional integrals have to be computed to form K b. Of course, KA~Sj is nonzero only b if both nodes A and B lie on the boundary ~ . Also, to compute KA~nj in practice, the infinite sum in (51) has to be truncated after a finite number of terms. Other computational aspects of the method are discussed in [8]. The matrix K b is symmetric. This follows from (48) and from the symmetry of the kernels m~ in (29). Thus the exact boundary condition preserves the symmetry of the standard finite element formulation. After the finite element solution in g2 has been found, the solution outside g2 can be determined from equations (20). The finite element solution Uhr(R,0') and uho(R, 0') can be used to replace u,(R, 0') and uo(R, 0') in those equations. The present finite element method can also be used to solve problems in large finite domains. To see this, we consider a finite two-dimensional domain ~ with an irregularly shaped outer boundary 0~. We suppose that the medium is homogeneous and isotropic, except perhaps in a region near a~. To obtain a small computational domain g~ we cut out of ~ a circular domain D within which the medium is homogeneous and isotropic. The boundary of D is denoted ~. Thus the computational domain ~ = ~ - D is b o u n d e d internally by ~ and externally by aR. By finding an analytic solution u in D, we derive an exact boundary condition on ~ of the form (4). Then we combine this boundary condition with the finite element method in £2. The finite element formulation is identical to the one given before. To get the exact boundary condition we go through a procedure analogous to the one outlined in Section 2. Then we find that equations (8)-(30) hold with the Hankel function H~1) replaced everywhere by the Bessel function .In. Since the boundary condition is real in this case, the whole problem is real.
4. Numerical experiments We now compare the results obtained by using the finite element scheme employing the exact boundary condition with results obtained by using the scheme employing local approximate boundary conditions
D. Givoli, J.B. Keller / Non-reflecting boundaries
269
and with the exact solution. We consider first the case where the domain is infinite. The local boundary conditions that we use are the one by Lysmer and Kuhlemeyer (1) with the parameters a,= b = 1, and two local conditions obtained in the appendix, namely the condition (A11) and its asymptotic counterpart (A14). The local boundary conditions (2) and (3) are not compatible with the finite element method, because they involve the N e u m a n n operator Ou~/On rather than the traction T~. We consider an infinite plate with a circular hole of radius a. Now from (15) it is easy to see that the following are two solutions of the elastodynamic equations:
Solution 1 (longitudinal mode): Ur(r, O) = k L H ~ ( kLr) c o s nO,
uo(r, O) =
n_ Hn(kLr) sin nO. r
(52)
Solution 2 (transverse mode): n
u,(r, O) = - - H , ( k T r ) sin nO, r
uo(r, O) = - k T H ' ( k T r ) cos nO.
(53)
In order to have (52) or (53) be the exact solution for the plate with a hole, we use as boundary condition on the hole boundary that obtained by setting r = a in (52) or (53). We choose a = 0.5 and we set the artificial b o u n d a r y ~ at r = R = 1. The other parameters are to = 1, p = 1, E = 1 and ~, = 0.3. The mesh is composed of 120 bilinear quadrilateral elements. It is shown in Figure 1. In the exact boundary condition (23) we use the first five terms. In fact, only the term corresponding to n in (52) or (53) contributes to the solution.
Fig. 1. Finite element mesh for the infinite plate with the circular hole.
We first solve the axisymmetric problem with the longitudinal mode. This case corresponds to n = 0 in equations (52). Figure 2 shows the real part o f the displacement u~ on ~. Five solutions are compared: the exact solution, the one obtained with the exact boundary condition, and the three solutions obtained by using the local boundary conditions (1), ( A l l ) and (A14). We see that the exact boundary condition and its local version ( A l l ) give results which coincide with the exact solution. The local condition (A11) is indeed expected to give identical results to those obtained with the non-local exact b o u n d a r y condition because ( A l l ) is exact in the axisymmetric case. The solution with the asymptotic condition (A14) is about 8% off the exact solution, whereas the L y s m e r - K u h l e m e y e r solution is almost 50% off at the peaks. Our next problem corresponds to n = 1 in equations (52), still in the longitudinal mode. Figures 3 and 4 present the real parts of u I and u2 respectively along ~ . In both figures the results obtained with the
D. Oivoli, J.B. Keller / Non-reflecting boundaries
270 Real
. ., . , ,. . ,.
ul
0.6"
...... y.' .....................
0.5" 0.4" 0.3 ~ 0.2 ~ 0.1"
...////'//'"
""''".....,,%.
0.0 : -0.1 -0.2 -0.3 -0.4 -0.5 -0.6
~
.....
-0.7 -0.8 -0.9 -1.0 '
.
.
.
.
.
.
.
'
120
I
'
'
'
'
.
.
.
.
'
I
130 node
LEGEND
'
'
'
'
'
'
'
'
'
140
.
.
.
.
.
.
.
.
.
160
number
Lys.-Kuh. exact b.c. local N=O
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
I
150
- ...............
asymp, b.c exact sol.
Fig. 2. The infinite plate problem: comparison of results obtained with various boundary conditions for longitudinal axisymmetric waves. The real part of the displacement u I is shown along the artificial boundary r = 1.
Reol
ul 2.2'
1.8 1.6
! 1.4
i
1.2
/
1.0
0.6
i
0.4
x ",
/
t/
,
\
/
\
/,
I
i
i
i
///"
//.
/
,, \
',.
"',. \
•
".
\\
.
.
.
.
.
.
.
.
.
.
///"
/ .,,"
., x
",,\ \
.
\
-...x.--.;i I
.
.
.
.
.
.
.
.
.
130
I
. . . . . . . . .
140 node
LEGEND
i
/ , ///
",,
"-.'2~- .>.." 120
/// "'\'k /" '\ .....'""...,.
/..,.'"'"'"',,,,, ;/ "...
2.0
.................................. L y s . - K u h . ........... exact, b.c. . . . . . . . . local 14=0
I
.
.
.
.
.
.
.
.
.
150
160
number - ...............
asymp, b.c exact sol.
Fig. 3. The infinite plate problem: comparison of results obtained with various boundary conditions for longitudinal 'cos 0 waves'. The real part of the displacement ul is shown along the artificial boundary r = 1.
D. Givoli, J.B. Keller / Non-reflecting boundaries Real
271
u2 1.0
.:..-'""'.,..
0.9 0.8
.......''"".
/."
/ /
)\_
0.7
0.6
/-
\
0.5
0.4 0.3 0.2 0.1 O.O -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7
\
/
/
i :
!"
x
~.,~
-~. -
"x
/ ~.~ ":/"
"~
.i/
"~:X
-\
./
'..
\~,~
/-/
/
\-\
/
.,, .~ /
x\
y/:
'"",.,,,,,
/." t.\.
". /.":: "~........."
".........f
-0,9 ~1.0 . . . .
120
'
'
'
'
'
1
'
'
'
.
.
.
.
.
.
130
.................................. L y s . - K u h . . . . . . . . . . . . e x a c t b.c.
........
I
,
'
'
,
. . . .
140 node
LEGEND
,,
/
/
\
-0.8
.
j':"
local
'
1
'
'
,
,
. . . .
150
,
160
number - ...............
asymp,
exact
b.c sol.
M=O
Fig. 4. The infinite plate problem: comparison of results obtained with various boundary conditions for longitudinal 'cos 0 waves', The real part of the displacement u 2 is shown along the artificial boundary r = 1.
exact b o u n d a r y condition deviate from the exact solution only slightly. On the other hand, the results using the three local approximate boundary conditions deviate considerably. The condition (1) gives unacceptable results. It is a coincidence that the local condition ( A l l ) is even worse than its asymptotic counterpart (A14); it so happens that the error due to the asymptotic approximation cancels part of the finite element error. N o w we solve the problem which corresponds to n = 1 in equations (53), in the transverse mode. In Figure 5 the real part of u2 on ~ is shown. The exact boundary condition again gives excellent results whereas the approximate conditions do not. We note that the results obtained from the conditions (A11) and (A14) tend to err more symmetrically with respect to the exact solution than condition (1). This is apparent both in Figure 3 and Figure 5. Next we move to a problem in a finite domain. We consider time-harmonic waves in a thin circular plate of radius a. The displacements on r = a are prescribed. We delete a circular domain of radius R. Then we impose the exact boundary condition on r = R. The computational domain is the ring R < r < a. Equations (52) and (53) are still solutions of the ¢lastodynamic equations if the Hankel functions H are replaced by the Bess¢l functions J. On r = a we impose the boundary condition obtained by setting r = a in (52) or (53) (with H replaced by J). We choose a = 1 and R = 0.5. The other parameters are the same as for the infinite domain problem. The mesh shown in Figure 1 is used again, but this time the non-reflecting boundary condition is imposed on the inner circle. We compare the exact solution with that obtained by using the exact boundary condition (A11) (with H replaced by J) and the local asymptotic condition (A25). Figure 6 depicts ul along the artificial b o u n d a r y r = 0.5 in the axisymmetric longitudinal mode. As in the infinite case, solutions with the exact boundary condition and the local condition (A11) coincide with the exact solution. The asymptotic condition gives bad results. In fact, the argument kLR of the Bessel
D. Givoli, J.B. Keller / Non-reflecting boundaries
272
Ree] u2 1.2 1.0
......'"""'",,...,,
0.8
.\
0.6
\
/" .;" ..."
/
/,..'"'"'",,,,.., ~\\ "..\
\
i
0.4
\
/
0.2
i / ":..... " • ",
"--'-:~" \ '"",.
0.0 -0.2
-"-
-0.4
,,.
~ " - ~ x " .x.- ~
i .,/'
x
,'"
~,.
i.~ , / , :
i ,.',"
~ . \ . ~ . . .. //'/z
-0.6 .....f,-
-0.8
".......,."
-1.0 .
.
.
.
.
'
'
'
'
120
I
'
'
'
'
'
.
.
.
.
I
130
.
.
.
.
.
.
.
.
.
I
140
.
.
.
.
.
.
'
1.50
'
'
160
node number ..................................Lys.-Kuh. ........... exact b.c. ........ l o c a l N-O
LEGENO
- ...............
asymp, b.c exact sol,
Fig. 5. The infinite plate problem: comparison of results obtained with various boundary conditions for transverse 'cos O waves'. The real part of the displacement u2 is shown along the artificial boundary r = 1.
ul 0.7 0.6
............................... .....",..
0.5 0.4
...""...... "",., "
0.3 0.1 0.0 -0.1 -0.2 -0.3
~
//"
~
/
"",.
.:/
""',.,",,.
-0.4
../""
...........................................
-0.6-0"P -0.7 -0.8
-0.92 -1.0 2
. . . . . . . . .
120
I
. . . . . . . . .
130
I
. . . . . . . . .
140
I
.
.
.
.
.
.
.
.
.
150
160
node number LEGEND
.................................. esymp, b.c exact s o l .
exact b.c. l o c a l H=O
Fig.&Thecircu~arp~r~b~m:c~m~ari~n~£resu~ts~btainedwithva~usb~undaryc~nditi~ns~r~n~tudina~axisymm~t~c waves. T h e d i s p l a c e m e n t u t i s s h o w n a i o n g t h e a ~ i f i d a l b o u n d a r y r = 0 . 5 .
D. Givoli, J.B. Keller I Non-reflecting boundaries
273
functions has the value 0.48 here, and for this value the asymptotic approximation invoked in (A25) is very bad. On the other hand, in the previous problem the argument was twice as large and, more important, the asymptotic approximation of the Hankel functions is quite good even for such small arguments. In Figure 7 the displacement u2 is shown in the case n = 1 of the longitudinal mode. The exact boundary condition is again superior by far to the approximate ones. We have measured the CPU times needed to obtain the results above with the various boundary conditions. After averaging the times used in all the cases considered, we estimate the running time in Table 1. The time corresponding to the Lysmer-Kuhlemeyer condition (1) is normalized to 1. We see that the computation using 5 terms of the exact boundary condition takes about 10% more time than that using the local conditions. This is a small price to pay for the tremendous gain in accuracy obtained in the preceding examples. Finally, to demonstrate a possible application of the method when an explicit analytic solution cannot be obtained, we consider the following problem. An infinite plate with a circular hole is given. Young's modulus of the plate material is Ep = 10. Into the circular hole an inclusion with a square hole is fitted, which is of a material ten times 'softer' than that of the plate, i.e. with Young's modulus of E I : 1. Periodic expansion and contraction of the square hole produces elastic waves in the inclusion and in the plate. u2 0.26" 0.21
/
0.16 0.11 .,i
"\ \
0.06
\ ,, / . . . .
.\
i,/
',.\
i
\",,
,/"
0.01
.,,'/
"\
/
"..'\
/ :"
-0.09
/
,, \
]'
/'~ , , / /
""'k
-0.14 -0.19 -0.24
//
'-'.,
/] -0.04
•
.........
120
k'.,X,,, / :/" ',.,,.,..,,./ I .........
I .........
130
LEGEND
..................................
asymp,
execL
"
'k,,,., // ",.,.,...,../' I .........
140 node
/ ' " ~--
150
160
number
b.e sol.
................ "..........
exact local
b.c. N=O
Fig. 7. The circular plate problem:comparison ofresults obtained with various boundary conditionsforlongitudinal'cos 0 waves'. The displacement u2isshown aiongthe a~ificialboundaryr=0.5. Table 1 Normalized CPU times for the finite element schemes employing various non-reflecting boundary conditions B.C.
Lys.-Kuh. B.C. (1)
Asymptotic B.C. (A14)
Local M = 0 B.C. (All)
Exact B.C. 2 terms
Exact B.C. 5 terms
Normalized time
1
1
1.02
1.07
1.11
D. Givoli, £B. Keller / Non-reflecting boundaries
274
The configuration is shown in Figure 8. The parameters are: v = 0.3 and p = 1 for both materials, to = 10, the radius of the inclusion is 2.5 and the edge of the square hole is 2. On the square boundary the normal displacement is prescribed to have the value 1, whereas the tangential displacment is unspecified. The artificial boundary ~ is set at r = R = 3. Thus of the plate material only a ring of thickness 0.5 is included in the computational domain. We solve this problem using the exact boundary condition on with 5 terms in the expansion (23)• The mesh includes 128 bilinear quadrilateral elements: the inclusion comprises three 'rings' of 32 elements each, and the domain of the plate is composed of an additional ring of 32 elements. Figure 9 shows the result for the real part of u~ along the radial line y = 0, starting from r = 1 at the middle of the side of the square hole and ending at r = 3 on the artificial boundary. The point r = 2.5 is •
.
°
.
.
.
.
•
.......
• " " : lZli kl'O:
: i i ! ! i i i i-:~'~": . . . . . . .
~: i ! i i i i
f..
.
. % . . . . . . .
::::::/: . . . . .
:,~:::::: J..
.
:::::': :::::r:
i iii i~:
.~. . . . . . .
:i::::: :l:::::
i,'iiiii
"::::::.~:..
.:/::::::"
Fig. 8. The set up of the plate-inclusion problem.
Real
ul 1
-1
-2 '
'
'
'
'
'
'
'
'
I
'
'
'
'
'
'
'
'
'
2 r
Fig. 9. The plate-inclusion problem. The real part of the displacement u~ is shown along the radial line y = 0.
D. Givoli, J.B. Keller / Non-reflecting boundaries
275
at the plate-inclusion interface. The displacement u2 is zero on that line due to symmetry. In Figure 10 the imaginary part of ul is shown along the artificial boundary in the first quadrant. Due to symmetry the distribution is similar in the other quadrants. In fact, the symmetry can be exploited by discretizing only a quarter of the plate and by using symmetry boundary conditions. If that is done the exact non-reflecting boundary condition must be modified according to the guidelines given in [17].
Appendix
Approximate local boundary conditions Sometimes it may be desirable to approximate the nonlocal boundary condition (4) by a local one. We shall now obtain local boundary conditions which are exact for waves that consist of the first M harmonics. For such waves, u,(R, 0) and u0(R, 0) can be expanded in the finite series M
M
u,(R,O)= ~.' (AmcosmO+BTsinmO),
uo(R,O)= Y/ (A'~cosmO+B'~sinmO).
m=O
(A1)
m~O
The kernels p~ in (25) can be written in the form
Z~icosn(O-O'),
i=j;
(~)
p~(0-0')= L~--R Z~j n sin n(O - 0'),
i Pj.
In (A2) the Z~ are functions of R. We now substitute (A1) and (A2) into the boundary condition (23). Imag u l 0.014 0.013 0.012 0.011 0.010 0.009 0.008 0.007
0.006 0.005 0.0042 0.003 0.002 0.001 0.000 -0.001 -0.002 -0.003 -0:004 -0.005
-0.006 I . . . . . .
129
' ' ' 1 ' '
130
. . . . . .
' 1 ' ' ~ ' ' ' ' ' ' 1
131
......
132
'''1
~'
.....
133
' ' 1 ' ' ' ' ' ' ' ' ' 1
134
. . . . . . .
135
¢ ' 1 ' ' '
136
. . . . . .
137
node number
Fig. 10. The plate-inclusion problem. The imaginary p a ~ o f the displacement u~ is shown along the a~iflcial boundary r = 3 in the fi~t quadrant.
276
D. Oivoli, J.B. Keller / Non-reflecting boundaries
The orthogonality of the trigonometric functions yields
M - Tr(R, 0) = ~,' [Z~'I(cos nOAh+sin nOB"r)+nZ~2(sin n O A h - c o s nOB~)],
n=0 M
-
(A3)
To(R, 0) = ~' [Z~2(cos nO A~ + sin nO B~) + nZ~'l(sin nO A~ - cos nO BT)].
n=0
Next we write Z~ as a finite sum of even powers of n, namely
M
n
Zij =
Z
m
m=0
Oto (n)2m;
n =0, 1 , . . . , M.
(A4)
For a fixed i and j, (A4) is a system of M + 1 linear algebraic equations for the coefficients ot~. More explicitly, this system has the form 1o
12
14
...
o M E M 4 ...
M2
ol b
Zb
La,"MJ
/z',T]
.
(A5)
Here 0° is defined to be 1. We now substitute the expansion (A4) into (A3). In doing so we make use of the following identities: (-1) m
(cos nO) = n 2m cos nO,
( - 1 ) m+l d02m+l d 2 ~ + ' (COS nO) = n 2m+' sin nO,
2~ (-1)"
(A6)
d2m+l (sinnO)=n2msinnO,
(--1)"~TX~--d~+l(sinnO)=n2"+~ cosnO. dlt
Then we obtain - T , ( R , 0) = Y.'
n=0 m = O
o~~ ( - 1)"
(A'~cosnO+B'~sinnO)
dEm+l
]
+ t~2(--1) m+l d02m+l (A~ cos nO+B~ sin nO) ,
M M[ - T o ( R , O ) = ~' Y. a ~ ( - 1 ) m n=O
(A7)
2~ (A~cosnO+B;sinnO)
m=O
d2m+l
]
+ a ~ ( - 1 ) m+~ d02m+~ (Am cos nO+B'~ sin nO) . But from (A1) this reduces to
d2mur
md2m÷%]
- T r = m=o ~ (--1)m a ~ 2dO ---~-a12 d--~J F ~ d2"Uo - T o = ,=0 ~ (-1)"L~22~-~21
'
.,
(A8) ~
J"
Equations (A8) are the desired local boundary condition. Since (A8) follows from the exact boundary condition, it is exact for waves which consist of the first M harmonics.
D. Givoli,J.B. Keller / Non-reflectingboundaries
277
We shall now derive explicit expressions for the first two boundary conditions in (A8). With M = 0, the system (A5) reduces to 0
O
(A9)
a o - Z u. From (A2), (25) and (18) we deduce E
Z°'-(l+~~)~Y
/' k H o L + l ~ L'~LOL "~'],
z%
=
z °, =
0,
(A10) E l Z°2 = (i ~-~v) (~kT ,--~,r+ I ~ Therefore, for M = 0 the boundary condition (A8) becomes E
(
HOL I \
E
- T , . = ( l + v ) TkL-'~OL+-~)Ur,
(! k HOT+I__~
- T o = ( l + v ) \ 2 THe)T RJUO.
(All)
This boundary condition is perfectly absorbing for axially symmetric cylindrical waves. For M = I, the system (AS) yields
o tou_- Z oo,
o ~l _u - Z u1 - Z ou.
(A12)
From (A8), the local boundary condition is d2ur - Tr = Z °, Ur -- ( Z l , - Z ° , ) - ~ +
t (Z,2
-
d3ur Z°2) dO3 , (A13)
-To=Z°2uo-(Z212-
z O ~ cl u e + : Z ~
d3uo
22, d02 t 2,-Z°,) dO3-
This condition is perfectly absorbing for axially symmetric waves as well as for those involving sin 0 and cos 0. When kLR and kTR are large, we can replace the Hankel functions by their asymptotic forms. For example, by using HoL/H~L~--i and HoT/H~T~--i in (All) wc get -T,=
(l+v)
u,,
TkLi-
-To=-(l+v-----j ~kTl----~ uo.
(A14)
To express the preceding boundary conditions in terms of Cartesian quantities we use vector transformations, as in Section 2. For example (All) and (A14) have the forms
-Tr(O)=fllur(O),
-To(O)=fl2uo(O).
(A15)
After transformation they become
- T~(O) = Aij(O) uj(O),
(A16)
where 2
Aq(O) = • J~k Vik(O) Vjk(O). k=l
The functions V0(0 ) are defined in (28).
(A17)
278
D. Givoli, J.B. Keller / Non-reflecting boundaries
For the sake o f completeness, we now use a similar procedure to localize the exact boundary condition derived in Keller and Givoli [8] for the reduced wave equation. This exact boundary condition is m " ( O - 0') u ( R , 0') R dO',
u~(R, O) = - Z'
(A18)
n=0
k H~I)'(kR) ~rg H ~ ) ( k g ) cos n ( O - 0').
m " ( O - 0') =
(A19)
Here u~ is the outward normal derivative of the unknown scattered field u, and k is the wave number. We write m " ( O - 0') =1_1_ Z "
~R
cos n(O - 0'),
H~)'(kR) Z " = - k H~1)(kR) .
(A20, A21)
Proceeding as above we obtain the local boundary condition M .
-m
~2 m U
--U~= ~ (--1) ama--~,
(A22)
m=0
where am is found from the solution of the system (A5), ignoring the indices ij. For M = 0 we have the boundary condition H~ol)'(kR) - u~ = - k H
(A23)
When k R is large we can use the asymptotic approximation H ' o / H o - i to obtain (A24)
- u~ = - iku.
This boundary condition has exactly the same form as the Sommerfeld radiation condition. In the case considered at the end of Section 3 where the original domain is finite, equations (A1)-(A13) hold with the Hankel function H~ ~) replaced everywhere by the Bessel function Jn. Thus, in the boundary condition ( A l l ) all the functions H are replaced by J. To obtain an asymptotic boundary condition analogous to (A14), we use the asymptotic approximations JoL/J~L~--cot(kLR--~/4 ) and JOT/J~T~ - c o t ( k T R -~r/4), which are valid for large values of kLR and kTR. Then we get the boundary condition T~-
-To
(l+v)
(, ootb
).
(1+,,)
We note that if a local boundary condition of the form (A16) is used instead of the exact boundary condition (4), the finite element formulation is exactly the same, except that the definition of b ( . , • ) becomes b(w, u) = J~ wiAijuj d ~ .
(A26)
Although (A26) has the same form as (39), the operator Aq in (A26) is local whereas M~j in (39) is not.
Acknowledgments This work was supported in part by ONR, AFOSR and NSF. The first author was also supported by the Technion's postdoctoral grant no. 27000827 and by the L. Kraus Research Fund.
D. Givoli, J.B. Keller / Non-reflecting boundaries
279
References [1] B. Engquist and A. Majda, "Absorbing boundary conditions for the numerical simulation of waves," Math. Comp. 31, 629-652 (1977). [2] A. Bayliss and E. Turkel, "Radiation boundary conditions for wave-like equations," Comm. Pure Appl. Math. 33, 707-725 (1980). [3] K. Feng, "Finite element method and natural boundary reduction," Proceedings of the International Congress of Mathematicians, Warsaw, 1439-1453 (1983). [4] L. Ting and M.J. Miksis, "Exact boundary conditions for scattering problems," J. Acoust. Soc. Am. 80, 1825-1827 (1986). [5] R.L. Higdon, "Numerical absorbing boundary conditions for the wave equation," Math. Comp. 49, 65-90 (1987). [6] R.C. MacCamy and S.P. Matin, "A finite element method for exterior interface problems," Int. J. Math. Math. Sci. 3, 311-350 (1980). [7] C.I. Goldstein, "A finite element method for solving Helmholtz type equations in wave guides and other unbounded domains," Math. Comp. 39, 309-324 (1982). [8] J.B. Keller and D. Givoli, "Exact non-reflecting boundary conditions," J. Comp. Phys. 82, 172-192 (1989). [9] J. Lysmer and R.L. Kuhlemeyer, "Finite dynamic model for infinite media," J. Eng. Mech. Div. ASCE 95, 859-877 (1969). [10] W.D. Smith, "A nonreflecting plane boundary for wave propagation problems," J. Comp. Phys. 15, 492-503 (1974). [11] C. Cerjan, D. Kosloff, R. Kosloff and M. Reshef, "A nonreflecting boundary condition for discrete acoustic and elastic wave equations," Geophysics 50, 705-708 (1985). [12] J. Sochacki, R. Kubicheck, J. George, W.R. Fletcher and S. Smithson, "Absorbing boundary conditions and surface waves," Geophysics 52, 60-71 (1987). [13] B. Engquist and A. Majda, "Radiation boundary conditions for acoustic and elastic wave calculations," Comm. Pure App. Math. 32, 313-357 (1979). [14] R.W. Clayton and B. Engquist, "Absorbing boundary conditions for acoustic and elastic wave equations," Bull. Seismol. Soc. Amer. 67, 1529-1540 (1977). [15] C.L. Scandrett, G. A. Kriegsmann and J.D. Achenbach, "Application of the limiting amplitude principle to elastodynamic scattering problems," S I A M J. Sci. Star. Comp. 7, 571-590 (1986). [16] R.L. Higdon, "Absorbing boundary conditions for elastic waves," SIAM J. Num. Anal., to appear. [17] D. Givoli and J.B. Keller, "A finite element method for large domains," Comp. Math. Appl. Mech. Eng. 76, 41-66 (1989).