Interface relaxation algorithms for BEM–BEM coupling and FEM–BEM coupling

Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992 www.elsevier.com/locate/cma

Interface relaxation algorithms for BEM–BEM coupling and FEM–BEM coupling Wael M. Elleithy, Masataka Tanaka

*

Department of Mechanical Systems Engineering, Faculty of Engineering, Shinshu University, 4-17-1 Wakasato, Nagano 380-8553, Japan Received 16 October 2002; received in revised form 6 February 2003

Abstract This paper presents several interface relaxation algorithms for boundary element–boundary element coupling (BEM–BEM) and for finite element–boundary element coupling (FEM–BEM). The domain of the original problem is sub-divided into sub-domains, which are modeled by the finite element or boundary element methods. The multidomain system is coupled using smoothing operators on the inter-domain boundaries. Separate computations for the BEM and FEM sub-domains and successive update of the boundary conditions at the interfaces are performed until convergence is achieved. The interface relaxation coupling algorithms preserve the nature of the FEM and BEM. Further, they do not require any access to the matrices generated by the FEM or BEM and make it easier to utilize different software in different sub-domains. Ó 2003 Elsevier B.V. All rights reserved. Keywords: Boundary element method; Finite element method; Interface relaxation; Coupling

1. Introduction The finite element method (FEM) and the boundary element method (BEM) are well known as powerful numerical techniques for solving wide range of problems in applied science and engineering. Both the FEM and the BEM have their own range of applications where they are most efficient and there are, undoubtedly, situations that favor FEM over BEM and vice versa. In general the FEM is best suited for solving problems with strong non-linearity in bodies of complex shapes but finite dimensions. The BEM on the other hand is very efficient for the analysis of infinite linear problems. It is attractive to decompose a computational domain into sub-domains and to use the most appropriate solution technique for each sub-domain (utilize different ‘‘legacy’’ software for solving the individual subproblems).

*

Corresponding author. Tel.: +81-26-269-5120; fax: +81-26-269-5124. E-mail addresses: [email protected] (W.M. Elleithy), [email protected] (M. Tanaka).

0045-7825/03/$ - see front matter Ó 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0045-7825(03)00312-8

2978

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

The idea of coupling the FEM and BEM is by now well known as an effective analysis tool, which makes use of their advantages. The first proposed FEM and BEM coupled formulation was presented by Zienkiewicz et al. [1]. An extensive literature survey on this topic can be found in [2–5]. The systems of equations produced by the two methods are expressed in different variables and cannot be linked as they stand. The conventional coupling methods alter the formulation of one of the methods to make it compatible with the other. The conventional coupling methods, however, may destroy the desirable features originally existing in the FEM matrices, namely, symmetry, sparsity and bandedness. Moreover, the implementation of the conventional coupling procedures requires access to the source codes and necessitates merging two different kinds of programs. In order to preserve the nature of both the FEM and BEM, the domain decomposition coupling methods have been developed [6–10]. In these methods, separate computation for the BEM and FEM sub-domains and successive update of the interface boundary conditions are performed until convergence is achieved. The domain decomposition coupling methods offer many advantages over the conventional coupling methods and appear to be promising. However, there are several challenging questions concerning the practical applications of such methods, e.g., finding the most suitable method for a particular application and establishing the convergence and optimal convergence conditions of the domain decomposition coupling methods. In boundary element analysis, sub-domain partition may be employed when the domains under consideration are governed by individual differential equations and/or constructed of different materials. Besides, in the case of domain with complicated boundary profile, the domain may be decomposed for better computational efficiency and accuracy. In the conventional BEM–BEM coupling formulation, subdomains equations are collected into a system of equations while satisfying the continuity and equilibrium conditions along the interfaces. However, the conventional methods require access to the source codes. Ref. [14] describes the Neumann–Neumann and Dirichlet–Neumann domain decomposition BEM–BEM coupling algorithms. In this paper, we present a number of interface relaxation coupling algorithms, which will allow the user to find the most suitable relaxer for a particular application. The algorithms have numerous attractive features, e.g., no preconditioning is needed, and they do not overlap sub-domains, are general enough and make it easy to utilize different software in different sub-domains. The paper is arranged as follows: Section 2 gives a brief presentation of the interface relaxation methods for the solution of composite PDEs. Sections 3 and 4 present the new interface relaxation algorithms for FEM–BEM coupling and BEM–BEM coupling, respectively. Numerical examples and guidelines for selecting the relaxation parameters are given in Sections 5 and 6, respectively.

