Symmetric Galerkin BEM for bodies with unconstrained contours

Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981 www.elsevier.com/locate/cma

Symmetric Galerkin BEM for bodies with unconstrained contours Francesco Freddi, Gianni Royer-Carfagni

*

University of Parma, Department of Civil–Environmental Engineering and Architecture, Parco Area delle Scienze 181/A, I-43100 Parma, Italy Received 31 July 2003; received in revised form 25 January 2005; accepted 14 February 2005

Abstract The solution via symmetric Galerkin boundary element method of problems in linear elasticity presents a major difficulty when a whole connected component of the boundary is under Neumann condition. In fact, apart from some special cases here discussed in detail, the resulting boundary integral equations are in general insensitive to the rigid body displacement of unconstrained connected boundary-components, even when the displacement is prescribed on some other boundary portions. A simple procedure is here proposed to solve the aforementioned indeterminacy while preserving the symmetry of the integral operator. Representative examples confirm the numerical effectiveness of the proposed method.
2005 Elsevier B.V. All rights reserved. Keywords: Boundary element method; Symmetric Galerkin; Periphractic bodies; Infinitesimal rigid body displacement; Uniqueness; Neumann boundary condition

1. Introduction For several years, the collocation method [1,2] has been the dominant numerical approach for solving boundary integral equations (BIEs). This technique, however, induces fully populated, non symmetric matrices in the corresponding linear system of equations. More recently [3,4], the symmetric Galerkin boundary element method (SG-BEM) approximation technique has been successfully applied to BIEs to furnish a symmetric integral formulation, ensuring uniqueness, stability and convergence of the numerical solution [5]. Compared with collocation methods, the SG-BEM is very expensive because of double *

Corresponding author. Tel.: +39 0521 905917; fax: +39 0521 905924. E-mail address: [email protected] (G. Royer-Carfagni).

0045-7825/$ - see front matter
2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.02.014

962

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

integrations, but the added symmetry makes this approach computationally efficient, both in the matrix construction and in the solution stage [6]. To our knowledge, a limited attention has been paid to the SG-BEM solution of problems in linear elasticity with unconstrained (Neumann) boundary. Since the elastic solution is defined up to an infinitesimal rigid body displacement (IRBD), the resulting pseudo-stiffness matrix is inevitably singular. A traditional countermeasure consists in restraining the boundary with statically-determined constraints, i.e., in such a way that the elastic deformation of the boundary does not induce any constraint reaction but the IRBD is restrained. Because of this property, such constraints are sometimes referred to as ‘‘dumb constraints’’. From a numerical point of view, the matrix singularity is removed by dropping in appropriate nodes of the discretized boundary a sufficient number of shape functions which, when set to zero, provide statically determinate constraints [4]. However, in some cases it may be difficult to determine which nodes to select to preserve the matrix symmetry. The situation is more complicated in the case of bodies which are periphractic,1 i.e., bodies whose boundary encloses one or more excluded holes. If the boundary conditions on a whole connected boundary-component (CBC), i.e., the external body surface or the surface of one of the internal excluded holes, are of the Neumann type, only comes the traction BIE into play for that contour in the SG-BEM formulation. But this equation remains satisfied by the addition of any IRBD to the displacement field of that CBC, even if some other part of the boundary is conveniently constrained. In other words, SG-BEM solutions can only be defined up to the IRBD of unconstrained CBCs, since the SG-BEM formulation itself is insensitive to the relative IRBD of disjoint contours. This inconvenience is not present in classical linear elasticity, where the compatibility with the material deformation uniquely determines the relative motion of the CBCs. Notice that with the collocation approach [2] this problem does not appear, because the corresponding equations, generated through the displacement BIE, are instead sensitive to the relative rigid displacement of disjoint boundaries. The aforementioned difficulty in the SG-BEM formulation is well-known, but it should be mentioned that there is one special case which, to our knowledge, has been only partially appreciated.2 Indeed, the IRBD of any Neumann CBC remains uniquely defined whenever it surrounds another contour which is, in turn, constrained. For example, consider a sphere with an inner hole. The IRBD of the internal CBC remains undetermined if this is subjected to conditions of the Neumann type. On the other hand, constraining a convenient portion of the inner cavity uniquely determines the displacement field on the external surface, whatever its boundary conditions are. Consequently, the IRBD of the external CBC under Neumann conditions remains undetermined only when the whole boundary (external and internal) is constraint-free. In the following, Neumann CBCs not surrounding constrained contours, whose IRBD remains undetermined in the SG-BEM formulation, will be referred to as ‘‘floating’’ contours or ‘‘floating’’ CBCs. To overcome this indeterminacy, in [8] the two-dimensional problem of a doubly-connected domain under Neumann conditions is solved by introducing an auxiliary problem, where each one of the CBCs (both internal and external) is subjected to additional dumb constraints. The traction BIE is used to solve the auxiliary problem and, successively, the displacement field for the original problem is obtained from weighted-averages using the displacement BIE. This procedure should be discussed on a purely theoretical ground. Because of the above-mentioned reasons, in order to prevent the floating of the outer boundary it would be sufficient to add dumb constrains to the inner boundary only. If both boundaries were

1 Milne-Thomson [7] appears to be the first to use the expression ÔperiphracticÕ in a topological sense. He noted its Greek derivative meaning Ôfenced aboutÕ and characterized a periphractic 3-dimensional region as one that is ‘‘bounded internally by one or more closed surfaces’’. Thus, a periphractic 3-dimensional region has embedded holes but not holes that pierce completely through it. The former is simply connected but the latter is not. A non-periphractic 3-dimensional region cannot have embedded holes which form internal boundaries and it need not be simply connected. 2 This point will be outlined in Section 3.

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

963

dumb-constrained, the body would be submitted, in general, to redundant constraints with non-null reactions and, consequently the original problem would be altered. However, the procedure presented in [8] is correct from a practical point of view. In fact, the dumb constrains referred to in [8] have a very special meaning, which reflects their practical numerical implementation. In particular, each constrained degree of freedom is not intended to fix the value of the corresponding displacement component, but it is implicitly made equivalent to dropping the corresponding equation in the linear system of equations obtained after discretization, so that the singularity of the pseudo-stiffness matrix is removed. With this procedure, the static role of the constraint is somehow softened and no redundant constraint reaction is generated even in the critical case. A more general approach to the same problem has been proposed by Aliabadi et al. [9], and consists in writing the displacement BIE on a limited portion of each unconstrained CBC. This procedure successfully solves the problem, but destroys the symmetry of the system matrix. In another paper by the same authors [10], an alternative method was proposed based on multi-zone formulation with interface conditions. This preserves the symmetry of the coefficient matrix, but implies the introduction of fictitious boundaries in order to render simply-connected a two-dimensional multi-connected domain. This technique increases the dimension of the problem and, what is more, the implementation of the matching conditions at the interface between contiguous domains considerably complicate the solution procedure. The aim of this paper is to propose an alternative approach to get through this problem, whose main advantage is that it removes the singularity of the resulting pseudo-stiffness matrix without introducing any artificial dumb constraint. Indeed, the adopted strategy relies upon the adaptation to the case of SG-BEM of a procedure proposed by Veˆrchery to regularize the equilibrium matrix of statically undetermined structures [11]. This approach is somehow similar to a computational technique used in finite elements based on ‘‘reduced-integration’’ [12] in order to stabilize the so called ‘‘kinematic modes’’ (a paradigmatic example of these is the hourglass mode in four-node bilinear elements using one-point quadrature [13]). In the proposed method for SG-BEM, automatically-generated additional terms are added to the original pseudo-stiffness matrix while preserving its symmetry. The resulting system matrix is not any more singular and its (unique) solution differs from the solution of the elastic problem only up to the IRBD of the floating contours. The true elastic solutions is determined a posteriori using a post-processing task based upon the displacement BIE. The advantages of the method are illustrated in representative examples.

