- Email: [email protected]
A variationally coupled finite element-boundary element method
004s7949/89 53.00 + 0.00 Q 1989Pergamon PM plc
Compunrs d Strucrures Vol. 33.No. I. pp. 17-20,1989 Printedin Great Britain.
A VA~IATIONALLY ELEMENT-BOUNDARY T. BELYWHKO,
COUPLED FINITE ELEMENT METHOD
H. S.
CHANG
and Y. Y. Lu
Departments of Civil and Mechanical Engineering, The Technological Institute, Northwestern University, Evanston, IL 60208, U.S.A. (Received 27 May 1988) Abstract-Avariationally based coupling is developed for the finite element and boundary element
methods. This is achieved by combining the variational forms for the boundary element and finite element subdomains to obtain a global variational form and then choosing a suitable set of test and trial functions. This method is illustrated in the context of the Laplace equation, and a numerical example is given.
#,j nj -=u,,=4
INTRODUCTION The coupling of finite element (FE) and boundary element (BE) methods is usually accomplished at the level of the discretized equations. Such coupling methods were first developed by Zienkiewicz et a/. [ 11;
a good review may be found in Brebbia et al. [2]. Even the recently proposed “global coupling” approach described by Li et al. [3] shares this characteristic of coupling at the level of the discretized equations. The purpose of this paper is to present a methodology whereby the discrete equations for a coupled FE-BE discretization are obtained directly from a single variational statement. This is accomplished by combining the variational forms which apply to the FE and BE discretization procedures and then using a suitable set of test and trial functions which span the entire domain. This leads to a more aesthetically appealing procedure and eliminates the need for treating interface nodes twice at the algebraic level. In dynamics, it may lead to improved numerical behavior as well, although this aspect is not discussed in this paper. In the following, the variationally coupled FE-BE procedure is described for the Laplace equation. The simplicity of this setting permits us to cIearly convey the essential features of the variationally coupled procedure. In order to demonstrate its soundness, some numerical results are also presented.
DEVELOPMENT OF
VARIATIONAL
onr,,
where the superposed bar denotes prescribed boundary values. Standard index notation is used, so repeated indices designate summations and a comma designates a partial derivative with respect to the indicated spatial variable; n, is the unit normal to the domain R. The domain a is subdivided into Q, and rZ,, as shown in Fig. 1, which are to be solved by finite elements and boundary elements, respectively. The interface between Cl, and f2, is r, ; the boundary of the remainder of QB is denoted by fe, and the boundary of the remainder of OF is denoted by F,. rF is subdivided into the prescribed function boundary rUF and the prescribed gradient boundary rrlF. The variational (or weak) forms of the governing equations in the two su~omains are given by the following:
s rLld)+rf
(%i - vu,,,) dT =
ui(,” s
r*F+rr
dr =
(@@* ii- u”,ti) dQ (2)
f Qe
u,iu*tdn*
I nF
(ue,, - vu,&) dI’ (L(t),,- OU,,)dT + I r: I r.9
We consider the following problem on the domain $2 bounded by r, in R
(la)
U-U
on r,
(lb)
(3)
where u are the test functions and u is the approximate solution. In order to more clearly bring out the terms on the interface, we write out (2) and (3) as
FORM
u,, = 0
(lc)
s rt 17
T.
18
BELYTSCHKOet
al.
where N;(x) are shape functions which are nonzero only in element e, and L;, is the Boolean connectivity matrix of element e; repeated upper-case subscripts are summed over the total number of nodes in the model m. The test functions are now defined as follows: v:(x)
ifx,Er,Ur,
v,(x) = { N,(x) Fig. 1.Nomenclature for domain subdivided into boundary element (B) and finite element (F) subdomains.
Now, noting that nf = --nP on r,, we eliminate this term from (4) by using (5), which gives
if x1 E S2, U rqF ’
where v?(x) = WF(x)R(x)
f I rs
(uo,ii-Uu,ii)dQ-
(13)
o,iu,idQ
(6)
I nF
s n,
where the third term has been modified by invoking the usual weak form of the natural boundary condition
s
v(g-u,.)dI-
=O.
(12)
1 in R,
(UC, -vu,.)dr.+~~,uu,.dT-~~~~ujdT
=
(11)
Note that the test functions are continuous across r,. Moreover, they have been constructed so that v: = WY on r, and so that VT vanishes at all nodes inside the finite element domain. The approximate solution (trial function) is constructed as follows:
(7)
u(x) = u,N,(x)
inCl,Ur,Ur,,
(14a)
For problems where the boundary element domain is infinite and encloses the finite element domain so that Te is absent, (6) can be simplified to
u(x) = &N,(x)
on ruf
(14b)
in RB
(14c)
on rq8
(W
on Tue.
(14e)
r,f
u,ii= 0 u(x) = u,N,(x) q(x) = &N,(x) >
Equation (6) represents a single variational form which holds for both the boundary element and finite element domains. For the purpose of implementing this variational form, it only remains to define appropriate test functions v and trial functions u for both domains. It is therefore useful to define the fundamental solution (Green’s function) for the infinite domain WF which satisfies
WY9 ii= 6tx- l,h
u(x) = C,N,(x) q(x) = q/N,(x) > Both the boundary and finite element solutions are represented by (14a) on the interface r,, so they are obviously continuous at the interface. The discrete equations for the coupled system are obtained by substituting eqns (11) and (14) into (6), which yields
(9)
where 6 is the Dirac delta function, and 4, is an arbitrary point in SL A set of M nodes is now placed on fs + r,+ r,+ R,. The nodes are numbered so that the first ms nodes are on Ts. the new m, nodes are on r,, and the remaining mF nodes are on RF and fF; the total number of nodes is m = mB + m, + mF. The global shape functions are denoted by N,(x), and they are related to the element shape functions by (10) N,(x) = 1 LLNXx),
Hu+Cq+Ku=f,
(15)
where for 1 < mg + m,,
(16)
s s
WY NJ dl-
G/J= -
KIJ
=
r.s
VlfiNJ,idR
G
(17)
(18)
19
A variationally coupled finite element-boundary element method
. 1.0
-
W&ii,N,dl-
i ru.9
+
u:qdT
(19)
I r#
h a ; 5
and, form,+m,~I~m,
-.t +
24n9l~csl Numerial ml.0 Numerical rd.9 Nummkd rd.8
+
Numuiul ~0.7
0.8
0.6
H,,=G,,=O
(20) 0.4
Ku =
h=
s nF
I r,f
N,,iN,,im
(21)
N,4 dr,
(22)
0.0 0.0
0.2
0.1
0.3
I
0.4
,3..5
w*@
where g = f on r,, U r, and g = 0 on TUB.The arrays u and q are both of length m, but as is conventional
in BE techniques, only m of these values are unknown. In eqn (18), the domain has been restricted to the elements in the domain contiguous to T,, which is denoted by a,,, since OF vanishes in the remainder of the domain, as can be seen from (12) and (13). The matrices H and G are the usual boundary element matrices; only (18) is different. The equations for the interior nodes of the finite element domain retain their usual symmetric, banded structure. The equations associated with the boundary element nodes are full and unsymmetric. The element counterparts of eqns (18) and (22) are obtained by replacing N with N’; thus (18) becomes K; =
0.2
vrfiN;,i dR.
(23)
s c)r
Note that this element matrix has contributions nodes on r,.
with analytical
Fig. 3. Results for example compared solution.
NUMERICAL EXAMPLES
A simple example has been included to verify the variationally coupled method. For this purpose, eqn (la) on a circular domain has been chosen. The inner domain is treated by BE with 24 nodes on the interface; the outer domain is treated by a mesh of 24 x 3 quadrilateral finite elements. The mesh is shown in Fig. 2. On the outer boundary, r = 1.0, u is prescribed, so it is a r, boundary. Two boundary conditions are considered:
(a)
u(r = 1, e) = 1.0
(W
u(r = 1,e) = sin 28.
to all
Y
The first boundary condition was considered simply to test whether the method satisfies the Laplacian counterpart of the rigid body motion test: do the products of G and K with any vector corresponding to a constant field vanish? Our results show that they do and that the solution for interior nodes is 1.0 to five significant digits on a 14-digit computer. This
Table I. Results for example with P = sin 26 on rU r
e
0.7
0 n/12 xl6 xl4 0 n/12 xl6 xl4 0 x/12 x/6 xl4
0.8
V2”.0
fg - 0.7
(1) u(1.B) = 1.0
*###. 1.0
(2) u(M)
MeshforFEM:24x3 Elammofw BBM: 24
= silts)
Fig. 2. Mesh and problem statement for example.
0.9
u
FEYBEexact
IA
error
0.83 x 1O-6 0.2409 0.4172 0.4817 0.2008 X 10-5 0.3172 0.5493 0.6343
i.245 0.424 0.49 0 0.32 0.554 0.64
0.85% 0.90% 0.89%
0.17 x 10-5 0.4035 0.6989 0.8070
0 0.405 0.701 0.81
0.37% 0.37% 0.37%
1.67% 1.69% 1.69%
T. BELYIXCHKO et al.
20
inci~~~lly provides a simple debugging test for the coupled method: the matrices K and H and their element constituents, G’ and H’, must be such that their product with any vector obtained from a uniform field must vanish. Results for the second problem are shown in Fig. 3. For a closer look, some numerical values and errors are given in Table 1. The results show excellent agreement with the closed-form
solution.
Acknowledgement-The support of the Army Research Otlice under Grant DAAL03-87-K-0035 to Northwestern University is gratefully acknowledged.
1. 0. C. Zienkiewicz, D. W. Kelly and P. Bettess, The coupling of the finite element method and boundary solution procedures. Inf. J. Numer. Merh. Engng 11, 355-375 (1977).
2. C. A. Brebbia, J. C. F. Telles and L. C. Wrobel, Boundary Element Techniques, Chapter 13. SpringerVerlag, Berlin (1984). 3. H.-B. Li, G.-M. Han, H. A. Mang and P. Torzicky, A new method for the coupting of finite element and boundary element discretized subdomains of elastic bodies. Comput. Meth. appl. Mech. Engng 54, 161-185 (1986).
Compunrs d Strucrures Vol. 33.No. I. pp. 17-20,1989 Printedin Great Britain.
A VA~IATIONALLY ELEMENT-BOUNDARY T. BELYWHKO,
COUPLED FINITE ELEMENT METHOD
H. S.
CHANG
and Y. Y. Lu
Departments of Civil and Mechanical Engineering, The Technological Institute, Northwestern University, Evanston, IL 60208, U.S.A. (Received 27 May 1988) Abstract-Avariationally based coupling is developed for the finite element and boundary element
methods. This is achieved by combining the variational forms for the boundary element and finite element subdomains to obtain a global variational form and then choosing a suitable set of test and trial functions. This method is illustrated in the context of the Laplace equation, and a numerical example is given.
#,j nj -=u,,=4
INTRODUCTION The coupling of finite element (FE) and boundary element (BE) methods is usually accomplished at the level of the discretized equations. Such coupling methods were first developed by Zienkiewicz et a/. [ 11;
a good review may be found in Brebbia et al. [2]. Even the recently proposed “global coupling” approach described by Li et al. [3] shares this characteristic of coupling at the level of the discretized equations. The purpose of this paper is to present a methodology whereby the discrete equations for a coupled FE-BE discretization are obtained directly from a single variational statement. This is accomplished by combining the variational forms which apply to the FE and BE discretization procedures and then using a suitable set of test and trial functions which span the entire domain. This leads to a more aesthetically appealing procedure and eliminates the need for treating interface nodes twice at the algebraic level. In dynamics, it may lead to improved numerical behavior as well, although this aspect is not discussed in this paper. In the following, the variationally coupled FE-BE procedure is described for the Laplace equation. The simplicity of this setting permits us to cIearly convey the essential features of the variationally coupled procedure. In order to demonstrate its soundness, some numerical results are also presented.
DEVELOPMENT OF
VARIATIONAL
onr,,
where the superposed bar denotes prescribed boundary values. Standard index notation is used, so repeated indices designate summations and a comma designates a partial derivative with respect to the indicated spatial variable; n, is the unit normal to the domain R. The domain a is subdivided into Q, and rZ,, as shown in Fig. 1, which are to be solved by finite elements and boundary elements, respectively. The interface between Cl, and f2, is r, ; the boundary of the remainder of QB is denoted by fe, and the boundary of the remainder of OF is denoted by F,. rF is subdivided into the prescribed function boundary rUF and the prescribed gradient boundary rrlF. The variational (or weak) forms of the governing equations in the two su~omains are given by the following:
s rLld)+rf
(%i - vu,,,) dT =
ui(,” s
r*F+rr
dr =
(@@* ii- u”,ti) dQ (2)
f Qe
u,iu*tdn*
I nF
(ue,, - vu,&) dI’ (L(t),,- OU,,)dT + I r: I r.9
We consider the following problem on the domain $2 bounded by r, in R
(la)
U-U
on r,
(lb)
(3)
where u are the test functions and u is the approximate solution. In order to more clearly bring out the terms on the interface, we write out (2) and (3) as
FORM
u,, = 0
(lc)
s rt 17
T.
18
BELYTSCHKOet
al.
where N;(x) are shape functions which are nonzero only in element e, and L;, is the Boolean connectivity matrix of element e; repeated upper-case subscripts are summed over the total number of nodes in the model m. The test functions are now defined as follows: v:(x)
ifx,Er,Ur,
v,(x) = { N,(x) Fig. 1.Nomenclature for domain subdivided into boundary element (B) and finite element (F) subdomains.
Now, noting that nf = --nP on r,, we eliminate this term from (4) by using (5), which gives
if x1 E S2, U rqF ’
where v?(x) = WF(x)R(x)
f I rs
(uo,ii-Uu,ii)dQ-
(13)
o,iu,idQ
(6)
I nF
s n,
where the third term has been modified by invoking the usual weak form of the natural boundary condition
s
v(g-u,.)dI-
=O.
(12)
1 in R,
(UC, -vu,.)dr.+~~,uu,.dT-~~~~ujdT
=
(11)
Note that the test functions are continuous across r,. Moreover, they have been constructed so that v: = WY on r, and so that VT vanishes at all nodes inside the finite element domain. The approximate solution (trial function) is constructed as follows:
(7)
u(x) = u,N,(x)
inCl,Ur,Ur,,
(14a)
For problems where the boundary element domain is infinite and encloses the finite element domain so that Te is absent, (6) can be simplified to
u(x) = &N,(x)
on ruf
(14b)
in RB
(14c)
on rq8
(W
on Tue.
(14e)
r,f
u,ii= 0 u(x) = u,N,(x) q(x) = &N,(x) >
Equation (6) represents a single variational form which holds for both the boundary element and finite element domains. For the purpose of implementing this variational form, it only remains to define appropriate test functions v and trial functions u for both domains. It is therefore useful to define the fundamental solution (Green’s function) for the infinite domain WF which satisfies
WY9 ii= 6tx- l,h
u(x) = C,N,(x) q(x) = q/N,(x) > Both the boundary and finite element solutions are represented by (14a) on the interface r,, so they are obviously continuous at the interface. The discrete equations for the coupled system are obtained by substituting eqns (11) and (14) into (6), which yields
(9)
where 6 is the Dirac delta function, and 4, is an arbitrary point in SL A set of M nodes is now placed on fs + r,+ r,+ R,. The nodes are numbered so that the first ms nodes are on Ts. the new m, nodes are on r,, and the remaining mF nodes are on RF and fF; the total number of nodes is m = mB + m, + mF. The global shape functions are denoted by N,(x), and they are related to the element shape functions by (10) N,(x) = 1 LLNXx),
Hu+Cq+Ku=f,
(15)
where for 1 < mg + m,,
(16)
s s
WY NJ dl-
G/J= -
KIJ
=
r.s
VlfiNJ,idR
G
(17)
(18)
19
A variationally coupled finite element-boundary element method
. 1.0
-
W&ii,N,dl-
i ru.9
+
u:qdT
(19)
I r#
h a ; 5
and, form,+m,~I~m,
-.t +
24n9l~csl Numerial ml.0 Numerical rd.9 Nummkd rd.8
+
Numuiul ~0.7
0.8
0.6
H,,=G,,=O
(20) 0.4
Ku =
h=
s nF
I r,f
N,,iN,,im
(21)
N,4 dr,
(22)
0.0 0.0
0.2
0.1
0.3
I
0.4
,3..5
w*@
where g = f on r,, U r, and g = 0 on TUB.The arrays u and q are both of length m, but as is conventional
in BE techniques, only m of these values are unknown. In eqn (18), the domain has been restricted to the elements in the domain contiguous to T,, which is denoted by a,,, since OF vanishes in the remainder of the domain, as can be seen from (12) and (13). The matrices H and G are the usual boundary element matrices; only (18) is different. The equations for the interior nodes of the finite element domain retain their usual symmetric, banded structure. The equations associated with the boundary element nodes are full and unsymmetric. The element counterparts of eqns (18) and (22) are obtained by replacing N with N’; thus (18) becomes K; =
0.2
vrfiN;,i dR.
(23)
s c)r
Note that this element matrix has contributions nodes on r,.
with analytical
Fig. 3. Results for example compared solution.
NUMERICAL EXAMPLES
A simple example has been included to verify the variationally coupled method. For this purpose, eqn (la) on a circular domain has been chosen. The inner domain is treated by BE with 24 nodes on the interface; the outer domain is treated by a mesh of 24 x 3 quadrilateral finite elements. The mesh is shown in Fig. 2. On the outer boundary, r = 1.0, u is prescribed, so it is a r, boundary. Two boundary conditions are considered:
(a)
u(r = 1, e) = 1.0
(W
u(r = 1,e) = sin 28.
to all
Y
The first boundary condition was considered simply to test whether the method satisfies the Laplacian counterpart of the rigid body motion test: do the products of G and K with any vector corresponding to a constant field vanish? Our results show that they do and that the solution for interior nodes is 1.0 to five significant digits on a 14-digit computer. This
Table I. Results for example with P = sin 26 on rU r
e
0.7
0 n/12 xl6 xl4 0 n/12 xl6 xl4 0 x/12 x/6 xl4
0.8
V2”.0
fg - 0.7
(1) u(1.B) = 1.0
*###. 1.0
(2) u(M)
MeshforFEM:24x3 Elammofw BBM: 24
= silts)
Fig. 2. Mesh and problem statement for example.
0.9
u
FEYBEexact
IA
error
0.83 x 1O-6 0.2409 0.4172 0.4817 0.2008 X 10-5 0.3172 0.5493 0.6343
i.245 0.424 0.49 0 0.32 0.554 0.64
0.85% 0.90% 0.89%
0.17 x 10-5 0.4035 0.6989 0.8070
0 0.405 0.701 0.81
0.37% 0.37% 0.37%
1.67% 1.69% 1.69%
T. BELYIXCHKO et al.
20
inci~~~lly provides a simple debugging test for the coupled method: the matrices K and H and their element constituents, G’ and H’, must be such that their product with any vector obtained from a uniform field must vanish. Results for the second problem are shown in Fig. 3. For a closer look, some numerical values and errors are given in Table 1. The results show excellent agreement with the closed-form
solution.
Acknowledgement-The support of the Army Research Otlice under Grant DAAL03-87-K-0035 to Northwestern University is gratefully acknowledged.
1. 0. C. Zienkiewicz, D. W. Kelly and P. Bettess, The coupling of the finite element method and boundary solution procedures. Inf. J. Numer. Merh. Engng 11, 355-375 (1977).
2. C. A. Brebbia, J. C. F. Telles and L. C. Wrobel, Boundary Element Techniques, Chapter 13. SpringerVerlag, Berlin (1984). 3. H.-B. Li, G.-M. Han, H. A. Mang and P. Torzicky, A new method for the coupting of finite element and boundary element discretized subdomains of elastic bodies. Comput. Meth. appl. Mech. Engng 54, 161-185 (1986).