2. Interface relaxation The domain decomposition methods can be classified as overlapping and non-overlapping. The overlapping, known also as Schwarz, methods were the first considered and have already proved themselves as efficient numerical procedures enjoying certain desirable convergence properties. However, the overlapping methods are not applicable to general composite PDE problems. Further, they may create serious complications in the Schwarz method even when the problem is of a simple geometry. The non-overlapping methods exhibit advantages compared to the overlapping ones and there are two principle viewpoints of the non-overlapping methods: preconditioning and interface relaxation. Interface relaxation is more general than the traditional domain decomposition methods, see Refs. [15– 19]. It follows SouthwellÕs relaxation, but at the PDE instead of the linear algebra level, to formulate relaxation as iterated interface smoothing procedures. Mu [15] presented a general framework for solving composite PDEs based on interface relaxation and it may be summarized as

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

2979

1. Denote the local PDEs by Li ui ¼ fi in Xi , for i ¼ 1; 2; . . . ; m, where Li is a differential operator and m is the number of sub-domains. 2. Define Cij as the interface between two adjacent sub-domains Xi and Xj . Let SðiÞ be the indices of those sub-domains that are neighbors of sub-domain Xi . Define the boundary value problem Pi;n that is solved on Xi at the nth relaxation step as Li ui;n ¼ fi

Bij ui;n ¼ bij;n Bi ui;n ¼ gi

ð1Þ

in Xi ; on Cij for j 2 SðiÞ;

ð2Þ

on Ci \ C;

ð3Þ

where Bi is a boundary condition operator, gi is a continuous function and Bij is a boundary condition operator such that Pi;n is well posed. The interface relaxation iteration is defined on the sub-domains independently. Details of the iteration are specified by an interface handler called a relaxer. It provides for the interface Cij to sub-domain Xi the right hand side data bij;n of the boundary condition according to certain relaxation procedures as well as the parameters in the definition of Bij . Mathematically this framework contains many existing domain decomposition methods and also allows the extension to a variety of new relaxers. Rice et al. [19] presented an interface relaxation algorithm for the solution of elliptic differential equations. The algorithm defines Bij as the Dirichlet operator and   oui;n1 ouj;n1 bij;n ¼ ui;n1  a þ on Cij for j 2 SðiÞ; on on where a is a relaxation parameter to ensure and/or accelerate convergence. Ref. [17] presents an interface relaxation algorithm that defines Bij as the Robin operator and bij;n ¼ 

ouj;n1 þ quj;n1 on

on Cij for j 2 SðiÞ;

where q is a relaxation parameter. The Dirichlet/Neumann average algorithm [19] consists of two PDE solving sweeps coupled with smoothing interface relaxation steps. It defines bij;n ¼ ð1  /2 Þui;n1 þ u2 uj;n1

on Cij for j 2 SðiÞ ðBij is the Dirichlet operatorÞ

and bij;nþ1 ¼ ð1  u1 Þ

oui;n ouj;n  u1 on on

on Cij for j 2 SðiÞ ðBij is the Neumann operatorÞ;

u1 and u2 are relaxation parameters.

3. FEM–BEM coupling Consider Fig. 1, where the domain of the original problem is decomposed into FEM and BEM subdomains. The corresponding boundary integral equation for the BEM sub-domain is given by ½H fug ¼ ½G fqg

on CB ;

ð4Þ

where u and q are column matrices containing the boundary nodal values for the potentials (displacements) and the fluxes (tractions). H and G are influence coefficient matrices. For the FEM sub-domain, the assembled element equations are given by ½K fug ¼ ff g

in XF ;

ð5Þ

2980

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

Fig. 1. (a) Domain decomposed into FEM and BEM sub-domains. (b) FEM modeling. (c) BEM modeling.

where K is the stiffness matrix for the system, and u and f are the nodal potentials (displacements) and integrated flux (force) vectors respectively. Now, let us define the following vectors: uIB : uBB : uIF : uFF :

interface potentials (displacements), approached from the BEM sub-domain, non-interface potentials (displacements) in the BEM sub-domain, interface potentials (displacements), approached from the FEM sub-domain, non-interface potentials (displacements) in the FEM sub-domain.

Similarly, one can define the flux (tractions) and the integrated flux (force) vectors for the BEM and FEM sub-domains, respectively. Eqs. (4) and (5) may be partitioned as follows:   B    B  uB qB H11 H12 G11 G12 ¼ ð6Þ H21 H22 G21 G22 uIB qIB and 

K11 K21

K12 K22



uFF uIF



 ¼

 fFF : fFI

ð7Þ

At the interface, the compatibility and equilibrium conditions should be satisfied, i.e., fuIB g ¼ fuIF g

on CI ;

ð8Þ

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