2. Synopsis of symmetric formulation of the boundary integral equations Let X
Rd , with d = 2 or 3 the dimension of the space where the problem is set, represent the undistorted natural reference configuration of a homogeneous isotropic linear elastic body, delimited by the (possibly disconnected) boundary C with outward unit normal l. Prescribed displacement and applied traction fields on C are indicated by u and p, respectively. The boundary integral formulation of a linear elastic problem has been proposed by Rizzo [1]. The approach is based upon two fundamental integral equations stemming from SomiglianaÕs identity, which correlate interior displacement with boundary quantities. The following notation is the same of [14]. The first equation is usually referred to as the ‘‘displacement representation formula’’. For isotropic linear elastic solid with zero volume forces, the displacement u(x) at x 62 C is expressed in terms of the boundary fields u(y) and p(y), y 2 C, by Z Z vX ðxÞuðxÞ þ Gup ðx; y; lðyÞÞuðyÞ dCy ¼ Guu ðx; yÞpðyÞ dCy ; x 62 C. ð1Þ C

C

where vX(x) is the characteristic function of X, i.e. vX(x) = 1 if x 2 X, vX(x) = 0 otherwise. Here the kernels Guu(x,y) and Gup(x, y; l(y)) are related to KelvinÕs classical solution for an elastic homogeneous isotropic

964

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

infinite body. In particular, [Guu(x, y)]ji and [Gup(x, y; l(y))]ji represent the ith component of the displacement vector at y and of the traction vector across a surface element at y with unit normal l(y), due to a concentrated unit force acting at x and parallel to the j coordinate axis. For the following consideration, it may be useful to recall that the representation (1) can be obtained through BettiÕs reciprocal theorem [4]. In fact, the Kelvin state for an infinite body also corresponds to the elastic solution for the body X subjected to a unit force acting in the jth direction at x 2 X, as well as to tractions [Gup(x, y; l(y))]ji acting at point y of the boundary C with outward unit normal l(y). In particular, the corresponding displacement of y 2 C is given by [Guu(x, y)]ji. Let now consider a second elastic state of X, for which [u(y)]i and [p(y)]i represent the ith displacement and traction field components on y 2 C and [u(x)]i the displacement components at x 2 X. An application of BettiÕs theorem to the body X subjected to the two aforementioned elastic states, directly furnishes (1). The second fundamental identity commonly used in BEM method is the ‘‘traction representation formula’’, i.e. Z Z vX ðxÞpðx; nðxÞÞ þ Gpp ðx; y; nðxÞ; lðyÞÞuðyÞ dCy ¼ Gpu ðx; y; nðxÞÞpðyÞ dCy ; x 62 C. ð2Þ C

C

The kernels Gpu(x, y; n(x)) and Gpp(x, y; n(x); l(y)), usually referred to as GebbiaÕs kernels, define the elastic state in an infinite body due to a unit concentrated relative displacement at x across a surface element of normal n(x) (unit layered distortion). In particular, [Gpu(x, y; n(x))]ji represents the ith component of the displacement at point y 2 C when the aforementioned discontinuity at x is in the jth displacement component. Similarly, [Gpp(x, y; n(x); l(y))]ji is the ith component of the traction vector at point y 2 C (of normal l(y)) due to the same displacement discontinuity at x in jth direction. Eq. (2) can be obtained using again BettiÕs theorem, by considering the static state of the body X subjected to a unit layered distortion at x 2 X and to surface tractions on y 2 C consistent with the definition of Gpp(x, y; n(x); l(y)). The displacement and the traction BIEs are derived from (1) and (2) by performing the limit to boundary, i.e. x 2 X ! x 2 C. In this limit process, the singularities of the kernels are activated and integralsRhave to be understood in a generalized sense. In particular, with standard notation [15], denoting with - the Cauchy principal value and with the Hadamard finite part of the closed surface integral, Eqs. (1) and (2) become [16,17] Z Z CðxÞuðxÞ þ - Gup ðx; y; lðyÞÞuðyÞ dCy ¼ Guu ðx; yÞpðyÞ dCy ; x 2 C; ð3Þ C

C

Z Gpp ðx; y; nðxÞ; lðyÞÞuðyÞ dCy ¼ - Gpu ðx; y; nðxÞÞpðyÞ dCy ;

DðxÞpðxÞ þ C

x 2 C.

ð4Þ

C

This representation follows from the order of the kernelsÕ singularity: Gup and Gpu are strongly singular kernels and Gpp is the hypersingular kernel [15]. The coefficients C(x) and D(x) depend upon the degree of regularity of the boundary and, in particular, they are equal to 12 I for smooth boundaries [17]. In the most general boundary value problem in elasticity, external tractions pðxÞ ¼
pðxÞ are assigned on Cp, whereas the displacements uðxÞ ¼
uðxÞ are prescribed on the complementary part Cu = CnCp. The common numerical solution approach consists in solving the displacement BIE (3) via collocation technique [2], where the unknowns are the values of p on Cu and of u on Cp. This formulation, however, is non-symmetric. The key to obtaining a symmetric system of solving equations is to exploit the symmetry properties of the kernel functions, combined with the Galerkin approximation scheme [3]. For, the shape functions /(x), employed to approximate the boundary unknowns, are also used as weight functions for the evaluations ‘‘on average’’ of (3) and (4). In particular, Eq. (3) is ‘‘averaged’’ only on that part of the boundary Cu with prescribed displacements, whereas Eq. (4) is ‘‘averaged’’ on the complementary boundary portion Cp, with prescribed tractions. The resulting system of integral equations thus becomes

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

Z

Z

965

Z /ðxÞ - Gup ðx; y; lðyÞÞuðyÞ dCy dCx

/ðxÞCðxÞuðxÞ dCx þ Cu C Z Z ¼ /ðxÞ Guu ðx; yÞpðyÞ dCy dCx ;

Cu

Cu

Z

/ðxÞDðxÞpðxÞ dCx þ

Cp

¼

Z Cp

ð5Þ

C

Z

Gpp ðx; y; nðxÞ; lðyÞÞuðyÞ dCy dCx

/ðxÞ Cp

C

Z /ðxÞ - Gpu ðx; y; nðxÞÞpðyÞ dCy dCx .

ð6Þ

C

We recall, in passing, that the symmetric formulation admits a variational approach, since the resulting equations correspond to saddle-points conditions for a certain functional [14].

3. The intrinsic indeterminacy of Neumann contours In most cases, the system (5) and (6) is a well-posed system of integral equations, i.e., modulo smoothness hypotheses, a unique solution exists and it depends continuously about the boundary data. There is however a congenital indeterminacy in the particular case that tractions are assigned on the entirety of one CBC. This spurious behavior is due to the intrinsic structure of the kernels. To illustrate, consider an infinite body in Rd , d = 2 or 3, and let c represent a smooth closed contour embracing the subregion domain x
Rd . Let urig : Rd ! Rd represent an infinitesimal rigid body displacement, i.e. urig is of the form urig ðxÞ ¼ u0 þ v0  x with u0 and v0 constant vectors. Then the following relationship holds: Z vx ðxÞurig ðxÞ ¼  Gup ðx; y; lðyÞÞurig ðyÞ dcy ;

ð7Þ

ð8Þ

c

