Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity – An SGBEM Approach for Nonconforming Discretizations of Curved Interfaces

Copyright c 2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity – An SGBEM Approach for Nonconforming Discretizations of Curved Interfaces R. Vodiˇ cka1,V.Manti ˇ c2and F. París2 Abstract: An original approach to solve domain decomposition problems by the symmetric Galerkin boundary element method is developed. The approach, based on a new variational principle for such problems, yields a fully symmetric system of equations. A natural property of the proposed approach is its capability to deal with nonconforming discretizations along straight and curved interfaces, allowing in this way an independent meshing of non-overlapping subdomains to be performed. Weak coupling conditions of equilibrium and compatibility at an interface are obtained from the critical point conditions of the energy functional. Equilibrium is imposed through local traction (Neumann) boundary conditions prescribed on a subdomain situated at one side of the interface, and compatibility is imposed through local displacement (Dirichlet) boundary conditions prescribed on the other subdomain situated at the opposite side of the interface. No additional unknowns such as Lagrange multipliers are introduced. An SGBEM code for 2D elastic domain decomposition problems has been implemented. The effectiveness of the approach developed is documented by numerical examples involving non-matching linear boundary element meshes at the interfaces, where the accuracy is analyzed by comparing the numerical results obtained versus the analytical solutions and by evaluating the convergence rate of the error in the (discretized) integral L2-norm and maximum-norm for h-refinements of boundary el- 1Technical University of Košice, Faculty of Civil Engineering, Department of Mathematics, Vysokoškolská 4, 042 00 Košice, Slovakia, on leave at2 2University of Seville, School of Engineering, Group of Elasticity and Strength of Materials, Camino de los Descubrimientos s/n, 41092 Seville, Spain ement meshes. Finally, a theoretical analysis of a problem with an interior and an exterior subdomain is introduced to explain the observed behaviour of numerical results. Keyword: boundary integral equations, boundary elements, SGBEM, domain decomposition, variational formulation, energy functional, nonmatching meshes. 1 Introduction In the solution of Boundary Value Problems (BVPs) of continuum mechanics it is sometimes useful or necessary to split the domain into two or more subdomains, such a procedure being referred to hereinafter as Domain Decomposition (DD). DD procedures are usually employed in solving multi-material and contact BVPs or BVPs including complex structures, sometimes with independent modeling of substructures, and in developing algorithms to solve large-scale BVPs on modern parallel computers. DD procedures are also used in coupling different discretization schemes, e.g. Finite Element Method (FEM) and Boundary Element Method (BEM). The development of numerical techniques for solution of BVPs via DD (DDBVPs) has substantially increased recently. There exist several ways of mathematical formulation and solution of DDBVPs [Quarteroni and Valli (1999); Steinbach (2003)]. In FEM, special attention is paid to developing techniques based on the mortar element concept; for recent works on application of this approach to elastic and contact problems see McDevitt and Laursen (2000); Puso (2004) and further references given therein. One of the 174 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 main problems is to find reliable methods tying nonconforming meshes along curved interfaces; in fact most of the works on mortar techniques have been restricted to straight/planar interfaces. The capability of a DD approach to deal with such discretizations represents an important advantage for engineering applications. Similar techniques are also under development in BEM [Blázquez et al. (1998); Hsiao et al. (2000a); González and Abascal (2000); Vodiˇcka (2000); Blázquez et al. (2002)], only the unknowns along the external boundary and along the interfaces of subdomains being directly computed. The application of DD in BEM is in fact required when solving multi-material BVPs. Symmetric Galerkin BEM (SGBEM) [Sirtori (1979); Bonnet et al. (1998)] represents a valuable alternative to the classical collocational BEM as it combines the boundary character of the collocational BEM with the symmetrical and energy based formulation typical for FEM. Several techniques have been developed by different authors [Maier et al. (1991); Gray and Paulino (1997); Ganguly et al. (1999); Vodiˇcka (2000); Panzeca et al. (2002)] for application of SGBEM to DD and contact BVPs. It is expected that SGBEM will also enable an easy and straightforward coupling with FEM [Hsiao et al. (2000b); Han and Atluri (2002); Ganguly et al. (2004); Springhetti et al. (2004)]. Some other advanced applications of SGBEM can be found in Aimi et al. (2003); Qian et al. (2004) and Duddeck (2006). A natural way to formulate SGBEM is a variational formulation, which provides SGBEM with some aspects analogous to FEM, e.g. conditions of a critical point of an energy functional represent the final system of equations to be solved and a convergence of error in energy norms can be relatively easily shown [Wendland and Hsiao (2004)]. Thus, it is worth an effort to deduce a variational formulation of SGBEM for DDBVPs. This has been done in the present work by generalizing to DDBVPs, through the development of an original approach, a variational formulation of one-domain SGBEM due to Bonnet (1995). A variational formulation of a DDBVP requires an adequate formulation of the coupling conditions along the interface. In the present work, it arises from a weak formulation of elastic contact conditionsproposed by Blázquez et al. (1998) (developing an original idea by Schnack (1987)), here adapted for DDBVP with perfect interfaces. In fact, the variational principle proposed produces the weak coupling conditions as a natural generalization of the classical strong point-wise coupling conditions, in such a way that the developed approach shares some ideas with the mortar approaches used in FEM [McDevitt and Laursen (2000); Puso (2004)]. An implementation of the above weak formulation of coupling conditions for non-matching boundary element meshes at curved interfaces has been included in the SGBEM code developed. Various methods of data transfer between nonmatching meshes via integrals computing over the discretized curved surfaces were introduced in Jiao and Heath (2004), one of these methods based on the concept of common refinement having been used in the present code. The structure of the present paper is as follows. Elastic BVP and DDBVP together with the two basic Boundary Integral Equations (BIEs) are introduced in the preliminary Section 2. In Section 3 the main theoretical results are presented. After introducing the functional of the energy for DDBVPs, it is shown that its critical point gives the solution of the original DDBVP. The resulting weak formulation of the BIE system, whose solution corresponds to this critical point, is first deduced in the integral form. Then, an operator notationis used to show easily its overall symmetrical structure. The discretization of the BIE system obtained by SGBEM, using linear continuous boundary elements, is introduced in Section 4. The accuracy and convergence for h-refinements of the present approach are documented in Section 5 by three numerical examples. In order to explain some aspects of the observed behaviour of the numerical results, a theoretical analysis of a simple DDBVP, with a bounded interior and an unbounded exterior subdomain, at the discretized level, is presented in the Appendix, where a reduction of the resulting system matrix to its Schur complement is carried out and the obtained matrix Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 175 operators, which include Galerkin discretizations of the local Steklov-Poincaré operators for subdomains [Quarteroni and Valli (1999); Hsiao et al. (2000b); Steinbach (2003)], are discussed. 2 Preliminaries. 2.1 Elastic BVP Let us consider a linear elastic body defined by a domain Ω⊂Rd(d=2or3)with a bounded Lipschitz boundary ∂ Ω=Γ(i.e. Γgiven locally as graphs of Lipschitz functionsin a finite number of appropriate cartesian coordinate systems, and Ωbeing locally on one side of Γ). Note that Γ may include corners but not cracks and cusps. Let ΓS⊂Γdenote the smooth part of Γ, i.e. excluding corners, edges, points of curvature jumps, etc. Let ndenote the unit outward normal vector defined at ΓS. Let xi(i=1,...,d) be a fixed cartesian coordinate system at the origin point O.Letu= (u1,...,ud)be the displacement solution of the following BVP for the Navier equation with zero body forces: cijkluk,lj (x)=cijkl ε kl,j(u(x)) = 0,x∈Ω,(1a) ui(x)=gi(x),x∈Γu,(1b) ti(x)=(Tn(x))i(u(x)) = hi(x),x∈Γt,(1c) with the fourth-order symmetric positively definite tensor of elastic stiffness cijkl (i,j,k,l= 1,...,d), strain tensor ε ij, traction operator Tn and the split of the boundary Γ=Γu∪Γt(/0= Γu∩Γt) in accordance with the boundary conditions. 2.2 Elastic BIEs The Somigliana displacement identity [Jaswon and Symm (1977); Baláš et al. (1989); París and Cañas (1997)] holds for the solution of the above BVP defined on a bounded homogenous Ω(interior BVP): χ Ω(x)ui(x)=Γ Uij(x,y)tj(y)dyS −Γ Tij(x,y)uj(y)dyS, x∈Rd\Γ,i,j=1,...,d, (2) where χ Ωis the characteristic function of Ω ( χ Ω(x)=1ifx∈Ωand χ Ω(x)=0ifx∈ Ω∪Γ), Uij is the fundamental solution in displacements of the Navier equation and Tij represents the fundamental tractions, obtained from the fundamental solution via the traction operator: T(x,y)= Tn(y)U(x,y)T,Tdenoting the transpose matrix. Note that (2) is also valid for a solution of an exterior BVP defined on an unbounded Ωwith a bounded Γ[Jaswon and Symm (1977); Baláš et al. (1989)] if the solution fulfills the well-known regularity condition at infinity: ui(x)=Uij(x,0)bj+O(x1−d), x→∞,bj=Γ tj(y)dyS.(3) Applying the traction operator to (2) the Somigliana traction identity is obtained: χ Ω(x)ti(x)=Γ T∗ ij(x,y)tj(y)dyS −Γ Sij(x,y)uj(y)dyS, x∈Rd\Γ, (4) where T∗(x,y)=Tn(x)U(x,y),S(x,y)= Tn(x)T(x,y)and n(x)is the unit vector normal to an auxiliary curve passing through the point x where the traction is to be evaluated. With reference to the above defined integral kernels, let us recall the following reciprocity relations [Bonnet et al. (1998)]: U(x,y)=UT(y,x), T∗(x,y)=TT(y,x)and S(x,y)=ST(y,x). Additionally, the fundamental solution is symmetric, i.e. U(x,y)=UT(x,y), and, according to Mantiˇc and París (1997a,b), also the hypersingular integral kernel Sis symmetric in 2D, i.e. S(x,y)= ST(x,y). Additionally, the kernels T∗(x,y)and T(x,y)change their sign when evaluated using the normals defined with respect to the opposite sides of a curve. After an asymptotic procedure applied to (2) and (4) for a boundary point x∈ΓS, the following BIEs, sometimes called u-BIE and t-BIE respectively, are obtained [Baláš et al. (1989); Guiggiani et al. (1992); Guiggiani (1995); MantiˇcandParís 176 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 (1995); Young (1996); París and Cañas (1997)]: 1 2ui(x)=Γ Uij(x,y)tj(y)dyS −− Γ Tij(x,y)uj(y)dyS,(5a) 1 2ti(x)=− Γ T∗ ij(x,y)tj(y)dyS −= Γ Sij(x,y)uj(y)dyS,(5b) x∈ΓS,i,j=1,...,d. the second and first integrals respectively on the right-hand sides of (5a) and (5b), with a strongly singular integral kernel, being evaluated in the sense of Cauchy principal value, and the second integral on the right-hand side of (5b), with a hypersingular integral kernel, being evaluated in the sense of Hadamard finite part. Both aforementioned BIEs, (5a) and (5b), will be used in the following derivation of the SGBEM formulation for a DDBVP. 2.3 Domain decomposition A DD approach starts with a split of Ωinto several subdomains. For the sake of simplicity, a split into two non-overlapping subdomains ΩA and ΩB, whose respective boundaries are denoted as ΓA= ∂ ΩAand ΓB= ∂ ΩB, will be considered hereinafter, Figure 1. The common part of ΓAand ΓBis denoted as Γc=ΓA∩ΓB. Hence, reconsidering the boundary conditions (1b) and (1c) and the corresponding split of Γ, we can write Γ η = Γ η u∪Γ η t∪Γc,where η =A,Band /0=Γ η u∩Γ η t= Γ η u∩Γc=Γ η t∩Γc. The DDBVP for the Navier equation (1) can be written in the form: c η ijklu η k,lj (x)=c η ijkl ε kl,j(u η (x)) = 0,(6a) x∈Ω η , η =A,B, u η i(x)=g η i(x),x∈Γ η u,(6b) t η i(x)=(Tn η (x))i(u η (x)) = h η i(x),x∈Γ η t,(6c) uA i(x)=uB i(x),x∈Γc,(6d) tA i(x)= −tB i(x),x∈Γc,(6e) where the coupling conditions between ΩAand ΩB, equations (6d) and (6e) respectively, enforce Ox1 x2 ΩAΩB ΓB u ΓA u ΓB t ΓA t Γc gA gB hA hB Figure 1: Domain decomposition problem. compatibility of the displacements uAand uBdefined on subdomains ΩAand ΩBand equilibrium of the tractions tAand tBdefined on ΓAand ΓB. These coupling conditions will be formulated in a weak form in the variational formulation introduced in the next section. 3 Variational formulation 3.1 Energy functional for a DDBVP and its variation Let us introduce an energy functional E(uA,uB)as a function of displacements uAand uB. The functional Ecan be expressed in the following form: E(uA,uB)=EA p(uA)+EB p(uB)+Ec(uA,uB),(7) where the functionals EA p(uA)and EB p(uB)represent the total energy functionals [Hartmann (1985)] associated to ΩAand ΩB, respectively, without the parts belonging to Γc. These functionals can be expressed in the form: E η p(u η )=1 2Ω ηε ij(u η )c η ijkl ε kl(u η )dV −Γ η t h η iu η idS−Γ η u t η iu η i−g η idS, η =A,B, (8) where the fufnctions t η represent the tractions of the displacement solutions u η calculated via the traction operator Tn,i.e.t η =Tn η (u η ). The last term in (7), which defines a form of interface energy associated to coupling conditions Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 177 across Γc, takes the form: Ec(uA,uB)=−Γc tA iuA i−uB idS,(9) which indicates that the solution of the DDBVP will cause this term, called interface constraint functional, to vanish. The condition of the critical point of the functional Ein (7) means that the first order variation of Evanishes, which can be written after a rearrangement as: 0= δ E(uA,uB; δ uA, δ uB)(10a) =ΩA ε ij( δ uA)cA ijkl ε kl(uA)dV−ΓA t hA i δ uA idS −ΓA u tA i δ uA idS−Γc tA i δ uA idS(10b) +ΩB ε ij( δ uB)cB ijkl ε kl(uB)dV−ΓB t hB i δ uB idS −ΓB u tB i δ uB idS−Γc tB i δ uB idS(10c) −ΓA u δ tA iuA i−gA idS+ΓB u δ tB iuB i−gB idS(10d) −Γc δ tA iuA i−uB idS−Γc δ uB itA i+tB idS. (10e) As will be explained in what follows, relation (10) implies that the solution of the DDBVP (6) is the critical point of the functional E.Ifthevolume integrals in expressions (10b) and (10c) are integrated by parts, giving the volume integrals with the virtual displacements δ u η without derivatives, these expressions render the Navier equations (6a) for functions uAand uBand the boundary conditions for tractions (6c). Expressions (10d) vanish only if the boundary conditions for displacements (6b) are fulfilled. Finally, the line (10e) enforces the equilibrium (6e) and compatibility (6d) conditions at Γc in a weak sense: the compatibility is imposed through the displacement (Dirichlet) boundary conditions on subdomain ΩAalong Γcwith ‘prescribed’ displacements uB, and the equilibrium is imposed through the traction (Neumann) boundary conditions on subdomain ΩBalong Γcwith ‘prescribed’ tractions tA. The relations between the functions defined along the boundary Γcare schematically presented in Figure 2. An arrow connects two functions which are equal to each other in the weak sense due to the relation whose number is written in the parentheses close to the arrow. We conclude that finding the critical point of the energy functional E(7) is equivalent to solving the DDBVP (6). Ox1 x2 ΩAΩB ΓB u ΓA u ΓB t ΓA t ΓcΓc gA gB hA hB uB,tB=TnB(uB) uB uA,tA=TnA(uA) tA (10e) (10e) Figure 2: Explanation of the variables in the functional E. In order to obtain a relation suitable for an SGBEM approach, equation (10) is modified integrating by parts the volume integrals, but now yielding displacements u η without derivatives. Recalling the meaning of δ t η as δ t η = Tn η ( δ u η ), the following relation is obtained: 0=−ΩA ∂ ∂ xl ε ij( δ uA)cA ijkluA kdV +ΓA t uA i δ tA idS+ΓA u gA i δ tA idS+Γc uA i δ tA idS −ΓA t hA i δ uA idS−ΓA u tA i δ uA idS−Γc tA i δ uA idS +−ΩB ∂ ∂ xl ε ij( δ uB)cB ijkluB kdV +ΓB t uB i δ tB idS+ΓB u gB i δ tB idS+Γc uB i δ tB idS −ΓB t hB i δ uB idS−ΓB u tB i δ uB idS−Γc tB i δ uB idS −Γc δ tA iuA i−uB idS−Γc δ uB itA i+tB idS. (11) 178 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 3.2 Symmetric variational formulation of BIEs for a DDBVP So far, no restrictions on the virtual displacements δ u η have been set. Let us consider only those δ u η which satisfy the Navier equation in Ω η , causing the volume integrals in (11) to vanish. Then, these displacements have the following boundary integral representation: δ u η i(x)=Γ η U η ij(x,y) ϕη j(y)dyS −Γ η T η ij (x,y) ψη j(y)dyS, x∈Ω η , (12) where the functions ϕη and ψη can be considered to have the form: ϕη = δ t η − δ ˜ t η , ψη = δ u η − δ ˜u η .(13) Note that the representation formula (12) can be obtained from the Somigliana identity (2) written for Ω η and for each connected subset of its complement in Rd. In fact, the function δ ˜u η represents the trace of a displacement solution of the Navier equation at each such connected subset, equipped with the radiation condition (3) if the pertinent subset is unbounded, and δ ˜ t η represents the corresponding tractions. The boundary conditions for these functions can be chosen arbitrarily. They are chosen here satisfying the following conditions: ϕη (x)=0,x∈Γ η t, η =A,B, ϕ A(x)= δ tA(x),x∈Γc,(14a) ψη (x)=0,x∈Γ η u, η =A,B, ψ B(x)= δ uB(x),x∈Γc,(14b) so that the boundary conditions for δ ˜ t η and δ ˜u η read: δ ˜ t η (x)= δ t η (x),x∈Γ η t, η =A,B, δ ˜ tA(x)=0,x∈Γc,(15a) δ ˜u η (x)= δ u η (x),x∈Γ η u, η =A,B, δ ˜uB(x)=0,x∈Γc.(15b) The integral representation (12) can be written with this particular choice of the functions ϕη and ψη in the form: δ u η i(x)=Γ η u U η ij(x,y) ϕη j(y)dyS −Γ η t T η ij(x,y) ψη j(y)dyS +Γc U η ij(x,y) ϕη j(y)dyS −Γc T η ij(x,y) ψη j(y)dyS, x∈Ω η , (16) which for the limit x→x∈Γ η (or directly utilizing the pertinent difference of equations (5a) for Ω η and the connected subsets of its complement in Rd) renders the relations: δ u η i(x)=Γ η u U η ij(x,y) ϕη j(y)dyS −Γ η t T η ij(x,y) ψη j(y)dyS +Γc U η ij(x,y) ϕη j(y)dyS −Γc T η ij(x,y) ψη j(y)dyS, x∈Γ η u, (17a) δ u η i(x)=Γ η u U η ij(x,y) ϕη j(y)dyS −−1 2 ψη i(x)+− Γ η t T η ij(x,y) ψη j(y)dyS +Γc U η ij(x,y) ϕη j(y)dyS −Γc T η ij(x,y) ψη j(y)dyS, x∈(Γ η t)S, (17b) δ u η i(x)=Γ η u U η ij(x,y) ϕη j(y)dyS −Γ η t T η ij(x,y) ψη j(y)dyS +Γc U η ij(x,y) ϕη j(y)dyS −−1 2 ψη i(x)+− Γc T η ij(x,y) ψη j(y)dyS, x∈(Γc)S. (17c) Analogously,the boundary integral representation of the virtual tractions δ t η (see also equation (5b)) Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 179 can be obtained: δ t η i(x)=1 2 ϕη i(x)+− Γ η u T η ∗ ij (x,y) ϕη j(y)dyS −Γ η t S η ij(x,y) ψη j(y)dyS +Γc T η ∗ ij (x,y) ϕη j(y)dyS −Γc S η ij(x,y) ψη j(y)dyS, x∈(Γ η u)S, (18a) δ t η i(x)=Γ η u T η ∗ ij (x,y) ϕη j(y)dyS −= Γ η t S η ij(x,y) ψη j(y)dyS +Γc T η ∗ ij (x,y) ϕη j(y)dyS −Γc S η ij(x,y) ψη j(y)dyS, x∈(Γ η t)S, (18b) δ t η i(x)=Γ η u T η ∗ ij (x,y) ϕη j(y)dyS −Γ η t S η ij(x,y) ψη j(y)dyS +1 2 ϕη i(x)+− Γc T η ∗ ij (x,y) ϕη j(y)dyS −= Γc S η ij(x,y) ψη j(y)dyS, x∈(Γc)S. (18c) Let the representations (17) and (18) be substituted into (11). Recall that the volume integrals in (11) vanish and the variations δ tAand δ uB along Γcare equal to ϕ Aand ψ B, respectively, due to (14). Hence, 0=ΓA t uA i(x)ΓA u TA∗ ij (x,y) ϕ A j(y)dyS −= ΓA t SA ij(x,y) ψ A j(y)dyS+Γc TA∗ ij (x,y) ϕ A j(y)dyS −Γc SA ij(x,y) ψ A j(y)dySdxS +ΓA u gA i(x)1 2 ϕ A i(x)+− ΓA u TA∗ ij (x,y) ϕ A j(y)dyS −ΓA t SA ij(x,y) ψ A j(y)dyS+Γc TA∗ ij (x,y) ϕ A j(y)dyS −Γc SA ij(x,y) ψ A j(y)dySdxS +Γc uA i(x)ΓA u TA∗ ij (x,y) ϕ A j(y)dyS −ΓA t SA ij(x,y) ψ A j(y)dyS +1 2 ϕ A i(x)+− Γc TA∗ ij (x,y) ϕ A j(y)dyS −= Γc SA ij(x,y) ψ A j(y)dySdxS −ΓA t hA i(x)ΓA u UA ij(x,y) ϕ A j(y)dyS −−1 2 ψ A i(x)+− ΓA t TA ij(x,y) ψ A j(y)dyS +Γc UA ij(x,y) ϕ A j(y)dyS −Γc TA ij(x,y) ψ A j(y)dySdxS −ΓA u tA i(x)ΓA u UA ij(x,y) ϕ A j(y)dyS −ΓA t TA ij(x,y) ψ A j(y)dyS+Γc UA ij(x,y) ϕ A j(y)dyS −Γc TA ij(x,y) ψ A j(y)dySdxS −Γc tA i(x)ΓA u UA ij(x,y) ϕ A j(y)dyS −ΓA t TA ij(x,y) ψ A j(y)dyS+Γc UA ij(x,y) ϕ A j(y)dyS −−1 2 ψ A i(x)+− Γc TA ij(x,y) ψ A j(y)dySdxS +ΓB t uB i(x)ΓB u TB∗ ij (x,y) ϕ B j(y)dyS −= ΓB t SB ij(x,y) ψ B j(y)dyS+Γc TB∗ ij (x,y) ϕ B j(y)dyS −Γc SB ij(x,y) ψ B j(y)dySdxS +ΓB u gB i(x)1 2 ϕ B i(x)+− ΓB u TB∗ ij (x,y) ϕ B j(y)dyS −ΓB t SB ij(x,y) ψ B j(y)dyS+Γc TB∗ ij (x,y) ϕ B j(y)dyS −Γc SB ij(x,y) ψ B j(y)dySdxS +Γc uB i(x)ΓB u TB∗ ij (x,y) ϕ B j(y)dyS −ΓB t SB ij(x,y) ψ B j(y)dyS +1 2 ϕ B i(x)+− Γc TB∗ ij (x,y) ϕ B j(y)dyS 180 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 −= Γc SB ij(x,y) ψ B j(y)dySdxS −ΓB t hB i(x)ΓB u UB ij(x,y) ϕ B j(y)dyS −−1 2 ψ B i(x)+− ΓB t TB ij(x,y) ψ B j(y)dyS +Γc UB ij(x,y) ϕ B j(y)dyS −Γc TB ij(x,y) ψ B j(y)dySdxS −ΓB u tB i(x)ΓB u UB ij(x,y) ϕ B j(y)dyS −ΓB t TB ij(x,y) ψ B j(y)dyS+Γc UB ij(x,y) ϕ B j(y)dyS −Γc TB ij(x,y) ψ B j(y)dySdxS −Γc tB i(x)ΓB u UB ij(x,y) ϕ B j(y)dyS −ΓB t TB ij(x,y) ψ B j(y)dyS+Γc UB ij(x,y) ϕ B j(y)dyS −−1 2 ψ B i(x)+− Γc TB ij(x,y) ψ B j(y)dySdxS −Γc ϕ A i(y)uA i(y)−uB i(y)dyS +Γc ψ B i(y)tA i(y)+tB i(y)dyS. (19) Finally, to obtain a variational formulation of the BIE system for SGBEM solution of DDBVPs, first the order of the integration in the double integrals is interchanged, and then the terms in the expression obtained are reordered to gather together the terms with the same function ϕη or ψη , utilizing the above mentioned reciprocity properties of the integral operators. The interchangeability of the order of the Cauchy principal value and Hadamard finite part integrals with the common integral has been assumed due to the results of the analysis presented in Bonnet (1995). Hence, 0=ΓA u ϕ A j(y)−ΓA u UA ji(y,x)tA i(x)dxS +ΓA t TA ji(y,x)uA i(x)dxS−Γc UA ji(y,x)tA i(x)dxS +Γc TA ji(y,x)uA i(x)dxS−ΓA t UA ji(y,x)hA i(x)dxS +1 2gA j(y)+− ΓA u TA ji(y,x)gA i(x)dxSdyS +ΓA t ψ A j(y)ΓA u TA∗ ji (y,x)tA i(x)dxS −= ΓA t SA ji(y,x)uA i(x)dxS+Γc TA∗ ji (y,x)tA i(x)dxS −Γc SA ji(y,x)uA i(x)dxS +−1 2hA j(y)+− ΓA t TA∗ ji (y,x)hA i(x)dxS −ΓA u SA ji(y,x)gA i(x)dxSdyS +Γc ϕ A j(y)uB j(y)−ΓA u UA ji(y,x)tA i(x)dxS +ΓA t TA ji(y,x)uA i(x)dxS−Γc UA ji(y,x)tA i(x)dxS +−1 2uA j(y)+− Γc TA ji(y,x)uA i(x)dxS −ΓA t UA ji(y,x)hA i(x)dxS+ΓA u TA ji(y,x)gA i(x)dxSdyS +Γc ψ A j(y)ΓA u TA∗ ji (y,x)tA i(x)dxS −ΓA t SA ji(y,x)uA i(x)dxS +−1 2tA j(y)+− Γc TA∗ ji (y,x)tA i(x)dxS −= Γc SA ji(y,x)uA i(x)dxS+ΓA t TA∗ ji (y,x)hA i(x)dxS −ΓA u SA ji(y,x)gA i(x)dxSdyS +ΓB u ϕ B j(y)−ΓB u UB ji(y,x)tB i(x)dxS +ΓB t TB ji(y,x)uB i(x)dxS−Γc UB ji(y,x)tB i(x)dxS +Γc TB ji(y,x)uB i(x)dxS−ΓB t UB ji(y,x)hB i(x)dxS +1 2gB j(y)+− ΓB u TB ji(y,x)gB i(x)dxSdyS +ΓB t ψ B j(y)ΓB u TB∗ ji (y,x)tB i(x)dxS −= ΓB t SB ji(y,x)uB i(x)dxS+Γc TB∗ ji (y,x)tB i(x)dxS −Γc SB ji(y,x)uB i(x)dxS +−1 2hB j(y)+− ΓB t TB∗ ji (y,x)hB i(x)dxS −ΓB u SB ji(y,x)gB i(x)dxSdyS +Γc ϕ B j(y)−ΓB u UB ji(y,x)tB i(x)dxS Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 181 +ΓB t TB ji(y,x)uB i(x)dxS−Γc UB ji(y,x)tB i(x)dxS +1 2uB j(y)+− Γc TB ji(y,x)uB i(x)dxS −ΓB t UB ji(y,x)hB i(x)dxS +ΓB u TB ji(y,x)gB i(x)dxSdyS +Γc ψ B j(y)tA j(y)+ΓB u TB∗ ji (y,x)tB i(x)dxS −ΓB t SB ji(y,x)uB i(x)dxS +1 2tB j(y)+− Γc TB∗ ji (y,x)tB i(x)dxS −= Γc SB ji(y,x)uB i(x)dxS+ΓB t TB∗ ji (y,x)hB i(x)dxS −ΓB u SB ji(y,x)gB i(x)dxSdyS. (20) Introducing the following operator notation will allow us to rewrite the variational formulation of the BIE system (20) in a compact and transparent form. Let ωη T rZ η rsw η s=Γ η r ωη j(y)Γ η s Z η ji(y,x)w η i(x)dxSdyS, (21) where ω stands for ϕ or ψ ,wstands for uor t,r and sstand for u,tor c,andZstands for U,T, T∗or S, and where the inner integral can be regular, weakly singular, Cauchy principal value or Hadamard finite part integral. Hereinafter, u η r= u η |Γ η rand t η r=t η |Γ η rfor r=u,tand u η c=u η |Γc and t η c=t η |Γcfor η =A,Bmean the restrictions of displacements and tractions on the respective boundary parts. Then (20) reads: 0= ϕ AT u−UA uutA u+TA utuA t−UA uctA c+TA ucuA c −UA uthA+1 2IA uu +TA uugA + ψ AT tTA∗ tu tA u−SA tt uA t+TA∗ tc tA c−SA tcuA c +−1 2IA tt +TA∗ tt hA−SA uugA + ϕ AT c−UA cutA u+TA ct uA t−UA cctA c +−1 2IA cc +TA ccuA c +IAB cc uB c−UA cthA+TA cugA + ψ AT cTA∗ cu tA u−SA ctuA t+−1 2IA cc +TA∗ cc tA c −SA ccuA c+TA∗ ct hA−SA cugA + ϕ BT u−UB uutB u+TB ut uB t−UB uctB c+TB ucuB c −UB ut hB+1 2IB uu +TB uugB + ψ BT tTB∗ tu tB u−SB ttuB t+TB∗ tc tB c−SB tcuB c +−1 2IB tt +TB∗ tt hB−SB uugB + ϕ BT c−UB cutB u+TB ctuB t−UB cctB c +1 2IB cc +TB ccuB c−UB ct hB+TB cugB + ψ BT cTB∗ cu tB u−SB ctuB t+1 2IB cc +TB∗ cc tB c −SB ccuB c+IBA cc tA c+TB∗ ct hB−SB cugB (22) or equivalently in matrix-operator form: ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ϕ A u ψ A t ϕ A c ψ A c ϕ B u ψ B t ϕ B c ψ B c ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ T MATRIX I ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ tA u uA t tA c uA c tB u uB t tB c uB c ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ϕ A u ψ A t ϕ A c ψ A c ϕ B u ψ B t ϕ B c ψ B c ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ T MATRIX II ⎛ ⎜ ⎜ ⎝ hA gA hB gB ⎞ ⎟ ⎟ ⎠ , (23) 188 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 -4.10-3 -3.10-3 -2.10-3 -1.10-3 0.100 1.10-3 2.10-3 0 0.2 0.4 0.6 0.8 1 -2.10-3 -1.10-3 0.100 1.10-3 2.10-3 3.10-3 4.10-3 (tn)1-(ta)1 (tn)2-(ta)2 s t1 t2 t1 t2 NA<NB (B) NA<NB (A) NA>NB (B) NA=NB (B) -1.10-2 -5.10-3 0.100 5.10-3 1.10-2 0 0.2 0.4 0.6 0.8 1 -1.10-2 -5.10-3 0.100 5.10-3 1.10-2 (tn)1-(ta)1 (tn)2-(ta)2 s t1 t2 NA<NB (B) NA>NB (A) Figure 9: Rectangular plate under exponential load. Traction errors, meshes 16–40, all meshes (top), finer meshes (bottom). As already mentioned, the reasons for the above observed strong oscillatory behaviour of tractions at the finer mesh when NA>NBare discussed, studying a simple DDBVP in more detail and analyzing the numerical procedure of the solution of the linear equation system, in Appendix and especially in part Appendix A.3 where an explanation of these oscillations is given. 5.3.2 Ring under circumferential shear Numerical results are presented for meshes 12– 28 along Γc. The displacements obtained at the finer meshes of the nonconforming discretizations shown in Figure 10 approximate very well the analytic solution. Due to the fact that the displacement components u1and u2are mutually symmetric with respect to the angle α = π 4(u1( α )= u2( π 2− α )), only results for u1will be presented in what follows. 0 0.5 1 1.5 2 2.5 3 3.5 0π/12 π/6 π/4 π/3 5π/12 π/2 0 0.5 1 1.5 2 2.5 3 3.5 u* 1 u* 2 α u1u2 u1 u2 NA<NB (B) NA>NB (A) analytical Figure 10: Ring under circumferential shear. Displacement distributions, finer meshes, meshes 12– 28. In Figure 11, where displacement errors are plotted, a strong oscillatory behaviour of the error obtained at the finer mesh in the case NA>NB is again observed, in contrast with a relatively smooth distribution of the error at other meshes. Nevertheless, in much the same way as was already observed in the previous example, these oscillations are rather regular. The pattern of the nonconforming discretization is 3–7 in the present case, each third node of the coarser mesh coinciding with each seventh node of the finer mesh. The same pattern is repeated in the displacement errors as can be seen in Figure 11 looking at the triangle marks, filled and unfilled. The explanation given in the previous example, see Figure 6, can be applied here as well. Let the consistency error, which represents the accuracy of the approximation of the interface conditions by nonconforming discretizations, be defined as the difference of the solutions, or equivalently of their absolute errors, at the adjacent sides of Γc. It is interesting to observe in Figure 11 that the consistency error of displacements of a similar magnitude as in the case NA>NBappears also in the case NA<NB. Nevertheless, in the latter case this consistency error is rather smooth in comparison with the former case. Notice that the maximum value of the absolute error is achieved at the coarser mesh in the case NA<NBand that the conforming discretization Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 189 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0π/12 π/6 π/4 π/3 5π/12 π/2 (u * n)1-(u * a)1 α u1 NA<NB (B) NA>NB (A) NA<NB (A) NA>NB (B) NA=NB Figure 11: Ring under circumferential shear. Displacement errors, meshes 12–28. provides substantially smaller errors than the nonconforming discretizations. When comparing displacement and traction solutions obtained at the finer meshes of the nonconforming discretizations, Figures 10 and 12, it is seen that the oscillatory behaviour of displacements in the case NA>NBis substantially magnified in the case of tractions, whereas traction oscillations in the case NA<NBare hardly visible. The substantial difference between the magnitude of errors in these two cases is clearly seen in Figure 13 (bottom), where, due to the above mentioned symmetry, only the component t1is plotted. 0 1 2 3 4 5 0π/12 π/6 π/4 π/3 5π/12 π/2 0 1 2 3 4 5 t1 t2 α t1t2 t1 t2 NA<NB (B) NA>NB (A) analytical Figure 12: Ring under circumferential shear. Traction distributions, finer meshes, meshes 12– 28. -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0π/12 π/6 π/4 π/3 5π/12 π/2 (tn)1-(ta)1 α t1 NA<NB (B) NA<NB (A) NA>NB (B) NA=NB -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0π/12 π/6 π/4 π/3 5π/12 π/2 (tn)1-(ta)1 α t1 NA<NB (B) NA>NB (A) Figure 13: Ring under circumferential shear. Traction errors, meshes 12–28, all meshes (top), finer meshes (bottom). Traction error distributions at all meshes, except for the finer mesh in the case NA>NB(which has been omitted due to being out of the applied scale), shown in Figure 13 (top), have similar magnitude. All these errors seem to have some oscillations, the smallest errors being achieved in the case of the conforming mesh, though close to the end-points of Γc(corners of the adjacent subdomains) there are significant oscillations also in this case. The character of the strong oscillations of traction errors in the case NA>NBis illustrated considering the finest mesh 96–224, with the pattern 3–7 shown in Figure 14 (top) for the coarsest case. It should be noted, and it is also visible in the picture, that there are gaps and overlappings along the interface separating the discretized subdomains. The inner nodes of the pattern are marked by triangles, the end-points be- 190 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 ing marked by circles. Surprisingly smooth lines in Figure 14 (bottom) represent traction errors obtained at the subsets of the pattern nodes (denoted by numbers 1 to 6 or letters aor baccording to the top drawing) repeated along Γc, the substantially smaller errors at the coarser mesh being multiplied by an ad-hoc factor of value 100. For example, the two lines of circles representing the end points of the pattern for both the finer and coarser meshes coincide, the character of the error at those points being the same but the magnitude being substantially different. 12 3 4 5 6 a b -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0π/12 π/6 π/4 π/3 5π/12 π/2 (tn)1-(ta)1 α 1 2 34 5 6 a b t1 fine (A) coarse (B) Figure 14: Ring under circumferential shear. Mesh pattern ’3–7’ (top), traction errors in the case NA>NB, mesh 96–224 (bottom). 5.3.3 Square plate under uniaxial compression This example differs from the previous ones in that the curved interface is closed, a possible influence of the interface end-points on the error behaviour thus being eliminated. The displacements along Γcobtained at the finer meshes of the nonconforming discretizations are shown in Figure 15, an excellent approximation of the analytical solution again being obtained. Notice that expressions of u1and u2along Γcin terms of α include, due to the curvature of Γc, trigonometric functions. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0π/4 π/2 3π/4 π5π/4 3π/2 7π/4 2π 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 u * 1 u* 2 α u1 u2 u1 u2 NA<NB (B) NA>NB (A) analytical Figure 15: Square plate under uniaxial compression. Displacement distributions, finer meshes, meshes 24–56. Figure 16 shows the displacement error distributions, a strong oscillatory behaviour of the finer mesh in the case NA>NBfollowingthe discretization pattern, 3–7 in the present case, being observed as in the previous examples. -3.0.10-3 -2.5.10-3 -2.0.10-3 -1.5.10-3 -1.0.10-3 -5.0.10-4 0.0.100 5.0.10-4 1.0.10-3 0π/4 π/2 3π/4 π5π/4 3π/2 7π/4 2π-2.0.10-3 -1.5.10-3 -1.0.10-3 -5.0.10-4 0.0.100 5.0.10-4 1.0.10-3 1.5.10-3 2.0.10-3 (u * n)1-(u * a)1 (u* n)2-(u* a)2 α u1 u2 u1 u2 NA<NB (B) NA>NB (A) NA<NB (A) NA>NB (B) NA=NB Figure 16: Square plate under uniaxial compression. Displacement errors, meshes 24–56. Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 191 The behaviour of the displacement errors in Figure 16 for the rest of the meshes is rather smooth. A small consistency error appears in the case NA<NB, the compatibility condition being only slightly violated. Although the error achieved by the conforming-mesh is apparently the smoothest one, the error at the coarser mesh in the case NA>NBis surprisingly somewhat smaller. Tractions obtained at finer meshes of the nonconforming discretizations, Figure 17, show similar behaviour to the previous examples, oscillations being visible in the case NA>NBin t2andalsoin t1(whose analytic value is zero), though smaller ones, as could be expected. -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 0π/4 π/2 3π/4 π5π/4 3π/2 7π/4 2π -1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 t1 t2 α t1 t2 t1 t2 NA<NB (B) NA>NB (A) analytical Figure 17: Square plate under uniaxial compression. Traction distributions, finer meshes, meshes 24–56. Figure 18 (bottom) shows traction errors at finer meshes of the nonconforming discretizations, typically strong oscillations being observed in the case NA>NB. Oscillationsof traction errors in the case NA<NB, almost invisible in Figure 18 (bottom), can be seen in Figure 18 (top), where the results for the finer mesh of the case NA>NBhave been omitted, being out of the present scale. The repeating pattern of the errors is now also observable for the case NA<NB. The absence of endpoints on a smooth Γcyields a very smooth traction error distribution in the conforming case. Finally, it should be mentioned that, whereas in the previous two examples the choice of the subdomains ΩAand ΩBwas arbitrary, with no relevance to the numerical results, in the present -1.4.10-2 -1.2.10-2 -1.0.10-2 -8.0.10-3 -6.0.10-3 -4.0.10-3 -2.0.10-3 0.0.100 2.0.10-3 0π/4 π/2 3π/4 π5π/4 3π/2 7π/4 2π-8.10-3 -6.10-3 -4.10-3 -2.10-3 0.100 2.10-3 4.10-3 6.10-3 8.10-3 ( tn ) 1- ( ta ) 1 (tn)2-(ta)2 α t1 t2 t1 t2 NA<NB (B) NA<NB (A) NA>NB (B) NA=NB -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 0π/4 π/2 3π/4 π5π/4 3π/2 7π/4 2π -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 (tn)1-(ta)1 (tn)2-(ta)2 α t1 t2 NA<NB (B) NA>NB (A) Figure 18: Square plate under uniaxial compression. Traction errors, meshes 24–56, all meshes (top), finer meshes (bottom). example subdomain ΩAshould be the circle as shown in Figure 4(c). This is related to the problem of the invertibility of the SGBEM system for a single domain BVP with cavities subjected to traction (Neumann) boundary conditions, studied in Vodiˇcka et al. (2006). To avoid the above restriction on the choice of the interior subdomain, some of the techniques developed in Vodiˇcka et al. (2006) could be used. 5.4 Error convergence Convergence of displacement and traction errors along Γcin the discretized L2-andMAX-norms for h-refinements of boundary element meshes is studied for all the above examples. Dependence of an error norm on the number N, defined as the minimum number of elements along Γcbelonging to one of the bodies, i.e. N=min(NA,NB), is plotted in log-log scale. For each refinement, 192 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 given by N, results for the following three discretizations are presented: N=NA<NB,NA> NB=Nand N=NA=NB. Each line in the plot corresponds to a discretization of Γcand either to the coarser or to the finer mesh, a letter in parenthesis indicating the corresponding subdomain. The convergence rates β obtained by a linear regression using the three finest meshes are given in the key of the plots for each particular line. 5.4.1 Rectangular plate under exponential load The displacement error norms computed are presented in Figure 19, the magnitude of both error norms in all cases being quite similar. All convergence rates β ≃2, the expected quadratic convergence of displacements for linear boundary elements thus being verified in both norms. 10-6 10-5 10-4 10-3 10-2 10-1 10+0 10+1 10+2 10+3 ||un-ua||L2/ ||ua||L2 N=min(NA,NB) NA<NB (B) β=1.98556 NA>NB (A) β=1.98993 NA<NB (A) β=2.00343 NA>NB (B) β=1.98617 NA=NB β=1.98615 10-6 10-5 10-4 10-3 10-2 10-1 10+0 10+1 10+2 10+3 ||un-ua||MAX/ ||ua||MAX N=min(NA,NB) NA<NB (B) β=2.00764 NA>NB (A) β=1.98644 NA<NB (A) β=2.01439 NA>NB (B) β=2.00610 NA=NB β=2.00582 Figure 19: Rectangular plate under exponential load. Convergence rates for displacements in the L2-norm (top) and MAX-norm (bottom). Figure 20 presents the traction error norms obtained. Without looking at β values, it is clearly seen that the worst convergence behaviour corresponds to the finer mesh solution for NA>NBand the best one corresponds to the coarser mesh solution for NA<NB, the other three cases having similar convergence behaviour to each other. When comparing rates of convergence in the L2-norm, the worst case has only the linear convergence ( β ≃1), whereas the best case is not far from the quadratic convergence ( β ≃1.85), in the rest of the cases β ≃1.5. Convergence in the MAX-norm is in all cases linear, excepting the coarser mesh solution for NA<NBwhere β ≃1.37. It should be noticed that in the latter case the convergence becomes progressively worse with increasing refinement. The above differences in the convergence behaviour are related to very large oscillations of the error at the finer mesh in the case NA>NBand also to an increase in oscillations close to the end-points of Γc(corners of the adjacent subdomains), which happen in all cases, see Figure 9. This increase in oscillations naturally has stronger influence in the MAX-norm of error than in the L2-norm. In particular, the conforming mesh solution has relatively large errors close to these corners, causing the rate of convergence in the MAX-norm to be the lowest one. 5.4.2 Ring under circumferential shear Figure 21 shows the convergence of the displacement error norms. Convergence rates β ≃2for all cases, the expected quadratic convergence of displacements again being verified in both norms. The differences between the L2-norms of the error obtained in different cases are still relatively small, although somewhat greater than in the previous example. These differences are also small in the MAX-norm except for the conforming mesh solution, whose MAX-norm of the error is substantially smaller than in the nonconforming discretizations. The convergence of the traction error norms is plotted in Figure 22, the worst results again being obtained by the finer mesh for NA>NB. Only linear convergence is obtained at this mesh in the L2- norm, whereas for the rest of the meshes β ≃1.5. These lower rates of convergence in comparison Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 193 10-5 10-4 10-3 10-2 10-1 10+0 10+0 10+1 10+2 10+3 ||tn-ta||L2/ ||ta||L2 N=min(NA,NB) NA<NB (B) β=1.55246 NA>NB (A) β=1.01243 NA<NB (A) β=1.84407 NA>NB (B) β=1.50885 NA=NB β=1.55354 10-5 10-4 10-3 10-2 10-1 10+0 10+0 10+1 10+2 10+3 ||tn-ta||MAX/ ||ta||MAX N=min(NA,NB) NA<NB (B) β=1.06685 NA>NB (A) β=0.99857 NA<NB (A) β=1.36526 NA>NB (B) β=1.08798 NA=NB β=0.96263 Figure 20: Rectangular plate under exponential load. Convergence rates for tractions in the L2- norm (top) and MAX-norm (bottom). with displacements might be, at least partially, related to an increase of oscillation in tractions at all meshes close to the end-points of Γc(corners of the adjacent subdomains). Actually, the linear rate of convergence is obtained in the MAX- norm of the error at all meshes, the greatest error norm again being obtained at the finer mesh for NA>NB. 5.4.3 Square plate under uniaxial compression This example allows us to study whether the absence of end-points at Γc(corners of the adjacent subdomains) has an influence on the solution convergence. Convergence of the displacement error norms is plotted in Figure 23. The quadratic convergence of displacements in the L2-norm is again verified, as all β ≃2. Notice that the smallest β is achieved by the conforming mesh. Nevertheless, evaluating β in the MAX-norm, the conform- 10-6 10-5 10-4 10-3 10-2 10 -1 10+0 10+1 10+ 2 ||un-ua||L2/ ||ua||L2 N=min(NA,NB) NA<NB (B) β=1.99059 NA>NB (A) β=1.97424 NA<NB (A) β=2.01274 NA>NB (B) β=1.98175 NA=NBβ=2.01480 10-6 10-5 10-4 10-3 10-2 10 -1 10+0 10+1 10+ 2 ||un-ua||MAX/ ||ua||MAX N=min(NA,NB) NA<NB (B) β=1.98928 NA>NB (A) β=1.98028 NA<NB (A) β=1.99383 NA>NB (B) β=1.98790 NA=NBβ=1.85998 Figure 21: Ring under circumferential shear. Convergence rates for displacements in the L2- norm (top) and MAX-norm (bottom). ing mesh is the best one, with a very high rate of convergence β ≃3, the rest of the meshes giving β ≃2. It is also interesting that the magnitude of this error norm for the coarser mesh in the case NA>NBis comparable with the conforming case, though its convergence rate is the smallest. The convergence of the traction error norms is plotted in Figure 24. The worst and best convergences, approximately linear and quadratic, in both error norms have been obtained respectively by the finer mesh solution in the case NA>NB and the conforming mesh. As can be observed in Figure 18, the error for the coarser mesh in the case NA>NBslightly oscillates around the conforming mesh error, and consequently it could be expected that the error norms will have similar values at the both meshes. Actually, this is observed in the L2-norm, where β ≃1.9forthe coarser mesh in the case NA>NB; whereas in the 194 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 10-4 10-3 10-2 10-1 10+0 10+0 10+1 10+2 ||tn-ta||L2/ ||ta||L2 N=min(NA,NB) NA<NB (B) β=1.43993 NA>NB (A) β=1.01406 NA<NB (A) β=1.47822 NA>NB (B) β=1.49997 NA=NB β=1.47076 10-4 10-3 10-2 10-1 10+0 10+0 10+1 10+2 ||tn-ta||MAX/ ||ta||MAX N=min(NA,NB) NA<NB (B) β=0.90373 NA>NB (A) β=0.99813 NA<NB (A) β=1.04390 NA>NB (B) β=1.01572 NA=NB β=0.98021 Figure 22: Ring under circumferential shear. Convergence rates for tractions in the L2-norm (top) and MAX-norm (bottom). MAX-norm a decrease in the rate of convergence is observed in this case, β ≃1.5, due to the above mentioned oscillations. Both meshes in the case NA<NBpresent similar convergence behaviour, only slightly worse than in the conforming case, β ≃1.6and β ≃1.8intheL2-norm and MAX- norm respectively. 6 Conclusions The Symmetric Galerkin BEM approach applied to a Domain Decomposition BVP has been enriched by a possibility to formulate it through a variational principle, based on a search for a critical point of an energy functional which results in a weak formulation of the coupling conditions between the subdomains. An advantage of the method developed is that no additional unknowns, different from the natural variables of the BEM (displacements and tractions defined on the 10-6 10-5 10-4 10-3 10-2 10-1 10+1 10+2 ||un-ua||L2/ ||ua||L2 N=min(NA,NB) NA<NB (B) β=1.97667 NA>NB (A) β=1.99520 NA<NB (A) β=1.99533 NA>NB (B) β=1.99689 NA=NB β=1.95360 10-6 10-5 10-4 10-3 10-2 10-1 10+1 10+2 ||un-ua||MAX/ ||ua||MAX N=min(NA,NB) NA<NB (B) β=1.97694 NA>NB (A) β=1.98157 NA<NB (A) β=2.03508 NA>NB (B) β=1.91372 NA=NB β=3.04765 Figure 23: Square plate under uniaxial compression. Convergence rates for displacements in the L2-norm (top) and MAX-norm (bottom). subdomain boundaries), are defined at interfaces. This method has been tested numerically using continuous linear boundary elements, very good results having been obtained for both matching and non-matching meshes at straight and curved interfaces. To the authors’ knowledge, this study is the first to present a convergence study of numerical results obtained by SGBEM applied to elastic DDBVPs with non-matching meshes along curved interfaces. The convergence of the method for h-refinements has been thoroughly studied by three examples using the discretized L2and maximum norms of the error. The quadratic rate of convergence, β ≃2, of displacements has been obtained in both norms in all the cases, as could be expected for linear elements. Nevertheless, the convergence rates for tractions have in general been found a bit lower (somewhere between linear and Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 195 10-4 10-3 10-2 10-1 10+0 10+1 10+2 ||tn-ta||L2/ ||ta||L2 N=min(NA,NB) NA<NB (B) β=1.56589 NA>NB (A) β=0.97476 NA<NB (A) β=1.59949 NA>NB (B) β=1.93261 NA=NB β=2.00296 10-4 10-3 10-2 10-1 10+0 10+1 10+2 ||tn-ta||MAX/ ||ta||MAX N=min(NA,NB) NA<NB (B) β=1.83562 NA>NB (A) β=0.99769 NA<NB (A) β=1.78866 NA>NB (B) β=1.54240 NA=NB β=2.02885 Figure 24: Square plate under uniaxial compression. Convergence rates for tractions in the L2- norm (top) and MAX-norm (bottom). quadratic, 1  β 2), and dependent, in the case of non-matching meshes along the interface, on the choice of the subdomains, ΩAand ΩB,which control the tractions and displacements along the interface respectively. Reliable results in tractions on both sides of the interface, with β ≃1.5inthe L2-norm, have been obtained when the more refined mesh is situated at the interface side controlling the displacements, i.e. NA<NB. An explanation for this observation, using a reduction of the resulting linear system to the Schur complement system, has been introduced in Appendix A. In presence of the end-points at the interface (corners of the adjacent subdomains) the somewhat worse convergence of traction error norms observed, including the case of conforming meshes, can be, at least partially, attributed to stronger oscillations of computed tractions in the vicinity of these corners. Note that this phenomenon is usually present in BEM solutions, and thus is not specific for the present coupling approach. The solution of a problem with a smooth closed interface (without corners of the adjacent subdomains) has confirmed this explanation, as an improvement of the convergence has been observed, this improvement being substantialconsidering the maximumnorm of traction errors. Acknowledgement: The work has been performed during the stay of R.V. at the University of Seville, which has been supported by a NATO Scientific Committee Fellowship (Grant number OTAF03-054). The other financial support has been given to V.M. and F.P. from the Spanish Ministry of Science and Technology (Project MAT2003-03315) and to R.V. from the Scientific Grant Agency of the Slovak Republic (Grant number 1/1089/04). References Aimi, A.; Diligenti, M.; Lunardini, F.; Salvadori, A. (2003): A new application of the panel clustering method for 3D SGBEM. CMES: Computer Modeling in Engineering and Sciences,vol. 4, pp. 31–50. Baláš, J.; Sládek, J.; Sládek, V. (1989): Stress Analysis by Boundary Element Method. Elsevier, Amsterdam. Blázquez, A.; París, F.; Mantiˇ c, V. (1998): BEM solution of two dimensional contact problems by weak application of contact conditions with nonconforming discretizations. Int. J. of Solids and Structures, vol. 35, pp. 3259–3278. Blázquez, A.; Vodiˇ cka, R.; París, F.; Mantiˇ c, V. (2002): Comparing the conventional displacement BIE and the BIE formulationsof the first and the second kind in frictionless contact problems. Eng. Anal. Bound. Elem., vol. 26, pp. 815–826. Bonnet, M. (1995): Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity. Eng. Anal. Bound. Elem., vol. 15, pp. 93–102. 196 Copyright c  2007 Tech Science Press CMES, vol.17, no.3, pp.173-203, 2007 Bonnet, M.; Maier, G.; Polizzotto, C. (1998): Symmetric Galerkin boundary element method. Appl. Mech. Rev., vol. 15, pp. 669–704. Duddeck, F.M.E. (2006): An alternative approach to boundary element methods via the Fourier transform. CMES: Computer Modeling in Engineering and Sciences, vol. 16, pp. 1–14. Ganguly, S.; Layton, J.; Balakrishna, C. (2004): A coupling of multi-zone curved Galerkin BEM with finite elements for independently modelled sub-domains with non-matching nodes in elasticity. Int. J. Num. Meth. Engrg.,vol. 59, pp. 1021 – 1038. Ganguly, S.; Layton, J.; Balakrishna, C.; Kane, J. (1999): A fully symmetric multizone Galerkin boundary element method. Int. J. Num. Meth. Engrg., vol. 44, pp. 991 – 1009. González, J.A.; Abascal, R. (2000): Solving rolling contact problems using boundary element method and mathematical programming algorithms. CMES: Computer Modeling in Engineering and Sciences, vol. 1, pp. 141–150. Gray, L.; Paulino, G. (1997): Symmetric Galerkin boundary integral formulation for interface and multi-zone problems. Int. J. Num. Meth. Engrg., vol. 40, pp. 3085 – 3101. Guiggiani, M. (1995): Hypersingular boundary integral equations have an additional free term. Comput. Mech., vol. 16, pp. 253 – 260. Guiggiani, M.; Krishnasamy, G.; Rudolphi, T.; Rizzo, F. (1992): A general algorithm for the numerical solution of hypersingular boundary integral equations. J. Appl. Mech., vol. 59, pp. 604 – 614. Han, Z. D.; Atluri, S. N. (2002): SGBEM (for cracked local subdomain) – FEM (for uncracked global structure) alternating method for analyzing 3D surface cracks and their fatigue-growth. CMES: Computer Modeling in Engineering and Sciences, vol. 3, pp. 699–716. Hartmann, F. (1985): The Mathematical Foundation of Structural Mechanics. Springer–Verlag, Berlin. Hsiao, G.; Steinbach, O.; Wendland, W. (2000a): Domain decomposition methods via boundary integral equations. J. Comp. Appl. Math., vol. 125, pp. 521–537. Hsiao, G.; Schnack, E.; Wendland, W. (2000b): Hybrid coupled finite-boundary element methods for elliptic systems of second order. Comput. Methods Appl. Mech. Engrg., vol. 190, pp. 431– 485. Jaswon,M.A.;Symm,G.T.(1977): Integral Equation Methods in Potential Theory and Elastostatics. Academic Press, London. Jiao, X.; M.T.Heath (2004): Commonrefinement based data transfer between nonmatching meshes in multiphysics simulations. Int. J. Num. Meth. Engrg., vol. 61, pp. 2402–2427. Maier, G.; Diligenti, M.; Carini, A. (1991): A variational approach to boundary element elastodynamic analysis and extension to multidomains problems. Comput. Methods Appl. Mech. Engrg., vol. 92, pp. 193 – 213. Mantiˇ c, V.; París, F. (1995): Existence and evaluation of the free terms in the hypersingular boundary integral equation of potential theory. Eng. Anal. Bound. Elem., vol. 16, pp. 253–260. Mantiˇ c, V.; París, F. (1997a): Symmetrical representation of stresses in the Stroh formalism and its application to a dislocation and a dislocation dipole in an anisotropic medium. J. Elast.,vol. 47, pp. 101–120. Mantiˇ c, V.; París, F. (1997b): Symmetry properties of the kernels of the hypersingular integral and the corresponding regularized integral in the 2D Somigliana stress identity for isotropic materials. Eng. Anal. Bound. Elem., vol. 20, pp. 63– 168. McDevitt, T.; Laursen, T. (2000): A mortarfinite element formulation for frictional contact problems. Int. J. Num. Meth. Engrg., vol. 48, pp. 1525–1547. Núñez, J.; Muñoz, D.; Mantiˇ c, V.; París, F. (2003): A weak application of transmission conditions in 2D BEM with non-conforming meshes Symmetric Variational Formulation of BIE for Domain Decomposition Problems in Elasticity 197 for isotropic and anisotropic potential problems. In Gallego, R.; Aliabadi, M. (Eds): Advances in Boundary Element Techniques IV, pp. 111–116. Panzeca, T.; Salerno, M.; Terravechia, S. (2002): Domain decomposition in the symmetric boundary element analysis. Comput. Mech., vol. 28, pp. 191 – 201. París, F.; Cañas, J. (1997): Boundary Element Method, Fundamentals and Applications. Oxford University Press, Oxford. Puso, M. (2004): A 3D mortar method for solid mechanics. Int. J. Num. Meth. Engrg., vol. 59, pp. 315–336. Qian, Z.; Han, Z.; Atluri, S. (2004): Directly derived non-hyper-singular boundary integral equations for acoustic problems, and their solution through Petrov-Galerkin schemes. CMES : Computer Modeling in Engineering and Sciences, vol. 5, pp. 541–562. Quarteroni, A.; Valli, A. (1999): Domain decomposition methods for partial differential equations. Oxford University Press, Oxford. Schnack, E. (1987): A hybrid BEM model. Int. J. Num. Meth. Engrg., vol. 24, pp. 1015–1025. Sirtori, S. (1979): General stress analysis by means of integral equations and boundary elements. Meccanica, vol. 14, pp. 210–218. Springhetti, R.; Novati, G.; Margonari, M. (2004): Weak coupling of the symmetric Galerkin BEM with FEM for potential and elastostatic problems. CMES: Computer Modeling in Engineering and Sciences, vol. 13, pp. 67–80. Steinbach, O. (2003): Stability estimates for hybrid coupled domain decomposition methods. Springer-Verlag, Berlin Heidelberg. Vodi ˇ cka, R. (2000): The first-kind and the second-kind boundary integral equation systems for solutionof frictionless contact problems. Eng. Anal. Bound. Elem., vol. 24, pp. 407–426. Vodi ˇ cka, R.; Mantiˇ c, V. (2004): On invertibility of elastic single layer potential operator. J. Elast., vol. 74, pp. 147–173. Vodi ˇ cka, R.; Mantiˇ c, V.; París, F. (2006): On the removal of the non-uniqueness in the solution of elastostatic problems by symmetric Galerkin BEM. Int. J. Num. Meth. Engrg., vol. 66, pp. 1884 – 1912. Wendland, W.; Hsiao, G. (2004): Boundary element methods: Foundation and error analysis. In E. Stein, R. De Borst, T.J.R. Hughes(Ed): Encyclopedia of Computational Mechanics, Vol. 1 Fundamentals. John Wiley, New York, pp. 339- 373. Young, A. (1996): A single-domain boundary element method for 3-D elastostatic crack analysis using continuous elements. Int. J. Num. Meth. Engrg., vol. 39, pp. 1265–1293. Appendix A Analysis of the SGBEM formulation for a DDBVP with a single boundary The aim of this Appendix is to show how the discretized system (25) reduces to its Schur complement and to analyze the reduced system obtained. Although, for the sake of clarity, a simple DDBVP with Ω=R2partitioned to an interior (bounded) subdomain ΩAand an exterior (unbounded) subdomain ΩBseparated by a bounded interface curve, Figure A1(a), is considered, it is expected that the conclusions deduced from the following analysis can be extended to other more complex problems. ΩA ΩB Γc Ω+ Γ n+ t Ω− Γ n− n+ t (a) (b) (c) Figure A1: (a) A domain decomposition. Interior (b) and exterior (c) domains and pertinent normal vectors.