ffFI g þ ½M fqIB g ¼ 0 on CI ;

2981

ð9Þ

where M is the converting matrix, which depends on the interpolation functions used to represent the fluxes (tractions) on the interface. In Section 3.1 we briefly review the existing domain decomposition FEM–BEM coupling algorithms. In Sections 3.2–3.4, we present several new interface relaxation algorithms for FEM–BEM coupling. 3.1. Domain decomposition FEM–BEM coupling algorithms Kamiya et al. [8] employed the Neumann–Neumann and Dirichlet–Neumann algorithms for the coupling of FEM and BEM. Lin et al. [9] and Feng and Owen [10] presented a sequential form of the Dirichlet– Neumann FEM–BEM coupling algorithm. 3.1.1. Sequential Dirichlet–Neumann FEM–BEM coupling algorithm The sequential Dirichlet–Neumann FEM–BEM coupling algorithm may be described as follows [9,10]: Set the initial guess fuIB;0 g. Do for n ¼ 0; 1; 2; . . . for the BEM sub-domain:  B    B  uB;n qB;n H11 H12 G11 G12 solve ¼ for fqIB;n g H21 H22 G21 G22 uIB;n qIB;n for the FEM sub-domain: I I g þ ½M fqIB;n g ¼ 0 for ffF;n g solve ffF;n   F   F  uF;n fF;n K11 K12 solve ¼ for fuIF;n g I K21 K22 uIF;n fF;n apply fuIB;nþ1 g ¼ ð1  hÞfuIB;n g þ hfuIF;n g where h is a relaxation parameter to ensure and/or accelerate convergence. Elleithy et al. [11,12] established the convergence and optimal convergence conditions for the sequential Dirichlet–Neumann FEM–BEM coupling algorithm. The algorithm, however, is not suitable for problems when the Neumann boundary conditions are specified on the entire real boundary of the FEM sub-domain. The application of the algorithm leads to the Neumann condition being specified on the whole boundary of the sub-domain and a system of equations cannot be solved without some additional conditions. 3.1.2. Parallel Neumann–Neumann FEM–BEM coupling algorithm The parallel Neumann–Neumann FEM–BEM coupling algorithm may be described as follows [8]: Set the initial guess fqIB;0 g and fqIF;0 g. Do for n ¼ 0; 1; 2; . . . until convergence for the BEM sub-domain:  B    B  uB;n qB;n H11 H12 G11 G12 solve ¼ for fuIB;n g I H21 H22 G21 G22 uB;n qIB;n for the FEM sub-domain: I I solve ffF;n g ¼ ½M fqIF;n g for ffF;n g   F   F  uF;n fF;n K11 K12 solve ¼ for fuIF;n g I K21 K22 uIF;n fF;n

2982

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

apply fqIB;nþ1 g ¼ fqIB;n g þ bðfuIF;n g  fuIB;n gÞ apply fqIF;nþ1 g ¼ fqIB;nþ1 g where b is a relaxation parameter to ensure and/or accelerate convergence. This algorithm is not suitable for problems when the Neumann boundary conditions are specified on the entire real boundary of the FEM or BEM sub-domains. 3.1.3. Parallel Dirichlet–Neumann FEM–BEM coupling algorithm The parallel Dirichlet–Neumann FEM–BEM coupling algorithm may be described as follows [8]: Set initial guess fuIB;0 g and fqIF;0 g. Do for n ¼ 0; 1; 2; . . . until convergence for the BEM sub-domain:   B    B  uB;n qB;n H11 H12 G11 G12 solve ¼ for fqIB;n g H21 H22 G21 G22 uIB;n qIB;n for the FEM sub-domain: I I solve ffF;n g ¼ ½M fqIF;n g for ffF;n g   F   F  uF;n fF;n K11 K12 solve ¼ for fuIF;n g I I K21 K22 uF;n fF;n apply fuIB;nþ1 g ¼ ð1  cÞfuIB;n g þ cfuIF;n g apply fqIF;nþ1 g ¼ fqIB;n g where c is a relaxation parameter to ensure and/or accelerate convergence. This coupling algorithm has the same drawback as that of the sequential Dirichlet–Neumann coupling method. 3.2. Geometric contraction based FEM–BEM coupling algorithm The geometric contraction based FEM–BEM coupling algorithm estimates new values of the potentials/ displacements by adding to the old ones, a geometrically weighted combination of the normal boundary derivatives (fluxes/tractions) of the adjacent sub-domains. Noting that q ¼ kðou=onÞ (k is a material property), the algorithm adds a correction to the vectors fuIF g and fuIB g so as to match the normal derivatives fqIF g and fqIB g at the interface. The geometric contraction based FEM–BEM coupling algorithm may be described as follows: Set initial guess fuIF;0 g and fuIB;0 g. For n ¼ 0; 1; 2; . . ., do until convergence for the BEM sub-domain:  B    B  uB;n qB;n H11 H12 G11 G12 solve ¼ for fqIB;n g H21 H22 G21 G22 uIB;n qIB;n for the FEM sub-domain:   F   F  uF;n fF;n K11 K12 I g ¼ ½M fqIF;n g for fqIF;n g solve ¼ and ffF;n I I K21 K22 uF;n fF;n apply fuIB;nþ1 g ¼ fuIB;n g  aðfqIB;n g þ fqIF;n gÞ and fuIF;nþ1 g ¼ fuIB;nþ1 g where a is a relaxation parameter to ensure and/or accelerate convergence.

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