where vx(x) = 1 if x 2 x, vx(x) = 0 if x 62 x. This formula could be verified analytically, but here we discuss its mechanical significance stemming from the principle of virtual work. Consider KelvinÕs solution for a unitary force acting at x 62 c in direction ej. For any y 2 c with outer unit normal l(y), [Gup(x, y, l(y))]ji ei represents,3 as recalled in Section 2, the traction vector acting on an infinitesimal portion of c in a neighborhood of y. Consider the equilibrium of x under the actions of such forces. Of course, if x 62 x, [Gup(x, y,l(y))]ji ei are equipollent to a null force with null moment and, consequently, the integral on the right hand side of (8) is obviously zero. On the other hand, if x 2 x, the sub-body x is in equilibrium under the action of the aforementioned tractions and the unit force in direction ej acting at x. It is simple to verify that principle of virtual work again gives (8). Let now x ! c from the interior of x. Performing the ‘‘limit to the boundary’’, as already done to obtain (3), it can be demonstrated that Z CðxÞurig ðxÞ ¼  - Gup ðx; y; lx ðyÞÞurig ðyÞ dcy ; ð9Þ c

where we emphasize that lx(y) is the outward unit normal to x. As in (3), CðxÞ ¼ 12 I for smooth c. The same result is obtained by letting x ! c from the exterior of x [16]. 3

From now on, we will assume the summation convention on repeated indexes.

966

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

Consider now two disjoint contours ce and ci embracing the regions xe and xi respectively, and consider the case when xe contains ci. It can be shown that
Z CðxÞurig ðxÞ if x 2 ce ;  - Gup ðx; y; lðyÞÞurig ðyÞ dcy ¼ ð10Þ urig ðxÞ if x 2 ci ; ce
Z 0 if x 2 ce ;  - Gup ðx; y; lðyÞÞurig ðyÞ dcy ¼ CðxÞurig ðxÞ if x 2 ci . ci

ð11Þ

Notice that here we have denoted with l(y) the outward unit normal to boundary of the region xn  xenxi, comprised between ce and ci. In particular, if y 2 ce then l(y) coincides with lxe ðyÞ, i.e. the outer unit normal to xe defined as in (9), but if y 2 ci one has that lðyÞ ¼ lxi ðyÞ, with lxi ðyÞ unit normal to xi. Using the property that Gup(x, y, l(y)) = Gup(x, y, l(y)), formulas (10) and (11) are easily obtained from (8) and (9). Of course, such formulas can be easily extended to an arbitrary number of contours. Let us now consider the Gebbia kernel Gpp(x, y, n(x), l(y)) appearing in (6). In this case, modulo some smoothness hypotheses for the contour c, the following identity holds Z 0 ¼  Gpp ðx; y; nðxÞ; lðyÞÞurig ðyÞ dcy ; 8x 62 c. ð12Þ c

To give a physical interpretation to this formula, recall that [Gpp(x, y,n(x),l(y))]ji ei represents the traction vector at a surface element at y of normal l(y) due to unit layered distortion, i.e. a unit discontinuity for the jth displacement component at x across a surface element of normal n(x). In order to produce such a distortion, we imagine to cut the body along an infinitesimal surface portion through x and normal to n(x) and, on the resulting free surfaces, we apply tractions dual in energy to the unit distortion we desire to produce. For example, if n(x) is parallel to the x-axis, we apply surface tractions rxx to produce a discontinuity in the x-component of the displacement, surface tractions sxy to achieve a discontinuity in the y-component and surface tractions sxz to generate a z-component discontinuity. But such tractions, being opposite in sign on each lips, have null resultant and null moment resultant. Consequently, since x 62 c and the sub-body x must be in equilibrium, the tractions along c have as well zero resultant and zero moment resultant whatever the source point x is located, be it inside or outside x. Therefore, from the principle of virtual work, identity (12) follows for either x 2 x or x 62 x, provided x does not belong to the interface c. By the same argument used to obtain (4), the limit-to-the-boundary process x ! x 2 c gives [16] 0 ¼  Gpp ðx; y; nðxÞ; lðyÞÞurig ðyÞ dcy ;

8x 2 c.

ð13Þ

c

The indeterminacy which is associated to considering Neumann conditions on the integrity of one of the body contours follows from the identities (10), (11) and (13). For example, consider a hollowed domain X with disconnected boundary, delimited by the external contour Ce and one interior contour Ci. Let u : oXp ! Rd and p : oXu ! Rd be a solution of the system of integral equations (5) and (6), corresponding to the elastic solution of the boundary value problem where tractions
p and displacements
u are assigned on oXp and oXu, respectively. Suppose first that Ci  oXp and Ce  oXu 5 B. Since (5) must be evaluated for x 2 oXu  Ce, we have from (11)1 that this equation remains satisfied if an arbitrary IRBD of Ci is added to u*. The same conclusion holds for (6) because of (13). On the other hand, assume now that Ce  oXp. From (13), Eq. (6) still remains satisfied if an arbitrary IRBD is added to Ce. If oXu 5 B, Eq. (5) must be evaluated at x 2 oXu  Ci, but from (10)1, (5) is not satisfied if an IRBD is added to Ce. In other words, the symmetric BIE formulation is insensitive to IRBD of the exterior contour if and only if no boundary portion is constrained (so that Eq. (5) is identically satisfied). Whenever some part of the

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

967

boundary is conveniently constrained, be it either on the external or the in internal contour, the motion of the external contour results uniquely determined. The extension to the case of bodies delimited by a certain number of CBCs is perfectly similar. The previous considerations may be summarized by saying that in the SG-BEM formulation • The displacement field is defined up to the IRBD of those interior CBCs under Neumann conditions. • The displacement field of the exterior CBC is uniquely determined whenever some portion of the boundary, be it external or internal, is properly constrained. This conclusion is similar to the one reported by Pe´rez-Gavila´n and Aliabadi [9], with one exception. Indeed, the authors do not seem to realize that the displacement field on the exterior contour is uniquely determined in the particular case when tractions are assigned on the whole exterior contour, but a non-null portion of one of the interior contours is constrained. 4. A symmetric boundary formulation for undetermined Neumann problems The SG-BEM is based upon a proper use of Eqs. (3) and (4). In particular, displacement integral equation (3) is written for x 2 Cu and the traction integral equation (4) for x 2 Cp. In this way, the unknowns become the displacements u on Cp and tractions p on Cu. The resulting system of integral equations may be written in the following compact form K½v ¼ f.

ð14Þ