2983

In practical applications, the value of fug and fqg may differ greatly in order of magnitude. An important aspect of this new coupling algorithm includes scaling of fug and fqg while selecting the dimensionless value of the relaxation parameter, a. 3.3. Robin relaxation FEM–BEM coupling algorithm The Robin relaxation FEM–BEM coupling algorithm predicts new values at the interface by making a convex combination of Dirichlet and Neumann data from the adjacent sub-domains. The coupling algorithm may be described as follows: I I Define fgB;nþ1 g ¼ fqIF;n g þ qfuIF;n g and fgF;nþ1 g ¼ fqIB;n g þ qfuIB;n g, where q is a relaxation parameter to ensure and/or accelerate convergence. Set initial guess fuIF;0 g and fuIB;0 g. For n ¼ 0; 1; 2; . . ., do until convergence for the BEM sub-domain and with the Robin boundary conditions on the interface:   B    B  uB;n qB;n H11 H12 G11 G12 solve ¼ for fuIB;n g and fqIB;n g I H21 H22 G21 G22 uB;n qIB;n for the FEM sub-domain  F  and  with the Robin boundary conditions on the interface: F u f K11 K12 F;n F;n I solve g ¼ ½M fqIF;n g for fuIF;n g and fqIF;n g ¼ and ffF;n I K21 K22 uIF;n fF;n I g ¼ fqIF;n g þ qfuIF;n g apply fgB;nþ1 I apply fgF;nþ1 g ¼ fqIB;n g þ qfuIB;n g:

3.4. Dirichlet/Neumann averaging FEM–BEM coupling algorithm The Dirichlet/Neumann averaging FEM–BEM coupling algorithm may be described as follows: Set initial guess fuIF;0 g and fuIB;0 g. For n ¼ 0; 2; 4; . . ., do until convergence for the BEM sub-domain:  B    B  uB;n qB;n H11 H12 G11 G12 solve ¼ for fqIB;n g I H21 H22 G21 G22 uB;n qIB;n for the FEM sub-domain:   F   F  uF;n fF;n K11 K12 I g ¼ ½M fqIF;n g for fqIF;n g solve ¼ and ffF;n I K21 K22 uIF;n fF;n apply fqIB;nþ1 g ¼ ð1  u1 ÞfqIB;n g  u1 fqIF;n g and fqIF;nþ1 g ¼ ð1  u1 ÞfqIF;n g  u1 fqIB;n g where u1 is a relaxation parameter to ensure and/or accelerate convergence for the BEM sub-domain:   B    B  uB;nþ1 qB;nþ1 H11 H12 G11 G12 solve ¼ for fuIB;nþ1 g H21 H22 G21 G22 uIB;nþ1 qIB;nþ1 for the FEM sub-domain:   F   F  uF;nþ1 fF;nþ1 K11 K12 I solve g ¼ ½M fqIF;nþ1 g for fuIF;nþ1 g ¼ and ffF;nþ1 I K21 K22 uIF;nþ1 fF;nþ1 apply fuIB;nþ2 g ¼ ð1  u2 ÞfuIB;nþ1 g þ u2 fuIF;nþ1 g and fuIF;nþ2 g ¼ ð1  u2 ÞfuIF;nþ1 g þ u2 fuIB;nþ1 g where u2 is a relaxation parameter to ensure and/or accelerate convergence.

2984

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

4. BEM–BEM coupling Consider Fig. 2, where the domain of the original problem is decomposed into m BEM sub-domains, X1 ; X2 ; . . . ; Xm . Let SðiÞ be the indices of those sub-domains that are neighbors of sub-domain Xi . Now, let us define the following vectors: uIi : potentials (displacements) in the sub-domain Xi at the interface with all neighboring sub-domains, i.e., IðijÞ uIi ¼ [ui for all j 2 SðiÞ, uBi : non-interface potentials (displacements) in the sub-domain Xi . Similarly, one can define the flux (traction) vectors for the BEM sub-domains. The discretized integral equations for the BEM sub-domains ði ¼ 1; . . . ; mÞ may be written in the following form: 

Hi;11 Hi;21

Hi;12 Hi;22



uBi uIi





Gi;11 ¼ Gi;21

Gi;12 Gi;22



 qBi : qIi

ð10Þ

At the interface ðCkl Þ between two adjacent sub-domains, Xk and Xl , the compatibility and equilibrium conditions should be satisfied, i.e., IðklÞ

g ¼ ful

IðklÞ

g ¼ fql

fuk fqk

IðlkÞ

g;

IðlkÞ

ð11Þ g:

ð12Þ

In Section 4.1 we briefly review the Neumann–Neumann BEM–BEM domain decomposition coupling algorithm presented by Kamiya et al. [14]. In Section 4.2 we present the interface relaxation BEM–BEM coupling algorithms.

Fig. 2. (a) Domain decomposed into m sub-domains. (b) BEM modeling for X1 and X2 .

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

4.1. Neumann–Neumann BEM–BEM coupling algorithm The Neumann–Neumann BEM–BEM coupling algorithm may be described as follows [14]: Set initial guess fqIi;0 g for the BEM sub-domains ði ¼ 1; . . . ; mÞ. Do for n ¼ 0; 1; 2; . . . until convergence for i ¼ 1; . . . ; m     B  uBi;n qi;n Hi;11 Hi;12 Gi;11 Gi;12 solve ¼ for fuIi;n g Hi;21 Hi;22 Gi;21 Gi;22 uIi;n qIi;n IðijÞ

IðijÞ

IðjiÞ

IðijÞ

apply fqi;nþ1 g ¼ fqi;n g þ bðfuj;n g  fui;n gÞ for all j 2 SðiÞ where b is a relaxation parameter to ensure and/or accelerate convergence. 4.2. Interface relaxation BEM–BEM coupling algorithms The geometric contraction BEM–BEM based coupling algorithm may be described as follows: Set initial guess fuIi;0 g for the BEM sub-domains ði ¼ 1; . . . ; mÞ. For n ¼ 0; 1; 2; . . ., do until convergence for i ¼ 1; . . . ; m     B  uBi;n qi;n Hi;11 Hi;12 Gi;11 Gi;12 solve ¼ for fqIi;n g Hi;21 Hi;22 Gi;21 Gi;22 uIi;n qIi;n IðijÞ

IðijÞ

IðijÞ

IðjiÞ

apply fui;nþ1 g ¼ fui;n g  aðfqi;n g þ fqj;n gÞ for all j 2 SðiÞ: The Robin relaxation BEM–BEM coupling algorithm may be described as follows: IðijÞ

IðjiÞ

IðjiÞ

Define fgi;nþ1 g ¼ fqj;n g þ qfuj;n g for all j 2 SðiÞ Set initial guess fuIi;0 g for the BEM sub-domains ði ¼ 1; . . . ; mÞ. For n ¼ 0; 1; 2; . . ., do until convergence for i ¼ 1; . . . ; m with the Robin boundary conditions on the interface   B    B  ui;n qi;n Hi;11 Hi;12 Gi;11 Gi;12 solve ¼ for fuIi;n g and fqIi;n g Hi;21 Hi;22 Gi;21 Gi;22 uIi;n qIi;n IðijÞ

IðjiÞ

IðjiÞ

apply fgi;nþ1 g ¼ fqj;n g þ qfuj;n g for all j 2 SðiÞ The Dirichlet/Neumann averaging coupling algorithm may be described as follows: Set initial guess fuIi;0 g for the BEM sub-domains ði ¼ 1; . . . ; mÞ. For n ¼ 0; 2; 4; . . ., do until convergence for i ¼ 1; . . . ; m     B  uBi;n qi;n Hi;11 Hi;12 Gi;11 Gi;12 solve ¼ for fqIi;n g I Hi;21 Hi;22 Gi;21 Gi;22 ui;n qIi;n IðijÞ

IðijÞ

IðijÞ

IðijÞ

IðjiÞ

apply fqi;nþ1 g ¼ ð1  u1 Þfqi;n g  u1 fqj;n g for all j 2 SðiÞ   B    B  ui;nþ1 qi;nþ1 Hi;11 Hi;12 Gi;11 Gi;12 solve ¼ for fuIi;nþ1 g I Hi;21 Hi;22 G G ui;nþ1 qIi;nþ1 i;21 i;22 IðjiÞ

apply fui;nþ2 g ¼ ð1  u2 Þfui;nþ1 g þ u2 fuj;nþ1 g for all j 2 SðiÞ

2985

2986

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

5. Numerical examples As a first FEM–BEM coupling example, we consider a potential flow problem (Fig. 3). The domain of the original problem is decomposed into XB and XF which are governed by Laplace equation, i.e., ki r2 u ¼ 0 in Xi , where ki is the material property in the sub-domain Xi . The boundary conditions are selected such that uð0; yÞ ¼ 0, uða; yÞ ¼ 200 and zero flux elsewhere. Utilizing the interface relaxation algorithms presented in Section 3, the problem is investigated using different values of aB =aF and kB =kF . For aB =aF ¼ 1, the domain is modeled by 18 linear boundary elements and 40 linear triangular elements (see, i.e., Fig. 3 shows the discretization for aB =aF ¼ 1). The results are in good agreement with the analytic solution. Fig. 4 shows the applicable range and optimal values of the relaxation parameters (aB =aF ¼ 1 and kB =kF ¼ 2). Note that for the Dirichlet/Neumann averaging algorithm we set u1 ¼ u2 ¼ u. Beyond the values shown in Fig. 4, the interface relaxation coupling algorithms will not converge. For the Robin interface relaxation FEM–BEM coupling algorithm the relaxation parameter might be assigned the value q ¼ 1:0, as this is determined as the optimal value. With the increase of q, the algorithm is also converging but with a higher number of iterations. Tables 1 and 2 show the applicable range and the optimal values of the relaxation parameters determined experimentally with different combinations of aB =aF and kB =kF . We verified the applicable and optimal values of b, c and h, using the convergence and optimal convergence conditions

Fig. 3. (a) Potential flow problem. (b) FEM–BEM discretization.

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

2987

Fig. 4. Applicable range and optimal values of the relaxation parameters for aB =aF ¼ 1 and kB =kF ¼ 2.

Table 1 Applicable range of the relaxation parameters for different values of aB =aF and kB =kF aB =aF

Relaxation parameter

kB =kF 0.50

1.0

2.0

0.2

Geometric contraction based a Parallel Neumann–Neumann b Parallel Dirichlet–Neumann c Sequential Dirichlet–Neumann h Robin relaxation q

0.02–0.56 0.02–1.40 0.02–0.36 0.02–0.56 0.02–

0.02–0.32 0.02–1.66 0.02–0.18 0.02–0.32 0.02–

0.02–0.18 0.02–1.80 0.02–0.08 0.02–0.18 0.02–

1.0

Geometric contraction based a Parallel Neumann–Neumann b Parallel Dirichlet–Neumann c Sequential Dirichlet–Neumann h Robin relaxation q

0.02–1.32 0.02–0.66 0.02–1.98 0.02–1.32 0.02–

0.02–0.98 0.02–0.98 0.02–0.98 0.02–0.98 0.02–

0.02–0.66 0.02–1.32 0.02–0.48 0.02–0.66 0.02–

4.0

Geometric contraction based a Parallel Neumann–Neumann b Parallel Dirichlet–Neumann c Sequential Dirichlet–Neumann h Robin relaxation q

0.02–1.76 0.02–0.22 0.02-2.28 0.02–1.76 0.02–

0.02–1.56 0.02–0.40 0.02-2.66 0.02–1.56 0.02–

0.02–1.32 0.02–0.66 0.02–1.98 0.02–1.32 0.02–

given in Refs. [11–13]. Fig. 5 shows the applicable range of the relaxation parameters for the Dirichlet/ Neumann averaging algorithm obtained using different combinations of u1 and u2 . As a BEM–BEM coupling example, let us reconsider the potential flow problem but with the domain decomposed into three sub-domains (Fig. 6). The Neumann–Neumann interface relaxation algorithm is not suitable for solving this problem. Utilizing the algorithm leads to the Neumann condition being specified on the whole boundary of the sub-domain X2 . The problem is solved using the geometric contraction based BEM–BEM coupling algorithm and the results are in good agreement with the analytic solution. Fig. 7 shows the applicable ranges and the optimal values of the relaxation parameter a for different combinations

2988

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

Table 2 Optimal values of the relaxation parameters aB =aF

Relaxation parameter

kB =kF 0.50

1.0

2.0

0.2

Geometric contraction based a Parallel Neumann–Neumann b Parallel Dirichlet–Neumann c Sequential Dirichlet–Neumann h

0.28 0.72 0.10 0.28

0.16 0.84 0.06 0.16

0.10 0.92 0.04 0.10

1.0

Geometric contraction based a Parallel Neumann–Neumann b Parallel Dirichlet–Neumann c Sequential Dirichlet–Neumann h

0.66 0.34 0.38 0.66

0.50 0.50 0.18 0.50

0.32 0.68 0.12 0.32

4.0

Geometric contraction based a Parallel Neumann–Neumann b Parallel Dirichlet–Neumann c Sequential Dirichlet–Neumann h

0.88 0.10 0.78 0.88

0.78 0.2 0.54 0.78

0.66 0.34 0.38 0.66

Fig. 5. Applicable range of the relaxation parameters––the Dirichlet/Neumann averaging algorithm for different combinations of u1 and u2 .

of k1 , k2 and k3 . For example, a should be within 0  0:54 to assure convergence and the optimal value is determined as 0.38, for k1 ¼ 1:0, k2 ¼ 1:0 and k3 ¼ 2:0. Now let us consider Fig. 8, where a steel cantilever beam is subjected to a uniform tensile loading of 20  103 units at its free end. The beam is considered to be in a state of plane stress with an elastic modulus, E ¼ 29  106 units, and PoisonÕs ratio m ¼ 0:3. The beam is 20 units long and 10 units high, and is assumed to be weightless. The problem is solved using the sequential Dirichlet–Neumann FEM–BEM and the parallel Dirichlet–Neumann FEM–BEM coupling algorithms with the discretization shown in Fig. 8. The problem is investigated for different relative values of modulus of elasticity for the BEM and FEM sub-domains, EB =EF and the results obtained using the two algorithms match very well with the analytic solutions. Fig. 9 shows the applicable range and optimal value of c and h. The results in Fig. 9 indicate that as EB =EF decreases, the ranges from which the relaxation parameters to be chosen increases. The applicable ranges for both c and h reduce to very narrow ones for higher values of EB =EF .

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

Fig. 6. (a) BEM–BEM coupling example. (b) Domain decomposed into three sub-domains.

Fig. 7. Applicable ranges and the optimal values of the relaxation parameter a for different combinations of k1 , k2 and k3 .

2989

2990

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

Fig. 8. (a) Cantilever beam subjected to uniform tension. (b) Discretization of the beam.

Fig. 9. Applicable ranges and optimal values for the cantilever beam problem.

The problem is reinvestigated using the geometric contraction based coupling algorithm. Again the results are in good agreement with the analytic solution. For EB =EF ¼ 5, the applicable range and the optimal value of the relaxation parameter are determined as 0  0:12 and 0.11, respectively. Note that fug and fqg are scaled while choosing a.

6. Selection of the relaxation parameters Sections 3 and 4 present a wide class of interface relaxation FEM–BEM and BEM–BEM coupling algorithms (decompose and relax the interface values). Some of them use ‘‘history’’ (the new value at the

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

2991

interface is explicitly set to be the old one plus a correction term) some do not. The algorithms exhibit different mathematical and computational properties. The speed of convergence of the interface relaxation methods can be of high, moderate or low. It is natural to expect that the rate of convergence of all the interface relaxation algorithms be affected by certain problem parameters. Some of most important of these parameters are the mesh or grid discretization and geometric and material characteristics of the FEM and BEM sub-domains (BEM sub-domains in case of BEM–BEM coupling). The convergence and optimal convergence of the interface relaxation coupling algorithms are ensured by properly set relaxation parameters, which may be obtained by experimenting with different values. Alternatively, and prior to calculations, static values of the relaxation parameters may be determined. Elleithy et al. [11,12] investigated the convergence of the sequential Dirichlet–Neumann FEM–BEM domain decomposition coupling algorithm. El-Gebeily et al. established the convergence and optimal convergence conditions for the parallel Dirichlet–Neumann and Neumann–Neumann FEM–BEM domain decomposition coupling algorithms [13]. The problem is still open for establishing the convergence conditions and for the selection of the relaxation parameters. For the geometric contraction based and Robin interface relaxation coupling algorithms, the convergence conditions are to be established. Further, an extensive and systematic performance evaluation study of all interface relaxation algorithms may be considered for future work. For all the interface relaxation algorithms presented in Sections 3 and 4, the initial guess of the variables on the interface (interfaces) has an insignificant effect on the speed of convergence. It is, i.e., reasonable to start with values of zeros for the initial potentials (displacements) on the interface, which seems convenient, as well as, appropriate from the physical realization.

7. Conclusions In this paper, we present several interface relaxation algorithms for FEM–BEM coupling and BEM– BEM coupling (decompose and relax the interface values). The interface relaxation coupling algorithms make it easy to utilize different software for different sub-domains. The algorithms exhibit different mathematical and computational properties. The speed of convergence of the interface relaxation methods can be of high, moderate or low depending on the particular application. The rate of convergence of all the interface relaxation algorithms is affected by certain problem parameters. Some of most important of these parameters are the geometric and material characteristics of the FEM and BEM sub-domains (BEM subdomains in case of BEM–BEM coupling). The interface relaxation coupling algorithms presented in this paper, have the potential to work effectively, however, there is still much to be learned about their behavior and how to choose among them or to choose their parameters.

Acknowledgement The authors gratefully acknowledge the support of the Japan Society for the Promotion of Science and Shinshu University, Japan.

References [1] O.C. Zienkiewicz, D.W. Kelly, P. Bettes, The coupling of the finite element method and boundary solution procedures, Int. J. Numer. Methods Engrg. 11 (1977) 355–375. [2] State of the art in BEM/FEM coupling, Int. J. Boundary Element Methods Commun. 4 (2) (1993) 58–67.

2992

W.M. Elleithy, M. Tanaka / Comput. Methods Appl. Mech. Engrg. 192 (2003) 2977–2992

[3] State of the art in BEM/FEM coupling, Int. J. Boundary Element Methods Commun. 4 (3) (1993) 94–104. [4] S. Ganguly, J.B. Layton, C. Balakrishma, Symmetric coupling of multi-zone curved Galerkin boundary elements with finite elements in elasticity, Int. J. Numer. Methods Engrg. 48 (2000) 633–654. [5] R.A. Bialecki, Z. Ostrowski, A.J. Kassab, Q. Yin, E. Sciubba, Coupling BEM, FEM and analytic solutions in steady-state potential problems, Engrg. Anal. Boundary Elem. 26 (2002) 597–611. [6] W.M. Elleithy, H.J. Al-Gahtani, An overlapping domain decomposition approach for coupling the finite and boundary element methods, Engrg. Anal. Boundary Elem. 24 (5) (2000) 391–398. [7] E. Stein, M. Kreienmeyer, Coupling of BEM and FEM by a multiplicative Schwarz method and its parallel implementation, Engrg. Comput. 15 (2) (1998) 173–198. [8] N. Kamiya, H. Iwase, E. Kita, Parallel computing for the combination method of BEM and FEM, Engrg. Anal. Boundary Elem. 18 (1996) 221–229. [9] C.-C. Lin, E.C. Lawton, J.A. Caliendo, L.R. Anderson, An iterative finite element–boundary element algorithm, Comput. Struct. 39 (5) (1996) 899–909. [10] Y.T. Feng, D.R.J. Owen, Iterative solution of coupled FE/BE discretization for plate-foundation interaction problems, Int. J. Numer. Methods Engrg. 39 (1996) 1889–1901. [11] W.M. Elleithy, H.J. Al-Gahtani, M. El-Gebeily, Convergence of the iterative coupling of BEM and FEM, in: 21st World Conference on the Boundary Element Method, BEM21, Oxford University, UK, August 1999, pp. 281–290. [12] W.M. Elleithy, H.J. Al-Gahtani, M. El-Gebeily, Iterative coupling of BE and FE methods in elastostatics, Engrg. Anal. Boundary Elem. 25 (8) (2001) 685–695. [13] M. El-Gebeily, W.M. Elleithy, H.J. Al-Gahtani, Convergence of the domain decomposition finite element–boundary element coupling methods, Comput. Methods Appl. Mech. Engrg. 191 (43) (2002) 4851–4867. [14] N. Kamiya, H. Iwase, E. Kita, Parallel implementation of boundary element method with domain decomposition, Engrg. Anal. Boundary Elem. 18 (1997) 209–216. [15] M. Mu, Solving composite problems with interface relaxation, SIAM J. Sci. Comput. 20 (4) (1999) 1394–1416. [16] D. Funaro, A. Quarteroni, P. Zanolli, An interface procedure with interface relaxation for domain decomposition methods, SIAM J. Numer. Anal. 25 (6) (1988) 1213–1236. [17] P.L. Lions, On the Schwarz alternating method: a variant for nonoverlapping subdomains, in: R. Glowinski, G.H. Golub, G.A. Meurant, J. Periaux (Eds.), Domain Decomposition Methods for Partial Differential Equations, SIAM, 1990, pp. 202–223. [18] M. Mu, J.R. Rice, Modeling with collaborating PDE solvers––theory and practice, Comput. Syst. Engrg. 6 (1995) 87–95. [19] J.R. Rice, P. Tsompanopoulou, E.A. Vavalis, Interface relaxation methods for elliptic differential equations, Appl. Numer. Math. 32 (1999) 219–245.