Here K denotes a symmetric linear integral operator, v is the vector of the unknowns (displacements u on Cp and tractions p on Cu), while f is a term which depends upon the assigned boundary data (displacements
p on Cp). In the numerical approach, problem (14) is discretized using the Galerkin u on Cu and tractions
approximation scheme as in (5) and (6), but for the following considerations it may be convenient to refer to the aforementioned continuum formulation. Recall that the symmetric boundary integral formulation is insensitive, as discussed in Section 3, to the IRBD of full-Neumann CBCs that do not surround another constrained CBC. Consequently, since the kernel of K may contain IRBD of unconstrained contours, in general the operator K is singular. In this case, the numerical problem derived from (14) is inevitably ill-conditioned. Here, we present a strategy to bypass the aforementioned indeterminacy that, while preserving the symmetry of the operator, can be automatically implemented, in the sense that the user does not need to introduce any artificial additional constraint. This approach is an adaptation to SG-BEM of an idea which, to our knowledge, has been first advanced by Veˆrchery [11] in order to regularize the system of equilibrium equations for discrete structures. Here, our approach will be presented in abstract form referring to problem (14). The corresponding numerical implementation presents no difficulty and may be done using a standard approximation scheme a` la Galerkin. 4.1. Bodies with connected boundary In order to present the method, consider first the simplest case of a body X with connected boundary C under full Neumann conditions. Thus, v in (14) coincides with the boundary displacement field, while f represents the applied tractions. Moreover, the kernel of K contains the IRBD of the (external) boundary C  Cp. A basis for Ker(K) is given by4 3(d  1) independent (infinitesimal) rigid body displacement of 4

Recall that d, d = 2 or 3, represents the dimension of the space where the problem is set.

968

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

C. For example, for d = 3, the basis formed by the three translations and the three rotations around the reference axes reads
v1rig ¼ i;
v2rig ¼ j;
v3rig ¼ k;
v4rig

¼ zj þ yk;


v5rig

¼ xk þ zi;

ð15Þ
v6rig

¼ yi þ xj.

We define the inner product of any two vector fields va ; vb : C ! Rd as Z hva ; vb i :¼ va vb dC.

ð16Þ

C

Of course, there are many equivalent ways to define a scalar product. Actually, the particular boundary displacement that will be uniquely determined using our proposed scheme depends upon such a choice. But there is no inconsistency in this from a physical point of view, since the displacement field of bodies with unconstrained connected boundary can only be determined up to a IRBD, and our proposed method simply selects one of the possible IRBDs. On the other hand, we will show in the next section that, when the boundary is not connected, an ad hoc post-processing is necessary in order to recover the compatibility with the elastic deformation of the surrounding material. With respect to the defined scalar product (16), an orthonormal basis can be constructed from (15) using Gram–Schmidt orthonormalization algorithm. Let then ^vjrig ; j ¼ 1 . . . 3ðd  1Þ, denote such an orthonormal basis, coinciding with a set of orthogonal eigenvectors of K associated to the null eigenvalues, i.e. K½^vjrig  ¼ 0. The system (14) admits a solution iff h^vjrig ; fi ¼ 0; j ¼ 1 . . . 3ðd  1Þ;

ð17Þ

which, from a physical point of view, is equivalent to the equilibrium of the external forces applied to the body. To verify (17), let us denote with V the admissible space of vectors v for problem (14). Setting V ¼ V0 KerðKÞ, decompose v in the form v = v0 + v1, with v0 2 V0 and v1 2 Ker(K). Scalarly multiplying both sides of (14) with ^vjrig and using the linearity of the operator, (17) follows. Our strategy consists in considering, instead of (14), the auxiliary ‘‘regularized’’ problem Kr ½v ¼ f; with Kr ¼ K þ aK0 ;

ð18Þ

where a is an arbitrary positive scalar and the operator K0 is defined as 3ðd1Þ X

ð^vjrig  ^vjrig Þ½v;

ð19Þ

ð^vjrig  ^vjrig Þ½v :¼ h^vjrig ; vi^vjrig .

ð20Þ

K0 ½v :¼

j¼1

where

Notice that the operator ð^vjrig  ^vjrig Þ is the orthogonal projector onto ^vjrig 2 KerðKÞ. The connection between the original problem (14) and the regularized problem (18) consists in the fact that the operators K and Kr, have the same eigenvectors, though associated to different eigenvalues on Ker(K), that is K½v ¼ Kr ½v ¼ kv; if ; v eigenvector 62 KerðKÞ;

ð21Þ

where k is the associated eigenvalue, whereas K½v ¼ 0 and Kr ½v ¼ av; if v 2 KerðKÞ.

ð22Þ

In other words, every v 2 Ker(K) is an eigenvector for Kr associated to the a eigenvalue. For any non-zero a the operator Kr is bijective, so that the system of integral equations Kr ½v ¼ f

ð23Þ

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

969

admits the unique solution vr ¼ K1 vjrig ; fi ¼ 0; j ¼ 1 . . . 3ðd  1Þ, then vr is orthor ½f. In particular, since h^ j gonal to Ker(K), i.e. h^vrig ; vr i ¼ 0; j ¼ 1 . . . 3ðd  1Þ. This can be shown using an argument similar to that used for (17). Of course, the unique solution vr of (23) is nothing but one of the possible solutions of (14). But from the observations above, vr is also that particular solution of (14) which simultaneously satisfies Z j  ^vjrig vr dC ¼ 0; j ¼ 1 . . . 3ðd  1Þ. ð24Þ h^vrig ; vr i ¼ C

In words, among all possible fields which verify (14), vr is that field whose IRBD is zero ‘‘on average’’, in the sense determined by (24). In general, vr depends upon the form of scalar product (16). In fact, by selecting another equivalent expression for such a product, the form of the operator K0 in (18)–(20) changes and, consequently, also the solution vr of (23) may accordingly vary. From a practical point of view, the evaluation of K0 from (19) and (20) presents no particular difficulty, but determining the orthonormal basis ^vjrig ; j ¼ 1 . . . 3ðd  1Þ may slightly complicate the numerical implementation. For, it is important to observe that the base ^vjrig of Ker(K), used to defined the operator Kr in (18), does not need to be orthonormal. In fact, the strategy of the method consists in defining Kr such that, for every non-null v 2 Ker(K), then Kr[v] 2 Ker(K) and Kr[v] 5 0. This property, together with the orthogonality of f to Ker(K), implies that any solution vr of Kr[v] = f is as well orthogonal to Ker(K). Actually, an operator Kr with the aforementioned properties may be defined in the same manner as in (18)–(20) provided that ^vjrig is a basis for the kernel of K, but not necessarily orthonormal. However, if ^vjrig is not orthonormal, a generic vector v 2 Ker(K) is not necessarily an eigenvector for Kr associated to the a eigenvalue, i.e. in general Kr[v] 5 av in (22). Nevertheless, the fundamental property that Ker(K) 3 Kr[v] 5 0 whenever Ker(K) 3 v 5 0 is maintained. As a result, the solution vr of the corresponding system (23) preserves a characterization similar to (24). In general, the easiest way to implement the method is to use the vectors
vjrig ; j ¼ 1 . . . 3ðd  1Þ of (15) for the definition of K0 according to (19) and (20). This procedure has been adopted in the BEM code used to solve all the examples reported in Section 6. 4.2. Bodies with disconnected boundary If the boundary C of the body X is disconnected, the insensitiveness to IRBD is shared by all the internal contours with Neumann boundary conditions. This may increase the degree of singularity of the operator K of (14) since, in general, there may be more than one floating contours. This indeterminacy can be solved using the same procedure of Section 4.1, i.e. defining the operator Kr as Kr ¼ K þ

m X

ai Ki0 ;

ð25Þ

i¼1

where m represents the number of floating contours, a > 0, and Ki0 is the operator for the ith contour defined as in (19). The system Krv = f has a unique solution in terms of displacement, where the IRBD of each contour satisfies a condition similar to (24). In the case of domains with more than one CBC there is however a further basic difficulty for this approach. The displacement of the ith floating contour, uniquely determined by Krv = f, is not necessarily in agreement with the elastic solution of the boundary value problem. In fact, such indeterminacy is intrinsic in the symmetric boundary formulation but is not present in the parent elastic problem, where compatibility with the material deformation controls the relative motion of the CBCs. It should be mentioned, however, that the aforementioned IRBD indeterminacy does not affect the state of stress in the body calculated within the BEM framework. In fact, since the stress field is obtained through

970

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

Eq. (2), it is clear from (12) that the IRBD of any contour does not alter the value of the traction vector p(x, n(x)). To get rid of the indeterminacy, we propose the following strategy articulated in two steps. (i) Define the operator Kr as in (25) and solve the corresponding system of integral equations Kr[v] = f. (ii) Use the displacement field evaluated at step (i) to determine the displacement field of the parent elastic solution a posteriori, using equations obtained from SomiglianaÕs identity. To illustrate the method, let vr represent the unique solution of the system Kr[v] = f obtained from (25). In general, for mixed boundary conditions, vr collects both displacements and traction fields. The displacement field will be indicated with ur . It is then convenient to distinguish two cases. Case 1. Suppose that a non-zero portion of the boundary is constrained so that, despite SG-BEM indeterminacy, in the parent elastic problem displacements are uniquely determined. Because of the considerations set forth in Section 3, in the SG-BEM formulation the IRBD of the external contour Ce is uniquely defined be it constrained at some portion or free, but displacements remain undefined on those internal contours, say CNi , where only tractions are prescribed. Let Xi
Rd represent the domain encircled by CNi , that is the ith ‘‘hole’’ of X. The field urig;i : CNi ! Rd , represents the IRBD to be added to ur determined at step (i), to obtain the displacement field u* corresponding to the solution of the parent elastic problem. In other words,5 urig;i ¼ u  ur . To find urig,i consider a point x 2 Xi inside the ith hole. Since x 62 X, the identity (1) gives Z Z Gup ðx; y; lðyÞÞu ðyÞ dCy ¼ Guu ðx; yÞpðyÞ dCy . ð26Þ C

C

On the other hand, application of (8) to the IRBD urig;i ¼ u  ur , denoting with l(y) the outer unit normal to the domain X at y 2 C, gives6 Z urig;i ðxÞ ¼ Gup ðx; y; lðyÞÞðu ðyÞ  ur ðyÞÞ dCy CN i

Z

Gup ðx; y; lðyÞÞðu ðyÞ  ur ðyÞÞ dCy Z Z ¼ Guu ðx; yÞpðyÞ dCy  Gup ðx; y; lðyÞÞur ðyÞ dCy . ¼

C

C

ð27Þ

C

The significance of the representation (27) where urig,i(x) is evaluated inside the ‘‘hole’’ at a point x 62 X, is that it represents the IRBD that such point would have attained if the hole was made rigid and its boundary CNi was subjected to the displacement field urig;i : CNi ! Rd . Notice that in deriving the second line of (27) we have used that u  ur on the external boundary Ce and that, because of (8), the integrals of Gup ðu  ur Þ on those internal contours not encircling x are identically zero. The third line of (27) is a simple consequence of (26). In order to uniquely determine the displacement field, it is necessary to determine for each one of the ‘‘floating’’ contours 3(d  1) parameters, defining the IRBD of the type (7). For, it is necessary to write 3(d  1) independent equations that may obtained by evaluating (27) at some particular points of the domain holes.

5

Without risk of ambiguity, we will not change notation to indicate the restriction of u* and ur : C ! Rd to the contour CNi . Notice from the above definition of l(y) that at y 2 CNi the normal l(y) points inside the hole Xi , so that the sign on the right-hand side of (8) has to be changed. 6

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

971

As a matter of fact, it would be possible to evaluate urig,i directly at CNi ! Rd by performing the limit-tothe-boundary of (27). Using (3) and (11) and reasoning as in (26) and (27), one reaches an expression of the type Z Z 2CðxÞ urig;i ðxÞ ¼ Guu ðx; yÞpðyÞ dCy  Gup ðx; y; lðyÞÞur ðyÞ dCy CnCN

C

i Z   - Gup ðx; y; lðyÞÞur ðyÞ dCy  CðxÞur ðxÞ;

CN i

x 2 CNi .

ð28Þ

However, this approach involves the evaluation of strongly singular integrals and, from a numerical point of view, it is less attractive than the formulation (27), where all integrals are non-singular. Case 2. Suppose now that the whole boundary (both internal and external) is unconstrained, so that also the IRBD of the external boundary Ce remains undetermined. In this case, also the displacement field u*, solution of the parent elastic problem, can only be determined up to an infinitesimal rigid displacement of the whole body. But once the displacement of one of the contours is assigned, the IRBD of the remaining contours can be uniquely obtained using the procedure sketched for case 1. In particular, one could consider the solution vr  ur of the system7 Kr[v] = f, retain as assigned such a displacement on one of the contours, either internal or external, and afterwards follow the same procedure of case 1. In substance, one has to apply regularizing term ai Ki0 in (25) to each one of the m body contours, and apply the displacement Eq. (27) to m  1 contours. It is somewhat arbitrary to select the contour for which the displacement field obtained from the solution of the system Kr[v] = f should be regarded as true. In particular, apart from computational reasons, there is no theoretical advantage in selecting an internal rather than an external contour. Indeed it should be noticed that applying the regularizing term to one contour is different from constraining that contour. In fact, in this second case the displacement BIE (3) comes into play, whereas in the first case we only use the traction BIE (4). Let us suppose that the displacement of one contour Ck, either external or internal, is retained as a datum. The solution u* obtained through the aforementioned procedure verifies Z ~vjrig u dCk ¼ 0; j ¼ 1 . . . 3ðd  1Þ. ð29Þ h~vjrig ; u i :¼ Ck

Here, ~vjrig ; j ¼ 1 . . . 3ðd  1Þ represent a basis (not necessarily orthonormal) of the type (15) for the IRBD of Ck, which has been introduced to define K0 as per (18)–(20). 4.3. Geometric interpretation Geometric visualization may perhaps be of help to clarify the procedure outlined in Sections 4.1 and 4.2. Let V represent the set of the admissible displacement fields for the boundary C of the body X. Consider first the case when C is connected as in Section 4.1. Let then D? represent the subset of V generated by all the IRBD of C, which will be conventionally represented by a vertical vector in Fig. 1. Now, since the SG-BEM formulation is insensitive to any arbitrary IRBD of the boundary, the vector K[u] is orthogonal to D? for any u 2 V. Consequently, the system (17) admits a solution if and only if f is as well orthogonal to D?, as represented in Fig. 1. On the other hand, the vector Kr[u], defined as in (23), has in general a non-zero component in D? and, more in particular, it is orthogonal to D? if and only if u is orthogonal to D?. It follows that field vr , solution of the regularized problem of Kr[v] = f of (4.1) must 7

Now, the unknowns are the displacements on the whole boundary, so that vr  ur .

972

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

D⊥ Kr [u]

u

f

K[u]

Fig. 1. Geometrical interpretation of the regularization procedure for a body with connected boundary.

B

Kr [u]

D0⊥

De⊥

So

mi

ia na

Di⊥

f

So

mi

De⊥

gl

u

D0⊥

⊥

lo cu s

Di

O

O

gli

an

al

A

oc

us

C

K[u]

a

Π⊥

Π⊥

b

Fig. 2. (a) Geometrical interpretation of the regularization procedure for a body with two disjoint contours; (b) orthogonal view of the plane P?.

be orthogonal to D? as well. In substance, the regularization procedure described in Section 4.1 selects, among all possible solutions of the system (14), the unique field that is orthogonal to D?. If C is not connected as in Section 4.2, a basis for Ker(K) comprehend all the independent IRBD of the floating contours. To illustrate, consider the simplest case when X is delimited only by two disjoint contours, one internal and one external, and supposed as in case 2 of Section 4.2 that both contours are subjected to Neumann boundary conditions. ? ? ? For such a body KerðKÞ  D? e Di , where De and Di represent the subspace generated by the IRBD of the exterior and interior contours respectively. Observe in passing that if D? 0 denotes the subspace generated ? ? by the IRBD of the body X considered as a whole, then D? $ D

D since the presence of the material 0 e i restrains the relative motion of the contours. In Fig. 2 the subspace D? is represented by a vertical vector, 0 ? ? ? ? while KerðKÞ  D?

D is indicated by the plane P passing through D . Subspaces D? e i 0 e and Di are as ? well indicated as vectors, defining a basis for P . Since the body is in equilibrium, one would expect that the virtual work of the applied forces for any IRBD of the whole body should be zero. Translated in the SG-BEM, the field f should be orthogonal to ? ? D? 0 . Remarkably, it can be demonstrated that f is orthogonal to both De and Di , even when the resultant and moment resultant of the applied forces on each one of the body contours are not zero.8 This is why in Fig. 2 the vector f as been drawn at right angle to P?. To show this, notice that since C  Cp, only the Eq. (4) enters the system K[v] = f of (14), which becomes Z Gpp ðx; y; nðxÞ; lðyÞÞvðyÞ dCy ¼ DðxÞpðxÞ þ - Gpu ðx; y; nðxÞÞpðyÞ dCy . ð30Þ C

C

8 In absence of body forces, equilibrium requires that the resultant and moment resultant of the forces applied to both contours is zero.

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

973

The term of the right-hand side of this equations represents the vector f. Let now vr represent the solution of the regularized problem Kr[v] = f of (23). On the other hand, we will denote with u* one of the solution of the parent elastic problem (where the relative motion of the contours is constrained by the material presence). Clearly, because of (13) and the considerations set forth in Section 3, both v = u* and v ¼ vr solves the system (30). Thus, we can write Gpp ðx; y; nðxÞ; lðyÞÞur ðyÞ dCy .

fðxÞ 

ð31Þ

C

Now, let ubrig : Cb ! Rd represent a generic IRBD of either the internal or the external boundary, here generically denoted with Cb. Using the symmetry of the kernel Gpp, i.e. Gpp ðx; y; nðxÞ; lðyÞÞ ¼ GTpp ðy; x; lðxÞ; nðyÞÞ;

ð32Þ

with a proper limit-to-the-boundary procedure we obtain from (31) and (32)   Z Z ubrig ðxÞ fðxÞ dCx ¼ ubrig ðxÞ Gpp ðx; y; nðxÞ; lðyÞÞur ðyÞ dCy dCx Cb

Cb

¼

Z C

"

C

Gpp ðy; x; lðyÞ; nðxÞÞubrig ðxÞ dCx

ur ðyÞ

# dCy .

ð33Þ

Cb

But the last term is evidently null because of (13). Consequently, from (16), hubrig ; fi ¼ 0, which proves ? that f is orthogonal to both D? e and Di , as drawn in Fig. 2a. The interpretation of the regularization procedure is then similar to that of Fig. 1. Since K½ubrig  ¼ 0 for any ubrig : Cb ! Rd , we have that the imagine of a generic vector u through the operator K is orthogonal to P?, i.e. K[u] is parallel to f. This is another indication that, for the existence of solutions, f must be orthogonal to P?. On the other hand, Kr[u] has a non-zero component on P? if and only if u has as well a nonzero component on the same plane. Consequently, following the procedure of Section 4.2, the unique solution vr of Kr[v] = f is orthogonal to P? as well. But the elastic solution has to comply with condition (27). Such a condition, obtained from SomiglianaÕs ? identity, may be geometrically visualized as an affine correlation between D? e and Di and, consequently, ? may be interpreted as a line on the plane P . In Fig. 2a and b such a line is referred to as ‘‘Somigliana locus’’. The solution vr is also a solution of the parent elastic problem only if the component of vr on P? belongs to the Somigliana locus. Since the component of vr on P? is null, this can only happen if the Somigliana locus passes through the origin ‘‘O’’, but this eventuality is very special. In general, the Somigliana locus does not passes through the origin so that vr must be modified. In particular, we add to vr a vector lying P?, which has to be calibrated so that the tip of its arrow touches the Somigliana locus. In general, there are many ways to do this. For example, if we imagine to fix the displacement field of the external contour, this is equivalent to say that the component of the IRBD on D? e is null. Consequently, the IRBD of the internal contour to added to vr to achieve a solution of the parent elastic problem is represented by point A in Fig. 2b. If, on the other hand, it is the interior contour to be considered fixed, the point representative of the IRBD to be added to the exterior contour coincides with B in Fig. 2b. In general, we may decide to impose an IRBD to the external and to the internal contour. Once the IRBD has been assigned on one of the two contours, the IRBD to be assigned to the other contours follows from the compatibility with the Somigliana locus. For example, if we prescribe that the IRBD of the whole body X has to be zero, the displacement component on P? must be orthogonal to D? 0 . The corresponding point on the Somigliana locus is represented by C in Fig. 2b. The consequent IRBD to be ! ? assigned to the internal and external contours is represented by the components of OC on D? i and De , respectively.

974

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

5. The SG-BEM numerical approach The discretization of the SG-BEM system of integral equations (5) and (6) leads to a system of linear equations of the form " #  P KU KUP K ð WÞ ¼ ¼ F. ð34Þ PU P U K K T

Here, KU ; KP are symmetric matrices, KUP ¼ ðKPU Þ and, consequently, K is symmetric, F is representative of the assigned data, while W is the vector of the nodal unknowns. In particular, W  (P; U)T, where U and P represent the Nu displacement-nodal-values on Cu and Np traction-nodal-values on Cp respectively. The indeterminacy associated with closed contours with Neumann conditions implies that, in such cases, K is singular and the system (34) is numerically ill-conditioned. As recalled in the Introduction, a common technique to remove the singular character of the matrix consists in constraining some boundary portions with statically-determinate constraints (‘‘restraint approach’’). From a numerical point of view, this technique has two main disadvantages: it destroys the symmetry of the system matrix and, in addition, it requires the somehow arbitrary selection of the degrees of freedom to be eliminated. The following alternative technique, based upon the results of Section 4, has some potential advantages, i.e. it can be applied routinely and it is numerically more attractive than the restrain approach, especially for problems without obvious symmetries. In general, consider a body with disconnected boundary and let M represent the number of ‘‘floating’’ contours. Let ujrig;b ; j ¼ 1 . . . 3ðd  1Þ, represent a basis for the IRBD of the bth of such contours, referred to as Cb. Consider then the regularized system Kr[v] = f, defined as per (18)–(20), (23) and (25). Noticing that ujrig;i ðyÞ  0 if y 62 Cp , the part of this system resulting from the traction BIE (4) becomes Z 3X ðd1Þ M X Gpp ðx; y; nðxÞ; lðyÞÞuðyÞ dCy þ ab ujrig;b ðxÞ ujrig;b ðyÞ uðyÞ dCy C

b¼1

Cb

j¼1

Z ¼ - Gpu ðx; y; nðxÞÞpðyÞ dCy  DðxÞpðxÞ;

x 2 Cp ;

ð35Þ

where u : C ! Rd collects all the boundary displacement field (unknown on Cp and given on Cu). Let      0 P pðyÞ Wp ðyÞ ¼ WðyÞW ¼ 0 Wu ðyÞ U uðyÞ

ð36Þ

be discrete approximations of the displacement field u(y) and of the traction p(y) on Cp and Cu respectively. Here Wu(y) and Wp(y) are d · Nu and d · Np shape-function matrixes modelling the unknown fields on the boundary. The SG-BEM approach is obtained by imposing the fulfilment of (3), on Cu, and (35), on Cp, in an ‘‘average sense’’, using the d · Nu and d · Np test-function matrixes Uu(x) and Up(x). The numerical implementation of the regularized problem (25) follows the same arguments and, in particular, one has to discretize Eq. (35). For, both sides of (35) are multiplied by UTu ðxÞ and integrated on Cp. The only non-standard term is second one on the l.h.s. of (35), which becomes M X

ab

b¼1

¼

3X ðd1Þ j¼1

M X i¼1

ab

Z

UTu ðxÞujrig;b ðxÞ

Cb 3X ðd1Þ j¼1

Z

ujrig;b ðyÞ ðWu ðyÞUÞ dCy dCx Cb

Z

UTu ðxÞujrig;b ðxÞ dCx Cb

Z Cb

! WTu ðyÞujrig;b ðyÞ dCy

U.

ð37Þ

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

975

For any b, let us define the Nu · 3(d  1) matrix hR i R R 3ðd1Þ Bb ¼ Cb WTu u1rig;b dC Cb WTu u2rig;b dC Cb WTu urig;b dC .

ð38Þ

Analogously, the Nu · 3(d  1) matrix Cb is defined as hR i R R 3ðd1Þ Cb ¼ Cb UTu u1rig;b dC Cb UTu u2rig;b dC Cb UTu urig;b dC .

ð39Þ

Consequently, the expression (37) may be written in the following compact form as M X

ab Cb BTb U.

ð40Þ

b¼1

In order to obtain a full-symmetric formulation we set [Uu(x);Up(x)]  [Wu(x);Wp(x)] so that Cb  Bb. The expression (40) can be easily calculated by observing that any IRBD basis vector ujrig;b ð Þ : Cb ! Rd is a linear function. In fact, observe that in general the shape functions matrix Uu(x) may be written in the form (for the 3-D case d = 3) 2 3 2 36 7 0 0 6 uk 7 uk ðxÞ 6 76 7 7 ð41Þ Uu U ¼ 4 0 uk ðxÞ 0 56 6 vk 7; 6 7 0 0 uk ðxÞ 4 wk 5 where (uk, vk, wk) denote the displacement components of the generic kth boundary node. Moreover, common forms of the function uk(x) are such that Z Z Z Z uk ðxÞ dC ¼ 1; uk ðxÞx dC ¼ xk ; uk ðxÞy dC ¼ y k ; uk ðxÞz dC ¼ zk ; ð42Þ C

C

C

C

where (xk, yk, zk) are now the coordinates of the kth node. Therefore, if one selects e.g. vectors of the type (15) as a basis ujrig;b ðxÞ for the IRBD of the bth contour, from (41) and (42) we have 3 2 2 3 7 6 1 0 0 0 z y 0 0 0 zk y k 7 6 1 Z 7 6 6 7 ð43Þ x 5 dC ¼ 6 UTu 4 0 1 0 z 0 1 0 zk 0 xk 7 7. 6 0 Cb 7 6 0 0 1 y x 0 0 0 1 y x 0 5 4 k k Consequently, the matrix Bb in (38) takes the simple form h i Bb ¼ U1rig;b U2rig;b U6rig;b ;

ð44Þ

where T

T

T

ðU1rig;b Þ ¼ ½. . . ; 1; 0; 0; . . .; ðU2rig;b Þ ¼ ½. . . ; 0; 1; 0; . . .; ðU3rig;b Þ ¼ ½. . . ; 0; 0; 1; . . .; ðU4rig;b ÞT ¼ ½. . . ; 0; zk ; y k ; . . .; ðU5rig;b ÞT ¼ ½. . . ; zk ; 0; xk ; . . .; ðU6rig;b ÞT ¼ ½. . . ; y k ; xk ; 0; . . ..

ð45Þ

Obviously, this is not an orthonormal base for the ‘‘discretized’’ IRBD of the bth contour. It could be easily normalized using the Gram–Schmidt algorithm, but we have already remarked in Section 4.1 that such a step is not necessary.

976

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

The numerical approximation of the regularized problem (23) or (25) takes the form 2 U 3 " #    KUP K UP  P P K KU r r M 4 PU 5 P K r ð WÞ ¼ ¼ ¼ F; T P PU P K K þ Bb Bb U U Kr Kr

ð46Þ

b¼1

where, for the reasons discussed at length in Section 4, the matrix Kr is symmetric and non-singular. In particular, the nodal displacement vector Ur , solution of (46), verifies the discretized version of condition (29), i.e. Ujrig;b Ur ¼ 0;

j ¼ 1 . . . 3ðd  1Þ; b ¼ 1; . . . ; M.

ð47Þ

It should be noticed that Kr of (46) is determined by Bb of (44) which, in turn, depends upon the particular basis adopted to define the IRBD of the contours and the type of mesh used for the discretization. In order to recover the real displacement field, we use a discretized version of Eq. (27), which is preferable than (28) for obvious computational reasons. Such a discretization is performed in a similar way and does not present particular difficulties.

6. Numerical examples and conclusions A few representative examples, obtained with a numerical code expressly implementing the proposed regularizing procedure, are now presented. In particular, two-dimensional problems are considered for multiply connected, homogenous, isotropic elastic bodies in generalized plane stress, under various boundary conditions. The numerical integration scheme adopted here for the BEM solution of the hypersingular integral equation is the one proposed in [18]. In the case of ‘‘floating’’ contours, the displacement field has been calculated in two steps. (1) First, the regularized problem (46) is solved. Here, the regularizing term is furnished by expression (40), with Cb  Bb. In particular, the IRBD of any floating contour is referred to the canonical (not-orthonormal) basis (45), which has been used to calculate Bb as per (44). (2) Secondly, the numerical counterpart of Eq. (27) is used to find the displacement of a representative number of points, sufficient to uniquely determine the IRBD of the floating contour. The first example is a hollow square plate whose external boundary is constrained at one side, while a portion of the internal hole is subjected to a pressure p, as represented in Fig. 3. In this case the internal

Fig. 3. Displacement field in a multi-connected body with constrained external contour and Neumann internal boundary.

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

977

contour is clearly a floating contour and, if the regularizing procedure is not adopted, the system of solving equations results irremediably ill-conditioned. With the proposed method, one finds the deformation represented, at a convenient scale, in Fig. 3, where the adopted mesh is indicated with bold face dots. Fig. 4 represents the same body where, now, it is the internal contour to be constrained on one side and the external contour to be under Neumann condition. For the reasons mentioned in Section 3, none of the body contours is a now a floating contours and, indeed, the body deformation has been calculated without any regularizing procedure. When no part of the body is constrained as in Fig. 5, the displacement field can only be determined up to the IRBD of the whole body, but the relative displacement of the contours is uniquely established in the parent linear elasticity problem. However, in the SGBEM both contours are floating contours. In order to use the regularizing procedure, a practical approach consists in regarding the displacement field obtained from the regularized system (46) as the true displacement field for one of the contours and, starting from this, calculating a posteriori the IRBD of the other contour using Eq. (27). Fig. 5a shows the deformation obtained from the solution of (46). Symmetry is clearly spoiled by both translational and rotational IRBDs. Retaining the displacement of the internal contour, the displacement of the external contour is recovered by applying (27), as represented in Fig. 5b. Similarly, without any appreciable difference, one

Fig. 4. Displacement field in a multi-connected body with constrained internal contour and Neumann external boundary.

p1

p1

p2

a

p2

b

Fig. 5. Multi-connected body with all Neumann contours: (a) displacement field after step 1, coinciding with the solution of the regularized system (46); (b) true displacement field after step 2, obtained through Eq. (27).

978

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

could have retained the displacement of the external contour from the solution of (46) and find the displacement of the internal contour using (27). Observe that there is no contradiction between this approach, in which the displacement of the internal contour is calculated from (46) and the IRBD of the external contour is found a posteriori from (27), with the fact that the displacement field can be determined without any need for a regularizing procedure whenever the internal contour is constrained, as in the case of (4). In fact, applying the regularizing term (40) to one contour is not equivalent to constraining that contour. Indeed, as already discussed in Section 4.2 (case 2), whenever some part of a contour is constrained the displacement BIE (3) comes into play, but when the regularizing procedure is applied, only the traction BIE (4) enters the resulting system of equations. In the numerical example of Fig. 5, it has been verified directly that the system (46) remains singular if the regularizing term (40) is applied on one contour only, be it external or internal. The pictures in Fig. 6 refer to the same problem of Fig. 5, with the exception of the number and locations of the nodes used to discretize the internal boundary. In both cases, the deformation has been calculated applying the regularizing term (40) on both contours, regarding as true the displacement of the internal contour and applying Eq. (27) to recover a posteriori the displacement of the external contour. By comparing Figs. 5b and 6b, it should be noted that, whereas the relative displacement between the internal and external contours are the same, the two deformations differ one another by an IRBD of the whole body. This is not surprising if one recalls, from Section 3, that the nodal displacement vector Ur , solution of (46), verifies condition (47), where the vectors Ujrig;b depend not only upon the form of the basis adopted to define the IRBD of the bth contour, but also upon the mesh used for its discretization. Indeed, condition (47) suggests that the IRBD calculated from the regularized system (46) depends, for each contour, upon the location of the centroid of the mesh nodes used to discretize that contour, as the comparison between Figs. 5b and 6b seems to confirm. For a more detailed analysis, the method is verified in a paradigmatic example: the LameÕs classical 2-D problem of an axisymmetrical loaded annular disk, as in Fig. 7. In the numerical approach, the external boundary is approximated with a uniform discretization, whereas the node coordinates on the internal hole are given by h i h i b b x1 ¼ a cos ði=nÞ p þ h ; x2 ¼ a sin ði=nÞ p þ h ; i ¼ 0; . . . ; n; ð48Þ where n is node number, b P 1 allows to bias the nodal intervals, and h rotates the mesh. Notice that by changing the values of the parameters b and h, the centroid of the mesh nodes is displaced.

p1

p1

p2

a

p2

b

Fig. 6. Same problem of Fig. 5 but with a different mesh: (a) displacement field after step 1, coinciding with the solution of the regularized system (46); (b) true displacement field after step 2, obtained through Eq. (27).

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

979

Fig. 7. Annular disk: Geometry and boundary conditions.

Table 1 Annular disk: infinitesimal rigid body displacement of the inner hole and coordinates of its discretizing mesh for various values of parameters b and h of (48) b

h

uxrig1

uxrig2

xg1

xg2

1.0 1.1 1.2 1.3 1.4 1.5 1.5

0 0 0 0 0 0 p/4

0 1.83e2 3.45e2 4.95e2 6.30e2 7.52e2 5.31e2

0 0 0 0 0 0 5.31e2

0 2.79e2 5.29e2 7.55e2 9.59e2 1.14e2 8.10e2

0 0 0 0 0 0 8.10e2

Again, both internal and external contours are floating contours and, in general, after step (1), their relative displacement differs for the well known axisymmetric analytical solution because of an IRBD. The numerical results are reported in Table 1 for various values of the parameter b and h. In particular, the third and fourth columns report the components uxrig1 and uxrig2 of the relative IRBD of the inner hole with respect to the outer boundary in the x1 and x2 direction, respectively, calculated at step (2). In other words, the values uxrig1 and uxrig2 , with their sign changed, represent the components of the IRBM of the inner hole obtained after step (1), which is annihilated at step (2). Besides, the fifth and sixth columns indicate the coordinates xg1 and xg2 of the centroid of the mesh nodes, defined through (48). From the case (b,h)  (0,0), reported on the first row of the table, it is clear that if the mesh centroid coincides with the hole center, the corresponding IRBM is null. On the other hand, if the mesh centroid is moved in the direction of the x1-axis, the hole remains consequently displaced after step (1), as represented in Fig. 8. In particular, notice that whenever uxrig1 6¼ 0, then uxrig1 =xg1 ¼ const. = 6.56e1. Thus, the IRBM of the floating hole calculated at step (1) remains proportional to the coordinates of the centroid of the mesh nodes. This finding is confirmed by the case reported in the last row of Table 1, where

980

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

Fig. 8. Annular disk: displacement field after step 1 (continuous line) and displacement field after step 2 (dashed line).

Fig. 9. Annular disk with rotated mesh. Displacement field after step 1 (continuous line) and displacement field after step 2 (dashed line).

the mesh centroid is displaced now both in the x1 and x2 directions. The displacement fields obtained after steps (1) and after step (2) are graphically compared in Fig. 9. Finally, in order to show a more challenging example, consider the deformation of one of the walls of the famous Notre Dame du Haut Church in Ronchamp (France), Le CorbusierÕs masterpiece. The wall is clamped at the bottom and subjected to distributed loads on some of its sides, as represented in Fig. 10. Now there are sixteen floating contours and the problem is somehow complicated by the presence of sharp corners, but there are no conceptual differences with the previous examples. The displacement field calculated using the aforementioned method is drawn, at a suitable scale, in the same Fig. 10. This latter example confirms, in conclusion, that the numerical approach here proposed provides a robust procedure for solving the most various elasticity problems for bodies with an arbitrary number of floating contours.

F. Freddi, G. Royer-Carfagni / Comput. Methods Appl. Mech. Engrg. 195 (2006) 961–981

981

p

Fig. 10. Displacement field for one of the walls of Le-CorbusierÕs Notre-Dame-du-Haut Church in Ronchamp, under mixed boundary conditions.

References [1] F.J. Rizzo, An integral equation approach to boundary value problems of classical elastostatics, Quart. Appl. Math. 40 (1967) 83–95. [2] C.A. Brebbia, J.C.F. Telles, L.C. Wrobel, Boundary Element Techniques, Springer-Verlag, Berlin, 1984. [3] F. Hartmann, C. Katz, B. Protopsaltis, Boundary elements and symmetry, Ing. Archiv. 55 (1985) 440–449. [4] S. Sirtori, G. Maier, G. Novati, S. Miccoli, A Galerkin symmetric boundary-element method in elasticity: formulation and implementation, Int. J. Numer. Methods Engrg. 35 (1992) 255–282. [5] W. Hackbusch, Integral Equations: Theory and Numerical Treatment, Birkhauser Verlag, Basel, 1995. [6] J.H. Kane, C. Balarishna, Symmetric Galerkin formulation employing curved elements, Int. J. Numer. Methods Engrg. 36 (1993) 2157–2187. [7] L.M. Milne-Thomson, Theoretical Hydrodynamics, second ed., The Macmillan Company, New York, 1950. [8] A. Frangi, G. Novati, Symmetric BE method in two-dimensional elasticity: evaluation of double integrals for curved elements, Comput. Mech. 19 (1996) 58–68. [9] J.J. Pe´rez-Gavila´n, M.H. Aliabadi, Symmetric Galerkin BEM for multi connected bodies, Commun. Numer. Methods Engrg. 17 (2001) 761–770. [10] J.J. Pe´rez-Gavila´n, M.H. Aliabadi, A symmetric Galerkin BEM for multi-connected bodies: a new approach, Engrg. Anal. Boundary Elem. 25 (2001) 633–638. [11] G. Veˆrchery, Re´gularisation du syste`me de lÕe´quilibre des structures e´lastiques discre`tes, C.R. Acad. Sci. Paris, 311, Se´rie II (1990) 585–589. [12] T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Prentice Hall, Englewood Cliffs, 1987. [13] D. Kosloff, G.A. Frazier, Treatment of hourglass patterns in low order finite element codes, Numer. Anal. Meth. Geomech. 2 (1978) 55–72. [14] C. Polizzotto, An energy approach to the boundary element method. Part I, Elastic solids, Comput. Methods Appl. Mech. Engrg. 69 (1988) 167–184. [15] M. Bonnet, G. Maier, C. Polizzotto, Symmetric Galerkin boundary element method, Appl. Mech. Rev. 51 (1998) 669–704. [16] M. Guiggiani, Hypersingular formulation for boundary stress evaluation, Engrg. Anal. Boundary Elem. 13 (1994) 169–179. [17] M. Guiggiani, Formulation and numerical treatment of boundary integral equations with hypersingular kernels, in: V. Sladek, J. Sladek (Eds.), Singular Integrals in Boundary Element Methods, Computational Mechanics Publications, Southampton, 1998. [18] A. Aimi, M. Diligenti, G. Monegato, New numerical integration schemes for the BEM solution of hypersingular integrals equations, Int. J. Numer. Methods Engrg. 45 (1999) 1807–1830.