- Email: [email protected]
Symmetric boundary element method versus finite element method
Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367 www.elsevier.com/locate/cma
Symmetric boundary element method versus finite element method T. Panzeca a
a,*
, F. Cucco 1, S. Terravecchia
b
Dipartimento di Ingegneria Strutturale e Geotectnica, Universit a degli Studi di Palermo, Viale delle Scienze, 90128 Palermo, Italy b Dipartimento di Strutture, Universit a della Calabria, Arcavacata di Rende, Cosenza I-87030, Italy Received 6 August 2001; received in revised form 14 January 2002
Abstract The paper examines the effectiveness of the symmetric boundary element formulation when the continuum body is subdivided into large elements called macro-elements. The approach proposed combines a strong reduction of variables with an elastic solution close to the real response. Indeed, if the displacement method is used, this approach permits one to determine for every macro-element a relationship connecting the weighted traction vector defined on the sides of the interface boundary with the node displacement vector of the same boundary and with the external action vector. Such a strategy is very similar to that followed through the finite element method, but with the advantages of having the punctual satisfaction of the equilibrium and compatibility by using the fundamental solution, and of employing large elements whose discretization is performed on the boundary only. Some examples are shown using a general computing programme developed by the present authors. 2002 Elsevier Science B.V. All rights reserved. Keywords: Symmetric boundary element method; Boundary element; Substructuring
1. Introduction In all fields of applied mechanics and engineering there is a need to work out an analysis having the capacity to combine the employment of a reduced number of variables with solutions which are closer and closer to the real response. The finite element method (FEM), using elements generically defined as having high performance, has this aim and formulations like assumed stress (Pian, Casciaro), assumed strain (Wilson, Bergan, Hughes, Felippa) and enhanced assumed strain (Simo) may be placed in this ambit. The symmetric boundary element method (SBEM) proves to be preferred for its characteristics in order to reach the two objectives (few variables, closer response) both because it performs using large elements having boundary variables only, and because the use of the fundamental solutions guarantees the punctual satisfaction of the field equilibrium and the compatibility equations. On the boundary only the statical and
*
Corresponding author. E-mail address: [email protected] (T. Panzeca). 1 External collaborator.
0045-7825/02/$ - see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 2 ) 0 0 2 3 9 - 6
3348
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
kinematical conditions are written in the weighted form, and as a consequence the goodness of the solution is the result of the discretization performed on the boundary and of the type of variable modelling. Some computational problems arising in some coefficients of the solving system because of the presence of double integrals having strong singularity or hypersingularity in the definition domain have been solved by using analytical or numerical techniques. As a consequence the SBEM has had a large development in several fields of mechanics. The pertinent references are found in paper [1]. Some difficulties have arisen in the analysis of a continuum body having zone-wise different physical and geometrical characteristics subjected to external actions, but recently have been overcome by suitable strategies used in an appropriate substructuring. The latter idea was introduced by Maier et al. [2] in the ambit of variational elastic analysis in dynamics and resumed by Steinbach and Wendland [3]. Subsequently Gray and Paulino [4] used this expedient for potential problems including an example regarding heat transfer. Layton et al. [5] at first proposed a formulation which divides the body into macro-zones, each of which is governed by boundary quantities, and which after a condensation process gives rise to a system having only few asymmetric blocks. In a subsequent paper the same authors [6] presented a truly symmetric method for a two-zone elastic problem characterized by interfacial and noninterfacial unknowns. In [7,8] Panzeca et al. introduced a procedure which defines for each macro-element, obtained by appropriate substructuring, a matricial equation connecting the interfacial boundary quantities only. A subsequent assembly by means of weighted regularity conditions written on the interfaces among macroelements gives rise to a reduced equation system obtained as the solution of a problem with displacement unknowns, force unknowns or both variables. Inside the displacement approach the strategy presented in this paper proves to obtain a basic working matrix for every macro-element of domain X thought of as embedded in the infinite domain X1 having the same characteristic as X. The use of this matrix permits us to generate for every macro-element both the algebraic operators depending on the boundary type and the known vectors depending on the external actions on the element. Indeed, once the boundary type is established, i.e. constrained C1 , free C2 , interfacial C0 among macro-elements, some possible stiffness matrices can be defined having as a boundary • • • •
C0 , where the boundary is totally embedded in contiguous macro-elements; C0 [ C1 , where a portion of the boundary is constrained on the ground; C0 [ C2 , where a portion of the boundary is free; C0 [ C1 [ C2 , where all types of boundary coexist.
These stiffness matrices connect the weighted tractions defined on the sides of the interface boundary C0 to the node displacements of the same boundary. In addition, the same basic matrix permits us to derive for every macro-element the load vector including the layer forces acting on the free boundary C2 and the displacements imposed on the constraint C1 . The load vector including the body forces and the linear and volume kinematic discontinuities are obtained in the classical way. This load vector is defined on the sides of the boundary C0 where the stiffness matrix coefficients are computed. This strategy permits us to define the macro-elements of SBEM in a similar way to the finite elements of FEM, but has very high performance. As a consequence, the usual assembly rules involve a very small number of variables, i.e. the displacements of the macro-element interface boundary, and give this approach strong effectiveness because of the reduced dimension of the analysis problem. A general computing code was worked out to show the benefits of the approach proposed and the results of the examples dealt with are discussed below.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3349
2. Linear elastic analysis by substructuring Let us consider a homogeneous two-dimensional solid body B with domain X and boundary S referred to a Cartesian coordinate system ðx; yÞ. This body is subjected to imposed displacements u on the boundary portion Su restraining rigid motion, to known forces f on the free complementary boundary portion Sf ¼ S Su and to body forces b and imposed strains # in X. In this body the physical characteristics (Young’s modulus E and Poisson coefficient m) and the geometry (thickness s) are zone–wise variables. In the infinitesimal displacement hypothesis, the elastic response to the known external actions can be obtained in terms of displacements u on Sf , of forces f on Su , of displacements u, stresses r and strain in X by utilizing the SBEM. To this aim an appropriate substructuring must be employed considering the different zone–wise characteristics in the body. In Fig. 1(a) the body is subjected to external actions, and in Fig. 1(b) a domain subdivision into macro-elements is performed. For each of them an appropriate relation between mechanical and kinematical quantities defined on the boundary may be written, whose form depends on the type of approach to be employed. Indeed, in [8] the authors have shown that a substructuring analysis can be made by using approaches based on the displacements, the forces, or the mixed-value quantities. In the present paper the displacement approach will be utilized, but the same results can easily be obtained using alternative approaches. The substructuring in Fig. 1(b) involves the characterization of macro-elements whose boundaries can be constrained on C1 and free on C2 or can be interfacial on C0 among neighbouring macro-elements. In this way the macro-element A is characterized by a boundary C0 , whereas the macro-elements B, C and D, show the mixed boundaries C1 [ C0 , C2 [ C0 and C1 [ C2 [ C0 , respectively. The specification of different types of macro-elements involves all the possible cases to be characterized in substructuring a body. Let us consider in Fig. 2 the ith macro-element of domain X embedded in the infinite domain X1 having the same characteristics ðE; m; sÞ as the macro-element considered. The domain boundary may be thought of as a boundary C of X or as a boundary Cþ of X1 n X. At this stage no distinction is made for the constrained, free or interfacial boundary type. This macro-element is subjected to the layered forces f and to the double-layered displacement discontinuities v ¼ u, both on the boundary C, and to the body forces b, to thermal-like strains # and cracks d inside the domain, the symbol bar specifying that the corresponding quantities are known.
Fig. 1. (a) Two-dimensional body subjected to external actions, (b) appropriate substructuring of the domain.
3350
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 2. Macro-element embedded in the infinite domain.
According to the symbols used in [1], the Somigliana Identities (SI) can be written in the following way: u ¼ u½f þ u½u þ ^ u½ b; #; d ; ð1aÞ t ¼ t½f þ t½u þ ^t½ b; #; d
ð1bÞ
which represent the response to all the actions in X1 . In the latter the following positions are set Z Z u½f ¼ Guu f; t½f ¼ Gtu f; C
u½u ¼
ð2a; bÞ
C
Z
Gut ðuÞ;
t½u ¼
Z
C
^ u½ b; #; d ¼
Z
Guu bþ
Z
X
^t½ b; #; d ¼
Gtt ðuÞ;
ð2c; dÞ
C
Z
Gur # þ
Z
X
Gtu bþ
X
Z
Gtr # þ X
Gut d;
ð2eÞ
X
Z
Gtt d:
ð2fÞ
X
If the response is computed on the boundary Cþ in Fig. 2, one obtains u½ b; #; d ; uþ ¼ u½f þ uPV ½u 1u þ ^ 2
b; #; d
tþ ¼ tPV ½f þ t½u 12f þ ^t½
ð3aÞ ð3bÞ
where the apex PV appended to a symbol means that the corresponding quantities are computed taking the related singular integral as the Cauchy principal value. Imposing the punctual condition uþ ¼ 0, tþ ¼ 0 on the boundary Cþ permits one to obtain the following mixed boundary integral equations of the macro-element: u½ b; #; d ¼ 0; ð4aÞ u½f þ uPV ½u 1u þ ^ 2
b; #; d ¼ 0: tPV ½f þ t½u 12f þ ^t½ f, ðuÞ being the unknown function vectors defined on all the boundary of the body.
ð4bÞ
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3351
Fig. 3. Boundary forces and displacements of the nodes.
Let us execute the boundary discretization of the macro-element into the boundary elements and introduce the shape functions Wf , Wu to model the layered and double-layered boundary quantities f ¼ Wf F;
u ¼ Wu U
ð5a; bÞ
where F and U are vectors of unknown nodal values, forces and displacements, described in Fig. 3. In every boundary element the force vector F is defined through four components modelling the forces along the boundary element, whereas the displacement vector U is made up of two components per node modelling the displacements of the boundary elements contiguous to the node considered. We now introduce the modelling assumed by Eqs. (5a,b) in Eqs. (4a) and (4b) and perform the weighting of the system coefficients using the same shape functions in accordance with the Galerkin strategy. The following equations, written in weighted form, are obtained Z Z Z Z Z Z Z 1 þ Wf u Wf Guu Wf F þ Wf Gut Wu ðUÞ þ Wf Wu ðUÞ þ Wf ^u ¼ 0; ð6aÞ 2 C C C C C C C Z Z Z Z Z Z Z 1 Wu tþ Wu Gtu Wf F Wu Wf F þ Wu Gtt Wu ðUÞ þ Wu^t ¼ 0: ð6bÞ 2 C C C C C C C The latter equations written in compact form become Wþ Auu F þ Auf ðUÞ þ Cuf ðUÞ þ c W X ¼ 0;
ð7aÞ
bX ¼ 0 Pþ Afu F Cfu ðFÞ þ Aff ðUÞ þ P
ð7bÞ
with obvious symbol meaning, or also bX ¼ 0 ðA þ CÞX þ L where one sets Auu A¼ Afu
Auf ; Aff
ð8Þ
0 C¼ Cfu
Cuf ; 0
F X¼ ; U
b¼ L
c WX : bX P
ð9a–dÞ
The matrix A is symmetric. Indeed, the following positions are valid: Auu ¼ ATuu ;
Aff ¼ ATff ;
Afu ¼ ATuf
ð10a–cÞ
in that the matrix coefficients are obtained by using the fundamental solutions Ghk with h; k ¼ u; f , employing the symmetry properties [1], completed using the shape functions in accordance with the Galerkin approach. Otherwise the matrix C is hemisymmetric ðCfu ¼ CTuf Þ, Cuf being a block diagonal matrix. The sum matrix ðA þ CÞ is nonsymmetric and singular. The singularity arises because in the column block related to the vector ðUÞ, three columns are a linear combination of the remaining columns, the macro-element being subjected to a rigid motion in its plane. This circumstance takes on particular
3352
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
significance because it can be employed as a check on the column-block coefficients. Indeed, when we impose in the infinite domain a displacement discontinuity vector ðUÞ on C such as to cause a rigid motion of the domain X, we obtain Auf þ Cuf 0 ðUÞ ¼ : ð11Þ Aff 0 If this equality is satisfied, then the matrix shows coefficients computed correctly. At this stage other types of check cannot be introduced. The matrices ðA þ CÞ, written for all the macro-elements arising from the substructuring employed, take on a fundamental role in the analysis of the entire body. They are obtained by assuming the boundaries as free, but subjected to layered forces F and double-layered displacement discontinuities ðUÞ, without specifying whether the quantities are known or not. Starting from Eq. (8) it is possible to obtain an equation connecting the boundary mechanical and kinematical quantities, known and unknown, for every macro-element. 3. Symmetric operator specification In the ambit of the displacement method via SBEM, the analysis is dealt with by substructuring the body after a preliminary study made on macro-elements like those marked by capital Latin letters A, B, C and D in Fig. 1(b). For each macro-element a stiffness matrix and a load vector will be defined and computed. 3.1. Macro-elements like A This case regards a macro-element totally delimited by a boundary like C0 , i.e. a boundary interfacing with the neighbouring macro-elements. Eq. (8) is reformulated introducing the pedix 0 in all the components of the system c Au0u0 Au0f 0 þ Cu0f 0 F0 0 Wþ W X0 0 þ b ¼ ð12Þ A C A U 0 Pþ f 0u0 f 0u0 f 0f 0 P 0 X0 0 which may be rewritten in the following way: c Wþ 0 Au0u0 F0 þ ðAu0f 0 þ Cu0f 0 ÞðU0 Þ þ W X0 ¼ 0;
ð13aÞ
b Pþ 0 ðAf 0u0 Cf 0u0 ÞF0 þ Af 0f 0 ðU0 Þ þ P X0 ¼ 0:
ð13bÞ
The latter equation may be substituted by the following: b X0 P0 ¼ ðAf 0u0 þ Cf 0u0 ÞF0 þ Af 0f 0 ðU0 Þ þ P
ð14Þ
which defines the weighted traction vector on the actual boundary of the macro-element. It is to be noted that this change involves a variation in sign for the free term matrix Cf 0u0 which, as a consequence, becomes positive. The system written in compact form becomes þ c X0 ðAu0f 0 þ Cu0f 0 Þ Au0u0 F0 on Cþ W0 0 0 ¼ þ W ð15Þ b X0 ðAf 0u0 þ Cf 0u0 Þ Af 0f 0 U0 on C0 P0 P where the matrix is symmetric. We obtain F0 by using Eq. (13a), i.e. h i c F0 ¼ A1 u0u0 ðAu0f 0 þ Cu0f 0 Þð U0 Þ þ W X0
ð16Þ
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3353
and substitute it in Eq. (14), thus obtaining the following equation: b X0 P0 ¼ D00 U0 þ L
ð17Þ
where D00 ¼ ðAf 0u0 þ Cf 0u0 ÞA1 u0u0 ðAu0f 0 þ Cu0f 0 Þ Af 0f 0 ;
ð18aÞ
b X0 ðAf 0u0 þ Cf 0u0 ÞA1 c b X0 ¼ P L u0u0 W X0 :
ð18bÞ
Eq. (17) is representative of any macro-element like A and connects the weighted tractions acting on the sides of the discretized interface boundary C0 to the node displacements of the same boundary and to the domain external actions. The matrix D00 is singular because the macro-element can be subjected to a rigid motion in its plane. If the components of the vector U0 are chosen in such a way as to define a rigid motion without domain loads, the matricial product in Eq. (17) must become zero, i.e. D00 U0 ¼ 0
ð19Þ
and this check is additional the one in Eq. (11). Indeed, the first check is used to verify the column block related to the displacement discontinuity, and the second check can be utilized to verify the matrix D00 built using the block matrices characterizing the free body embedded in X1 . 3.2. Macro-elements like B The macro-element is delimited by a constrained boundary C1 and by an interface boundary C0 among neighbouring macro-elements. Eq. (8) is reformulated after rearrangement of the row blocks and the column blocks in relation to the boundary type. In this case one has 2 32 3 2 3 2 3 c Au1u1 Au1u0 Hu1f 1 Au1f 0 F1 0 Wþ W 1X 1 6 Au0u1 Au0u0 Au0f 1 Hu0f 0 76 F0 7 6 c 7 607 Wþ W 0 6 76 7 6 0X 7 ¼ 6 7 ð20aÞ 4 Ef 1u1 Af 1u0 Af 1f 1 Af 1f 0 54 U1 5 þ 4 P b 1X 5 4 0 5 Pþ 1 þ b 0X 0 U0 P0 Af 0u1 Ef 0u0 Af 0f 1 Af 0f 0 P where on sets Ef 1u1 ¼ Af 1u1 Cf 1u1 ;
Hu1f 1 ¼ Au1f 1 þ Cu1f 1 ;
ð20b; cÞ
Ef 0u0 ¼ Af 0u0 Cf 0u0 ;
Hu0f 0 ¼ Au0f 0 þ Cu0f 0
ð20d; eÞ
and where among the components of the vector which defines the boundary quantities the displacement discontinuity subvector ðU1 Þ is assigned on C1 . The latter equations may be rewritten in explicit form as follows: c Wþ 1 Au1u1 F1 þ Au1u0 F0 þ Hu1f 1 ðU1 Þ þ Au1f 0 ðU0 Þ þ W 1X ¼ 0;
ð21aÞ
c Wþ 0 Au0u1 F1 þ Au0u0 F0 þ Au0f 1 ðU1 Þ þ Hu0f 0 ðU0 Þ þ W 0X ¼ 0;
ð21bÞ
b Pþ 1 Ef 1u1 F1 þ Af 1u0 F0 þ Af 1f 1 ðU1 Þ þ Af 1f 0 ðU0 Þ þ P 1X ¼ 0;
ð21cÞ
b Pþ 0 Af 0u1 F1 þ Ef 0u0 F0 þ Af 0f 1 ðU1 Þ þ Af 0f 0 ðU0 Þ þ P 0X ¼ 0:
ð21dÞ
3354
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Eq. (21d) may be substituted by the following equation: b 0X P0 ¼ Af 0u1 F1 þ ðAf 0u0 þ Cf 0u0 ÞF0 þ Af 0f 1 ðU1 Þ þ Af 0f 0 ðU0 Þ þ P
ð22Þ
which defines the weighted traction vector on the boundary C0 of the macro-element. It is to be noted that here too this change involves a variation in sign of the free term matrix Cf 0u0 which, as a consequence, becomes positive. A block rearrangement of the matrix and of the vectors of Eq. (20a) is performed, taking into account of Eq. (22) and of the known displacement discontinuity vector (U1 ), so the system (20a) may be reformulated as follows: ð23Þ with the addition of the following equation: b Pþ 1 0 ¼ Ef 1u1 F1 þ Af 1u0 F0 þ Af 1f 0 ðU0 Þ þ Af 1f 1 ðU1 Þ þ P 1X
ð24Þ
and where Hf 0u0 ¼ HTu0f 0 . It is to be noted that the matrix in Eq. (23) is symmetric. Let us rewrite Eq. (23) in compact form, i.e. 0 ¼ BX þ B0 ðU0 Þ þ c Wu þ c WX;
ð25aÞ
bu þ P bX P0 ¼ BT0 X þ B00 ðU0 Þ þ P ð25bÞ T T T where the vector X ¼ F1 F0 has been introduced and where the remaining terms take on an obvious meaning. The vector X is obtained by using Eq. (25a) X ¼ B1 ½B0 ðU0 Þ þ c Wu þ c WX
ð26Þ
and is substituted in Eq. (25b) giving the following equation bu þ L bX P0 ¼ D00 U0 þ L ð01Þ
ð27Þ
where ð01Þ
D00 ¼ BT0 B1 B0 B00 ;
ð28aÞ
bu ¼ P b u BT B1 c L Wu; 0
ð28bÞ
b X BT B1 c bX ¼ P WX: L 0
ð28cÞ
In this case too, Eq. (27) connects the weighted traction vector defined on the sides of the discretized interface boundary C0 to the node displacement vector of the same boundary, to the domain external actions and to the load vector (U1 ) related to the constraint, the two latter vectors also being computed on the boundary elements of C0 . ð01Þ In this case the matrix D00 is nonsingular: indeed, a portion of the boundary of the macro-element is constrained on the ground. If the components of the vectors U0 and U1 are chosen in such a way as to define a rigid motion, without domain external actions, we have
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367 ð01Þ bu ¼ 0 D00 U0 þ L
3355
ð29Þ
b u in Eqs. (28a)– and L employed as second check of the coefficients of the submatrices which define (28c). Eq. (24) written on Cþ 0 is not employed to obtain the relation connecting the weighted tractions to the displacements of the boundary C0 . It may be of use as a check of the approximation of the numerical solution: the weighted traction value would have to be zero on each boundary element of Cþ 1. ð01Þ D00
3.3. Macro-elements like C The macro-element is delimited by a free boundary C2 and by an interface boundary C0 with neighbouring macro-elements. The previous Eq. (8) is reformulated after rearrangement the row blocks and the column blocks according to the boundary type. One has 2 32 3 2 3 2 3 c Au2u2 Au2u0 Hu2f 2 Au2f 0 0 Wþ F2 W 2X 2 7 6 7 6 7 6 7 þ6 W0 6 Au0u2 Au0u0 Au0f 2 Hu0f 0 76 F0 7 6 c W 7 607 ð30aÞ 76 7 þ 6 0X 7 ¼ 6 7 þ 6 b 2X 5 4 0 5 P2 4 Ef 2u2 Af 2u0 Af 2f 2 Af 2f 0 54 U2 5 4 P b 0X Af 0u2 Ef 0u0 Af 0f 2 Af 0f 0 0 Pþ U0 P 0 where, in addition to the position (20d,e), one sets Ef 2u2 ¼ Af 2u2 Cf 2u2 ;
Hu2f 2 ¼ Au2f 2 þ Cu2f 2
ð30b; cÞ
and where among the components of the vector which defines the boundary quantities the force vector F2 is assigned on the C2 . The latter equations may be rewritten in explicit form as follows: c Wþ 2 Au2u2 F2 þ Au2u0 F0 þ Hu2f 2 ðU2 Þ þ Au2f 0 ðU0 Þ þ W 2X ¼ 0;
ð31aÞ
c Wþ 0 Au0u2 F2 þ Au0u0 F0 þ Au0f 2 ðU2 Þ þ Hu0f 0 ðU0 Þ þ W 0X ¼ 0;
ð31bÞ
b Pþ 2 Ef 2u2 F2 þ Af 2u0 F0 þ Af 2f 2 ðU2 Þ þ Af 2f 0 ðU0 Þ þ P 2X ¼ 0;
ð31cÞ
b Pþ 0 Af 0u2 F2 þ Ef 0u0 F0 þ Af 0f 2 ðU2 Þ þ Af 0f 0 ðU0 Þ þ P 0X ¼ 0:
ð31dÞ
The use of the procedure adopted for the macro-element like C, and therefore the writing of Eq. (31d) on C0 instead of Cþ 0 and the arrangement of the block matrices and of the block vectors in Eq. (30a) permit us to reformulate the problem as follows: ð32Þ with the following equation: c Wþ 2 Au2u0 F0 þ Hu2f 2 ðU2 Þ þ Au2f 0 ðU0 Þ þ Au2u2 ðF2 Þ þ W 2X ¼ 0:
ð33Þ
The matrix in Eq. (32) is symmetric. The block equation may be written in the following compact form: 0 ¼ BX þ B0 ðU0 Þ þ c Wf þ c WX;
ð34aÞ
bf þ P bX P0 ¼ BT0 X þ B00 ðU0 Þ þ P
ð34bÞ
3356
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
where the vector XT ¼ FT0 UT2 has been introduced and where the remaining terms take on an obvious meaning. Obtaining the vector X by using Eq. (34a) and substituting it in Eq. (34b) the following equation can be attained: ð02Þ bf þ L bX P0 ¼ D00 U0 þ L
ð35Þ
where ð02Þ
D00 ¼ BT0 B1 B0 B00 ;
ð36aÞ
b f BT B1 c bf ¼ P Wf ; L 0
ð36bÞ
b X BT B1 c bX ¼ P WX: L 0
ð36cÞ
In this case too the final expression (35) connects the weighted traction vector on the discretized boundary elements of C0 to the node displacement vector, whose components are defined on the same boundary, to the domain external actions and to the force vector F2 related to the free boundary, both latter vectors also being computed on C0 . ð02Þ The matrix D00 is symmetric and singular: indeed, no portion of the boundary of the macro-element is constrained on the ground. As a consequence, when the vector U0 is chosen in such a way as to define a rigid motion, we must have ð02Þ
D00 U0 ¼ 0:
ð37Þ
3.4. Macro-elements like D All the boundary types are present, macro-elements. The block form of Eq. 2 Au1u1 Au1u2 Au1u0 Hu1f 1 Wþ 1 þ6 W2 6 Au2u1 Au2u2 Au2u0 Au2f 1 6 6 Au0u2 Au0u0 Au0f 1 Wþ 0 6 Au0u1 þ 6 P1 6 Ef 1u1 Af 1u2 Af 1u0 Af 1f 1 6 4 Af 2u1 Ef 2u2 Af 2u0 Af 2f 1 Pþ 2 þ P0 Af 0u1 Af 0u2 Ef 0u0 Af 0f 1
i.e. C1 constrained, C2 free, C0 interfacial among neighbouring (8) becomes the following: 32 3 2c 3 2 3 Au1f 2 Au1f 0 0 F1 W 1X 7 6 7 6 7 7 6 c 7 6 Hu2f 2 Au2f 0 76 F2 7 W 2X 7 6 0 7 76 7 6 7 607 c 7 Au0f 2 Hu0f 0 76 F0 7 6 W 0X 7 6 7 76 7þ6 ð38Þ 7¼6 6 6 7 7 6 b 1X 7 6 0 7 Af 1f 2 Af 1f 0 76 U1 7 6 P 7 7 6 6 7 76 7 b 2X 7 Af 2f 2 Af 2f 0 54 U2 5 6 5 405 4P b 0X Af 0f 2 Af 0f 0 0 U0 P
where the positions (20b,c), (20d,e) and (30b,c) are valid. The evaluation of the vector P0 on C0 instead of þ the vector Pþ 0 on C0 and an appropriate arrangement produce the following matricial equation:
ð39Þ
with the equations c Wþ 2 0 ¼ Au2u1 F1 þ Au2u0 F0 þ Hu2f 2 ðU2 Þ þ Au2f 0 ðU0 Þ þ Au2u2 F2 þ Au2f 1 ðU1 Þ þ W 2X ;
ð40aÞ
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
b Pþ 1 0 ¼ Ef 1u1 F1 þ Af 1u0 F0 þ Af 1f 2 ðU2 Þ þ Af 1f 0 ðU0 Þ þ Af 1u2 F2 þ Af 1f 1 ðU1 Þ þ P 1X :
3357
ð40bÞ
Eq. (39) is rewritten in compact form with an obvious meaning of the symbols 0 ¼ BX þ B0 ðU0 Þ þ c Wf þ c Wu þ c WX;
ð41aÞ
bf þ P bu þ P bX P0 ¼ BT0 X þ B00 ðU0 Þ þ P
ð41bÞ
and therefore ð012Þ bb þ L bX P0 ¼ D00 U0 þ L
ð42Þ
where ð012Þ
D00
¼ BT0 B1 B0 B00 ;
ð43aÞ
bf þ L b u ¼ ðP b f BT B1 c b u BT B1 c bb ¼ L Wf Þ þ ðP W u Þ; L 0 0
ð43bÞ
b X BT B1 c bX ¼ P WX: L 0
ð43cÞ
ð012Þ D00
The matrix is symmetric and nonsingular: indeed, a portion of the boundary is constrained on the ground. A rigid motion involves the following check: b u ¼ 0: D00 U0 þ L ð012Þ
ð44Þ
4. Solving system In order to obtain the elastic response of a structure subjected to known external actions, a two-dimensional continuum body serving as an explanation is considered in Fig. 4. As a consequence of the different physical ðE; mÞ and geometrical (s) characteristics it is subdivided into three macro-zones specified by the capital Latin letters A, B and C. The macro-elements A, B and C belong to macro-elements of type D, C and C respectively, and each of them is characterized by a relationship between the weighted tractions and the node displacements, both quantities defined on C0 , as can be seen in Eqs. (17), (27), (35) and (42) in the previous section. The macro-element interfaces are called by the small Latin letters a, b and c and for simplicity the double indices 00 characterizing the interface and the apices characterizing the macro-element type are suppressed. In each stiffness matrix D employed a rearrangement of the block rows and the block columns is made in order to assemble the quantities associated with every boundary dividing the vector P in function of the connection points among macro-elements. For the element A we will have
Fig. 4. Body subdivided into three macro-zones A, B and C.
3358
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
2
32 A Daa PAa 6 A 76 A 4 Pc 54 Dca PA0 DA0a
DAac DAcc DA0c
32 A 3 2 A 3 b P Ua DAa0 a 7 6 7 6 bA 7 DAc0 54 UAc 5 þ 4 P c 5 b DA00 UA0 PA
ð45Þ
0
or in compact form b A: P A ¼ DA UA þ P
ð46Þ A
UAa
UAc
The displacement vector U collects the node displacements of the boundary a, those of the boundary c, and UA0 those of the common node with both the boundary elements a and b. An equation like (46) can be written for the macro-element B and C. In this way one obtains a global relation connecting the weighted tractions to the displacements of the nodes concerning the same interface boundary b P ¼ DU þ P
ð47Þ
where the following definitions are valid: bT ¼ P b AT PT ¼ PAT PBT PCT ; P UT ¼ UAT
UBT
UCT ;
b CT ; P
b BT P
D ¼ diag DA
DB
ð48a; bÞ
DC :
ð48c; dÞ
The system compatibility equation is obtained by introducing the displacement vector of the interface nodes and that of the nodes which are common vertices with the macro-elements, i.e. uT ¼ uTa uTb uTc uT0 : ð49Þ As a consequence U ¼ Zu
ð50Þ
where Z is the topological matrix. The regularity condition of the weighted tractions takes on the form ZT P ¼ 0:
ð51Þ
The use of Eqs. (47), (50), and (51) gives rise to the following solving system: Ku þ ^f ¼ 0
ð52Þ
where for the structure being examined
ð53aÞ
and h ^f T ¼ ðZT P b ÞT ¼ ð P bA þ P b B ÞT a a
bB þ P b C ÞT ðP b b
bA þ P b C ÞT ðP c c
i bA þ P bB þ P b C ÞT : ðP 0 0 0
ð53bÞ
Here too the rigid motion technique can be employed in order to check the correctness of the assembly performed, but it is to be noted that this technique must involve both the vector u and the equivalent force vector ^f through the displacements imposed on the constraints present.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3359
Using Eq. (52) gives the displacement vector u and in succession we obtain all the boundary quantities of each macro-element and therefore through the SI the elastic response in the domains.
5. Numerical applications The first example to be shown includes two types of check, one to prove the convergence, the other to verify the goodness of the solution when the substructuring changes its geometry. In both the checks the solution chosen for comparison is obtained by considering the body as a single macro-element and performing a dense discretization on the boundary or, where it is possible, by making reference to close form solutions. The check is made on both the displacement and the stress at meaningful points in the body. The convergence check is obtained by • subdividing the body into macro-elements utilizing a sparse discretization of the boundary, and • preserving the geometry of the macro-elements, but increasing the number of the boundary elements on the sides in the same ratio ð2; 3; 4; . . .Þ in each macro-element. The check on the goodness of the solution contemplating the change in the substructure geometry of the macro-elements is obtained by • subdividing the body into macro-elements whose shapes show several unevennesses and performing a pre-established boundary discretization, and • preserving the number of boundary elements at the sides, but changing the geometry of the macro-elements inside the domain towards more regular shapes. In the first two examples the physical and geometrical characteristics of each macro-element are assumed equal, and this is done to permit the checks, but the approach employed and the computing programme make it possible to have different values for the Young modulus E, the Poisson coefficient m and for the thickness s of each macro-element. The modelling of the displacements and of the tractions on the boundary elements is assumed linear. 5.1. First example: cantilever beam Let us consider a beam (l h ¼ 10 4) constrained at the first extreme against rigid motion and subjected to a unit shear force distributed by a parabolic law at the tip (Fig. 5). Let us assume E ¼ 1, m ¼ 0:3, and the thickness s ¼ 1. The vertical displacement w of the point in the corner of the tip (x ¼ 10, y ¼ 2) and the axial stress rx at the upper extreme of the constrained element 54–55 (x ¼ 0, y ¼ 1) are chosen as quantities to be checked. At first the beam is considered as a single domain whose boundary is discretized by introducing a denser discretization made up by 56 boundary elements and characterized by 47 free nodes. The values obtained at the points chosen by the classic symmetric BE analysis using the discretization mentioned are taken as comparison quantities in the computations below. 5.1.1. Convergence check Let us assume a domain subdivided by six macro-elements, as in Fig. 6a, characterized by nodes, these being the vertices of the elements. This domain, thus subdivided, is assumed unchanged in the test. The
3360
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 5. Cantilever beam subjected to a unit shear force at the tip.
Fig. 6. Cantilever beam subdivided into six MEs: (a) with a single boundary element for every side, (b) with four boundary elements for every side.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3361
constrained and the tip boundaries will have an unchanged discretization, each made up by eight elements, whereas each side of the macro-element boundaries will be subdivided into 2, 3 and 4 elements (see Fig. 6b). The analysis is performed using the displacement method through the interface node displacements, therefore the progressive subdivision of the boundary sides of each macro-element leads to the following numbers of the free nodes: ME ME ME ME ME ME
1: 2: 3: 4: 5: 6:
13–18–23–28. 4–8–12–16. 7–14–21–28. 6–12–18–24. 9–18–27–36. 13–18–23–28.
When the macro-element boundary discretization increases, the number of the interfacial nodes involved in the equation system increases in the following way: 15–29–43–57. It is to be noted that the analysis of the beam subdivided into six macro-elements using few interfacial variables permits us to obtain displacements and tractions which are very close to the values obtained using the classic symmetric approach in the beam thought of as a single element (see Table 1). The reference closed form values in Table 1 are obtained by considering as a hypothesis a tangential stress having a parabolic form in the constrained portion. The true tangential stress distribution shows a different shape and as a consequence the response of both the displacements and the stresses is different from the theoretical one. A variable response is obtained for the stress: indeed, the value is a floating one around the reference solution. 5.1.2. Check of the goodness of the solution The substructuring is made up by six macro-elements each having very irregular geometry and the boundary sides each constituting a single boundary element, with the exception of the constraint and tip boundaries where a dense discretization (eight elements) is made. A subsequent modification of the macro-elements towards more regular geometries is performed and for every change it is to be noted that the solution in terms of displacements becomes close to the values obtained by the classic SBEM without substructuring chosen as a reference. In this case too the stress value is a floating one around the reference solution (Fig. 7, Table 2).
Table 1 Comparison between the displacement and stress values obtained by substructuring SBE analysis and the reference values
SBEM with subdomains
Classic SBEM Closed form
Number of interfacial nodes
wð10; 2Þ
Error % w respect to the classic SBEM
rx (0,1)
Error % rx respect to the classic SBEM
15 29 43 57
157:603 165:543 166:776 167:070
5.79 1.04 0.31 0.13
4:300 4:125 4:311 4:271
2.21 1.95 2.47 1.52
wð10; 2Þ ¼ 167:291 wð10; 2Þ ¼ 166:400
rx ð0; 1Þ ¼ 4:207 rx ð0; 1Þ ¼ 5:00
3362
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 7. Cantilever beam subdivided into six MEs: with (a) irregular, (b) more regular boundary shape.
Table 2 Comparison of the displacement and stress values obtained by substructuring SBE analysis considering more regular geometries of the macro-elements
SBEM with subdomains
Regularity
wð10; 2Þ
Error % w respect to the classic SBEM
rx (0, 1)
Error % rx respect to the classic SBEM
Lower
Better
157.603 159.850 162.051 163.834
5.79 4.45 3.13 2.07
4:300 4:293 4.175 4:068
2.21 2.04 0.76 3.30
Classic SBEM
wð10; 2Þ ¼ 167:291
rx ð0; 1Þ ¼ 4:207
5.2. Second example: frame The frame of Fig. 8a (E ¼ 1, m ¼ 0:3 and s ¼ 1) is subjected to a linearly variable force distribution having a unit maximum value. 5.2.1. Convergence check Let us subdivide the structure into two macro-elements having regular geometry and perform a dense boundary discretization of the two domains. In these conditions the displacement uðT Þ at the higher corner and the vertical stress ry ðSÞ at the left extreme of the constrained element RS are computed. These quantities are chosen as reference values.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3363
Fig. 8. Frame subdivided into: (a) two regular MEs, (b) 11 irregular MEs.
Table 3 Comparison of the displacement and stress values obtained by substructuring SBE analysis considering the substructurings of Fig. 8(a) and (b) the first being chosen as reference system
SBEM with subdomains
Number of interfacial nodes
uðT Þ
Error % u respect to the classic SBEM
ry ðSÞ of BEs RS
Error % ry respect to the classic SBEM
22 33 44 55 66
284.135 619.923 765.700 821.344 843.015
66.90 27.79 10.81 4.33 1.80
13:15 20:94 25:23 27:24 28:23
54.50 27.54 12.70 5.74 2.32
Classic SBEM with 2 MEs
uðT Þ ¼ 858:502
ry ðSÞ ¼ 28:90
The previous subdivision is necessary because the frame is a structural system having a boundary enclosed inside another boundary. This is the case of not simply connected bodies in which the symmetric formulation is unable to produce a unique solution due to the singularity of the matrix of the solving system. The solution could become single by giving rise to nonsymmetric matrices. In order to prove the convergence check the frame is subdivided into 11 macro-elements as in Fig. 8(b) and the boundary of each side of the subdomains at first is discretized by sparse way and successively by denser and denser way through a subdivision of the sides of the boundary elements by 2, 3 and 4. Some results of this process related in Table 3 show that, when the discretization of the macro-element boundary is denser and therefore the interfacial displacement unknowns increase, the displacements converge to the displacement values obtained by using a very dense discretization made up on the boundary of two macro-elements.
3364
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 9. Symmetric frame subjected to symmetric mechanical loads.
5.3. Other examples dealt with by displacement method inside the BE formulation • Geometrical (l ¼ 105, h ¼ 40, s ¼ 1), physical (ME 1: E ¼ 100,000, m ¼ 0:3; ME 2: E ¼ 100,000, m ¼ 0:3; ME 3: E ¼ 100,000, m ¼ 0:3) and mechanical (p ¼ 10; q ¼ 10) data. • Boundary elements: ME 1 ¼ 52; ME 2 ¼ 12; ME 3 ¼ 12. • Total interfacial nodes: n ¼ 4. • Displacement uðT Þ ¼ 0:129. • Displacement wðT Þ ¼ 1:222. • Normal stress ry ðSÞ ¼ þ436:527 (Fig. 9). • Geometrical (l ¼ 105, h ¼ 40, s ¼ 1, w1 ¼ 0:25, w2 ¼ 0:5) and physical (ME 1: E ¼ 10,000, m ¼ 0:3; ME 2: E ¼ 10,000, m ¼ 0:3; ME 3: E ¼ 10,000, m ¼ 0:3) data. • Boundary elements: ME 1 ¼ 52; ME 2 ¼ 12; ME 3 ¼ 12. • Total interfacial nodes: n ¼ 4. • Displacement uðT Þ ¼ 0:667. • Displacement wðTÞ ¼ 0:424. • Normal stress ry ðSÞ ¼ þ49:175 (Fig. 10). • Geometrical (l ¼ 4, h ¼ 4, s ¼ 1), physical (ME 1: E ¼ 100,000, m ¼ 0:2; ME 2: E ¼ 100,000, m ¼ 0:2; ME 3: E ¼ 200,000, m ¼ 0:3; ME 4: E ¼ 200,000, m ¼ 0:3) and mechanical (p ¼ 1) data. • Boundary elements: ME 1 ¼ 20, ME 2 ¼ 20, ME 3 ¼ 24, ME 4 ¼ 40. • Total interfacial nodes: n ¼ 25. • Displacement uðT Þ ¼ 4E 6. • Displacement wðT Þ ¼ 24E 6. • Displacement wðA1Þ ¼ 1E 6. • Displacement wðA2Þ ¼ 14E 6. • Normal stress ry ðSÞ ¼ 1:03 (Fig. 11). • Geometrical (l ¼ 8:25, h ¼ 7:5, s ¼ 1, w1 ¼ 1, w2 ¼ 0:5), physical (ME 1: E ¼ 8,750, m ¼ 0:3; ME 2: E ¼ 6,250, m ¼ 0:22; ME 3: E ¼ 7,500, m ¼ 0:25) and mechanical (q1 ¼ 0:5, q2 ¼ 1) data. • Boundary elements: ME 1 ¼ 54, ME 2 ¼ 24, ME 3 ¼ 43.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 10. Symmetric frame subjected to constraint subsiding of the right pillar.
Fig. 11. Symmetric plate with crack made up by two types of material subjected to symmetric mechanical loads.
• • • •
Total interfacial nodes: n ¼ 18. Displacement uðT Þ ¼ 2:040. Displacement wðT Þ ¼ 1:008. Normal stress ry ðSÞ ¼ 1590:495 (Fig. 12).
3365
3366
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 12. Frame subjected to mechanical and kinematical actions on the free and constraint boundaries, respectively.
6. Conclusions In the ambit of the boundary element theory applied to in-plane bi-dimensional systems, the symmetric formulation is taking on a bigger and bigger role in the analysis of continuum bodies. Inside this formulation a substructuring process permits us to obtain solutions for bodies made up by zone-wise different physical and/or geometrical characteristics. This goal is reached with the help of some researches permitting us to give stronger visibility to the symmetric formulation in the analysis processes. The approach proposed in the present paper employs macro-elements having large dimensions inside which the equilibrium and the compatibility are satisfied point by point by using the fundamental solutions. For each macro-element a relationship is written on the interfaces with other macro-elements connecting the weighed tractions at the boundary sides to the node displacements of the same boundary and to the external actions, after eliminating in the analysis process the unknowns defined in the remaining portions of the constrained and free boundaries. The regularity equations written on the interface sides among macro-elements give rise to a reduced number of unknowns constituted, in this case, by interface node displacements, thus making it possible to improve the response through a dense discretization. Through the use of a computing programme developed by the present authors the solution changes very little when the substructuring geometry changes, the number of boundary sides being equal, and in addition the solution converges rapidly when the interface boundary discretization is denser.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3367
This approach shows characteristics which are not far from the continuum theory and the presence of a computing programme makes it utilizable in order to study continuum problems as interface problems (crack, contact), composite material structures, coupling of SBEM-FEM, and, in any case, all the problems defined in large spaces for which a reduction of problem variables is absolutely necessary.
References [1] M. Bonnet, G. Maier, C. Polizzotto, Symmetric Galerkin boundary element method, Appl. Mech. Rev. 51 (1998) 669–704. [2] G. Maier, M. Diligenti, A. Carini, A variational approach to boundary element elastodynamic analysis and extension to multidomain problems, Comp. Meth. Appl. Mech. Engrg. 92 (1991) 192–213. [3] O. Steinbach, W.L. Wendland, Domain decomposition and preconditioning techniques in boundary element methods, in: W.L. Wendland (Ed.), Boundary Element Topics, Proceeding of the Final Conference of the Priority Research Programme Boundary Element Methods 1989–1995 of the German Research Foundation, Stuttgart, October 2–4, 1995, Springer, Berlin, 1997, pp. 471– 490. [4] L.J. Gray, G.H. Paulino, Symmetric Galerkin boundary integral formulation for interface and multi-zone problems, Int. J. Num. Meth. Engrg. 40 (1997) 3085–3101. [5] J.B. Layton, S. Ganguly, C. Balakrishna, J.H. Kane, A symmetric Galerkin multi-zone boundary element formulation, Int. J. Num. Meth. Engrg. 40 (1997) 2913–2931. [6] S. Ganguly, J.B. Layton, C. Balakrishna, J.H. Kane, A fully symmetric multi-zone Galerkin boundary element method, Int. J. Num. Meth. Engrg. 44 (1999) 991–1099. [7] T. Panzeca, M. Salerno, Macro-elements in the mixed boundary value problems, Comp. Mech. 26 (2000) 437–446. [8] T. Panzeca, M. Salerno, S. Terravecchia, Domain decomposition in the symmetric boundary element analysis, IABEM 2000 Simposium, Brescia, July 4–7, 2000. Special issue of the Comp. Mech. 25 (2002), in press.
Symmetric boundary element method versus finite element method T. Panzeca a
a,*
, F. Cucco 1, S. Terravecchia
b
Dipartimento di Ingegneria Strutturale e Geotectnica, Universit a degli Studi di Palermo, Viale delle Scienze, 90128 Palermo, Italy b Dipartimento di Strutture, Universit a della Calabria, Arcavacata di Rende, Cosenza I-87030, Italy Received 6 August 2001; received in revised form 14 January 2002
Abstract The paper examines the effectiveness of the symmetric boundary element formulation when the continuum body is subdivided into large elements called macro-elements. The approach proposed combines a strong reduction of variables with an elastic solution close to the real response. Indeed, if the displacement method is used, this approach permits one to determine for every macro-element a relationship connecting the weighted traction vector defined on the sides of the interface boundary with the node displacement vector of the same boundary and with the external action vector. Such a strategy is very similar to that followed through the finite element method, but with the advantages of having the punctual satisfaction of the equilibrium and compatibility by using the fundamental solution, and of employing large elements whose discretization is performed on the boundary only. Some examples are shown using a general computing programme developed by the present authors. 2002 Elsevier Science B.V. All rights reserved. Keywords: Symmetric boundary element method; Boundary element; Substructuring
1. Introduction In all fields of applied mechanics and engineering there is a need to work out an analysis having the capacity to combine the employment of a reduced number of variables with solutions which are closer and closer to the real response. The finite element method (FEM), using elements generically defined as having high performance, has this aim and formulations like assumed stress (Pian, Casciaro), assumed strain (Wilson, Bergan, Hughes, Felippa) and enhanced assumed strain (Simo) may be placed in this ambit. The symmetric boundary element method (SBEM) proves to be preferred for its characteristics in order to reach the two objectives (few variables, closer response) both because it performs using large elements having boundary variables only, and because the use of the fundamental solutions guarantees the punctual satisfaction of the field equilibrium and the compatibility equations. On the boundary only the statical and
*
Corresponding author. E-mail address: [email protected] (T. Panzeca). 1 External collaborator.
0045-7825/02/$ - see front matter 2002 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 2 ) 0 0 2 3 9 - 6
3348
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
kinematical conditions are written in the weighted form, and as a consequence the goodness of the solution is the result of the discretization performed on the boundary and of the type of variable modelling. Some computational problems arising in some coefficients of the solving system because of the presence of double integrals having strong singularity or hypersingularity in the definition domain have been solved by using analytical or numerical techniques. As a consequence the SBEM has had a large development in several fields of mechanics. The pertinent references are found in paper [1]. Some difficulties have arisen in the analysis of a continuum body having zone-wise different physical and geometrical characteristics subjected to external actions, but recently have been overcome by suitable strategies used in an appropriate substructuring. The latter idea was introduced by Maier et al. [2] in the ambit of variational elastic analysis in dynamics and resumed by Steinbach and Wendland [3]. Subsequently Gray and Paulino [4] used this expedient for potential problems including an example regarding heat transfer. Layton et al. [5] at first proposed a formulation which divides the body into macro-zones, each of which is governed by boundary quantities, and which after a condensation process gives rise to a system having only few asymmetric blocks. In a subsequent paper the same authors [6] presented a truly symmetric method for a two-zone elastic problem characterized by interfacial and noninterfacial unknowns. In [7,8] Panzeca et al. introduced a procedure which defines for each macro-element, obtained by appropriate substructuring, a matricial equation connecting the interfacial boundary quantities only. A subsequent assembly by means of weighted regularity conditions written on the interfaces among macroelements gives rise to a reduced equation system obtained as the solution of a problem with displacement unknowns, force unknowns or both variables. Inside the displacement approach the strategy presented in this paper proves to obtain a basic working matrix for every macro-element of domain X thought of as embedded in the infinite domain X1 having the same characteristic as X. The use of this matrix permits us to generate for every macro-element both the algebraic operators depending on the boundary type and the known vectors depending on the external actions on the element. Indeed, once the boundary type is established, i.e. constrained C1 , free C2 , interfacial C0 among macro-elements, some possible stiffness matrices can be defined having as a boundary • • • •
C0 , where the boundary is totally embedded in contiguous macro-elements; C0 [ C1 , where a portion of the boundary is constrained on the ground; C0 [ C2 , where a portion of the boundary is free; C0 [ C1 [ C2 , where all types of boundary coexist.
These stiffness matrices connect the weighted tractions defined on the sides of the interface boundary C0 to the node displacements of the same boundary. In addition, the same basic matrix permits us to derive for every macro-element the load vector including the layer forces acting on the free boundary C2 and the displacements imposed on the constraint C1 . The load vector including the body forces and the linear and volume kinematic discontinuities are obtained in the classical way. This load vector is defined on the sides of the boundary C0 where the stiffness matrix coefficients are computed. This strategy permits us to define the macro-elements of SBEM in a similar way to the finite elements of FEM, but has very high performance. As a consequence, the usual assembly rules involve a very small number of variables, i.e. the displacements of the macro-element interface boundary, and give this approach strong effectiveness because of the reduced dimension of the analysis problem. A general computing code was worked out to show the benefits of the approach proposed and the results of the examples dealt with are discussed below.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3349
2. Linear elastic analysis by substructuring Let us consider a homogeneous two-dimensional solid body B with domain X and boundary S referred to a Cartesian coordinate system ðx; yÞ. This body is subjected to imposed displacements u on the boundary portion Su restraining rigid motion, to known forces f on the free complementary boundary portion Sf ¼ S Su and to body forces b and imposed strains # in X. In this body the physical characteristics (Young’s modulus E and Poisson coefficient m) and the geometry (thickness s) are zone–wise variables. In the infinitesimal displacement hypothesis, the elastic response to the known external actions can be obtained in terms of displacements u on Sf , of forces f on Su , of displacements u, stresses r and strain in X by utilizing the SBEM. To this aim an appropriate substructuring must be employed considering the different zone–wise characteristics in the body. In Fig. 1(a) the body is subjected to external actions, and in Fig. 1(b) a domain subdivision into macro-elements is performed. For each of them an appropriate relation between mechanical and kinematical quantities defined on the boundary may be written, whose form depends on the type of approach to be employed. Indeed, in [8] the authors have shown that a substructuring analysis can be made by using approaches based on the displacements, the forces, or the mixed-value quantities. In the present paper the displacement approach will be utilized, but the same results can easily be obtained using alternative approaches. The substructuring in Fig. 1(b) involves the characterization of macro-elements whose boundaries can be constrained on C1 and free on C2 or can be interfacial on C0 among neighbouring macro-elements. In this way the macro-element A is characterized by a boundary C0 , whereas the macro-elements B, C and D, show the mixed boundaries C1 [ C0 , C2 [ C0 and C1 [ C2 [ C0 , respectively. The specification of different types of macro-elements involves all the possible cases to be characterized in substructuring a body. Let us consider in Fig. 2 the ith macro-element of domain X embedded in the infinite domain X1 having the same characteristics ðE; m; sÞ as the macro-element considered. The domain boundary may be thought of as a boundary C of X or as a boundary Cþ of X1 n X. At this stage no distinction is made for the constrained, free or interfacial boundary type. This macro-element is subjected to the layered forces f and to the double-layered displacement discontinuities v ¼ u, both on the boundary C, and to the body forces b, to thermal-like strains # and cracks d inside the domain, the symbol bar specifying that the corresponding quantities are known.
Fig. 1. (a) Two-dimensional body subjected to external actions, (b) appropriate substructuring of the domain.
3350
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 2. Macro-element embedded in the infinite domain.
According to the symbols used in [1], the Somigliana Identities (SI) can be written in the following way: u ¼ u½f þ u½u þ ^ u½ b; #; d ; ð1aÞ t ¼ t½f þ t½u þ ^t½ b; #; d
ð1bÞ
which represent the response to all the actions in X1 . In the latter the following positions are set Z Z u½f ¼ Guu f; t½f ¼ Gtu f; C
u½u ¼
ð2a; bÞ
C
Z
Gut ðuÞ;
t½u ¼
Z
C
^ u½ b; #; d ¼
Z
Guu bþ
Z
X
^t½ b; #; d ¼
Gtt ðuÞ;
ð2c; dÞ
C
Z
Gur # þ
Z
X
Gtu bþ
X
Z
Gtr # þ X
Gut d;
ð2eÞ
X
Z
Gtt d:
ð2fÞ
X
If the response is computed on the boundary Cþ in Fig. 2, one obtains u½ b; #; d ; uþ ¼ u½f þ uPV ½u 1u þ ^ 2
b; #; d
tþ ¼ tPV ½f þ t½u 12f þ ^t½
ð3aÞ ð3bÞ
where the apex PV appended to a symbol means that the corresponding quantities are computed taking the related singular integral as the Cauchy principal value. Imposing the punctual condition uþ ¼ 0, tþ ¼ 0 on the boundary Cþ permits one to obtain the following mixed boundary integral equations of the macro-element: u½ b; #; d ¼ 0; ð4aÞ u½f þ uPV ½u 1u þ ^ 2
b; #; d ¼ 0: tPV ½f þ t½u 12f þ ^t½ f, ðuÞ being the unknown function vectors defined on all the boundary of the body.
ð4bÞ
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3351
Fig. 3. Boundary forces and displacements of the nodes.
Let us execute the boundary discretization of the macro-element into the boundary elements and introduce the shape functions Wf , Wu to model the layered and double-layered boundary quantities f ¼ Wf F;
u ¼ Wu U
ð5a; bÞ
where F and U are vectors of unknown nodal values, forces and displacements, described in Fig. 3. In every boundary element the force vector F is defined through four components modelling the forces along the boundary element, whereas the displacement vector U is made up of two components per node modelling the displacements of the boundary elements contiguous to the node considered. We now introduce the modelling assumed by Eqs. (5a,b) in Eqs. (4a) and (4b) and perform the weighting of the system coefficients using the same shape functions in accordance with the Galerkin strategy. The following equations, written in weighted form, are obtained Z Z Z Z Z Z Z 1 þ Wf u Wf Guu Wf F þ Wf Gut Wu ðUÞ þ Wf Wu ðUÞ þ Wf ^u ¼ 0; ð6aÞ 2 C C C C C C C Z Z Z Z Z Z Z 1 Wu tþ Wu Gtu Wf F Wu Wf F þ Wu Gtt Wu ðUÞ þ Wu^t ¼ 0: ð6bÞ 2 C C C C C C C The latter equations written in compact form become Wþ Auu F þ Auf ðUÞ þ Cuf ðUÞ þ c W X ¼ 0;
ð7aÞ
bX ¼ 0 Pþ Afu F Cfu ðFÞ þ Aff ðUÞ þ P
ð7bÞ
with obvious symbol meaning, or also bX ¼ 0 ðA þ CÞX þ L where one sets Auu A¼ Afu
Auf ; Aff
ð8Þ
0 C¼ Cfu
Cuf ; 0
F X¼ ; U
b¼ L
c WX : bX P
ð9a–dÞ
The matrix A is symmetric. Indeed, the following positions are valid: Auu ¼ ATuu ;
Aff ¼ ATff ;
Afu ¼ ATuf
ð10a–cÞ
in that the matrix coefficients are obtained by using the fundamental solutions Ghk with h; k ¼ u; f , employing the symmetry properties [1], completed using the shape functions in accordance with the Galerkin approach. Otherwise the matrix C is hemisymmetric ðCfu ¼ CTuf Þ, Cuf being a block diagonal matrix. The sum matrix ðA þ CÞ is nonsymmetric and singular. The singularity arises because in the column block related to the vector ðUÞ, three columns are a linear combination of the remaining columns, the macro-element being subjected to a rigid motion in its plane. This circumstance takes on particular
3352
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
significance because it can be employed as a check on the column-block coefficients. Indeed, when we impose in the infinite domain a displacement discontinuity vector ðUÞ on C such as to cause a rigid motion of the domain X, we obtain Auf þ Cuf 0 ðUÞ ¼ : ð11Þ Aff 0 If this equality is satisfied, then the matrix shows coefficients computed correctly. At this stage other types of check cannot be introduced. The matrices ðA þ CÞ, written for all the macro-elements arising from the substructuring employed, take on a fundamental role in the analysis of the entire body. They are obtained by assuming the boundaries as free, but subjected to layered forces F and double-layered displacement discontinuities ðUÞ, without specifying whether the quantities are known or not. Starting from Eq. (8) it is possible to obtain an equation connecting the boundary mechanical and kinematical quantities, known and unknown, for every macro-element. 3. Symmetric operator specification In the ambit of the displacement method via SBEM, the analysis is dealt with by substructuring the body after a preliminary study made on macro-elements like those marked by capital Latin letters A, B, C and D in Fig. 1(b). For each macro-element a stiffness matrix and a load vector will be defined and computed. 3.1. Macro-elements like A This case regards a macro-element totally delimited by a boundary like C0 , i.e. a boundary interfacing with the neighbouring macro-elements. Eq. (8) is reformulated introducing the pedix 0 in all the components of the system c Au0u0 Au0f 0 þ Cu0f 0 F0 0 Wþ W X0 0 þ b ¼ ð12Þ A C A U 0 Pþ f 0u0 f 0u0 f 0f 0 P 0 X0 0 which may be rewritten in the following way: c Wþ 0 Au0u0 F0 þ ðAu0f 0 þ Cu0f 0 ÞðU0 Þ þ W X0 ¼ 0;
ð13aÞ
b Pþ 0 ðAf 0u0 Cf 0u0 ÞF0 þ Af 0f 0 ðU0 Þ þ P X0 ¼ 0:
ð13bÞ
The latter equation may be substituted by the following: b X0 P0 ¼ ðAf 0u0 þ Cf 0u0 ÞF0 þ Af 0f 0 ðU0 Þ þ P
ð14Þ
which defines the weighted traction vector on the actual boundary of the macro-element. It is to be noted that this change involves a variation in sign for the free term matrix Cf 0u0 which, as a consequence, becomes positive. The system written in compact form becomes þ c X0 ðAu0f 0 þ Cu0f 0 Þ Au0u0 F0 on Cþ W0 0 0 ¼ þ W ð15Þ b X0 ðAf 0u0 þ Cf 0u0 Þ Af 0f 0 U0 on C0 P0 P where the matrix is symmetric. We obtain F0 by using Eq. (13a), i.e. h i c F0 ¼ A1 u0u0 ðAu0f 0 þ Cu0f 0 Þð U0 Þ þ W X0
ð16Þ
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3353
and substitute it in Eq. (14), thus obtaining the following equation: b X0 P0 ¼ D00 U0 þ L
ð17Þ
where D00 ¼ ðAf 0u0 þ Cf 0u0 ÞA1 u0u0 ðAu0f 0 þ Cu0f 0 Þ Af 0f 0 ;
ð18aÞ
b X0 ðAf 0u0 þ Cf 0u0 ÞA1 c b X0 ¼ P L u0u0 W X0 :
ð18bÞ
Eq. (17) is representative of any macro-element like A and connects the weighted tractions acting on the sides of the discretized interface boundary C0 to the node displacements of the same boundary and to the domain external actions. The matrix D00 is singular because the macro-element can be subjected to a rigid motion in its plane. If the components of the vector U0 are chosen in such a way as to define a rigid motion without domain loads, the matricial product in Eq. (17) must become zero, i.e. D00 U0 ¼ 0
ð19Þ
and this check is additional the one in Eq. (11). Indeed, the first check is used to verify the column block related to the displacement discontinuity, and the second check can be utilized to verify the matrix D00 built using the block matrices characterizing the free body embedded in X1 . 3.2. Macro-elements like B The macro-element is delimited by a constrained boundary C1 and by an interface boundary C0 among neighbouring macro-elements. Eq. (8) is reformulated after rearrangement of the row blocks and the column blocks in relation to the boundary type. In this case one has 2 32 3 2 3 2 3 c Au1u1 Au1u0 Hu1f 1 Au1f 0 F1 0 Wþ W 1X 1 6 Au0u1 Au0u0 Au0f 1 Hu0f 0 76 F0 7 6 c 7 607 Wþ W 0 6 76 7 6 0X 7 ¼ 6 7 ð20aÞ 4 Ef 1u1 Af 1u0 Af 1f 1 Af 1f 0 54 U1 5 þ 4 P b 1X 5 4 0 5 Pþ 1 þ b 0X 0 U0 P0 Af 0u1 Ef 0u0 Af 0f 1 Af 0f 0 P where on sets Ef 1u1 ¼ Af 1u1 Cf 1u1 ;
Hu1f 1 ¼ Au1f 1 þ Cu1f 1 ;
ð20b; cÞ
Ef 0u0 ¼ Af 0u0 Cf 0u0 ;
Hu0f 0 ¼ Au0f 0 þ Cu0f 0
ð20d; eÞ
and where among the components of the vector which defines the boundary quantities the displacement discontinuity subvector ðU1 Þ is assigned on C1 . The latter equations may be rewritten in explicit form as follows: c Wþ 1 Au1u1 F1 þ Au1u0 F0 þ Hu1f 1 ðU1 Þ þ Au1f 0 ðU0 Þ þ W 1X ¼ 0;
ð21aÞ
c Wþ 0 Au0u1 F1 þ Au0u0 F0 þ Au0f 1 ðU1 Þ þ Hu0f 0 ðU0 Þ þ W 0X ¼ 0;
ð21bÞ
b Pþ 1 Ef 1u1 F1 þ Af 1u0 F0 þ Af 1f 1 ðU1 Þ þ Af 1f 0 ðU0 Þ þ P 1X ¼ 0;
ð21cÞ
b Pþ 0 Af 0u1 F1 þ Ef 0u0 F0 þ Af 0f 1 ðU1 Þ þ Af 0f 0 ðU0 Þ þ P 0X ¼ 0:
ð21dÞ
3354
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Eq. (21d) may be substituted by the following equation: b 0X P0 ¼ Af 0u1 F1 þ ðAf 0u0 þ Cf 0u0 ÞF0 þ Af 0f 1 ðU1 Þ þ Af 0f 0 ðU0 Þ þ P
ð22Þ
which defines the weighted traction vector on the boundary C0 of the macro-element. It is to be noted that here too this change involves a variation in sign of the free term matrix Cf 0u0 which, as a consequence, becomes positive. A block rearrangement of the matrix and of the vectors of Eq. (20a) is performed, taking into account of Eq. (22) and of the known displacement discontinuity vector (U1 ), so the system (20a) may be reformulated as follows: ð23Þ with the addition of the following equation: b Pþ 1 0 ¼ Ef 1u1 F1 þ Af 1u0 F0 þ Af 1f 0 ðU0 Þ þ Af 1f 1 ðU1 Þ þ P 1X
ð24Þ
and where Hf 0u0 ¼ HTu0f 0 . It is to be noted that the matrix in Eq. (23) is symmetric. Let us rewrite Eq. (23) in compact form, i.e. 0 ¼ BX þ B0 ðU0 Þ þ c Wu þ c WX;
ð25aÞ
bu þ P bX P0 ¼ BT0 X þ B00 ðU0 Þ þ P ð25bÞ T T T where the vector X ¼ F1 F0 has been introduced and where the remaining terms take on an obvious meaning. The vector X is obtained by using Eq. (25a) X ¼ B1 ½B0 ðU0 Þ þ c Wu þ c WX
ð26Þ
and is substituted in Eq. (25b) giving the following equation bu þ L bX P0 ¼ D00 U0 þ L ð01Þ
ð27Þ
where ð01Þ
D00 ¼ BT0 B1 B0 B00 ;
ð28aÞ
bu ¼ P b u BT B1 c L Wu; 0
ð28bÞ
b X BT B1 c bX ¼ P WX: L 0
ð28cÞ
In this case too, Eq. (27) connects the weighted traction vector defined on the sides of the discretized interface boundary C0 to the node displacement vector of the same boundary, to the domain external actions and to the load vector (U1 ) related to the constraint, the two latter vectors also being computed on the boundary elements of C0 . ð01Þ In this case the matrix D00 is nonsingular: indeed, a portion of the boundary of the macro-element is constrained on the ground. If the components of the vectors U0 and U1 are chosen in such a way as to define a rigid motion, without domain external actions, we have
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367 ð01Þ bu ¼ 0 D00 U0 þ L
3355
ð29Þ
b u in Eqs. (28a)– and L employed as second check of the coefficients of the submatrices which define (28c). Eq. (24) written on Cþ 0 is not employed to obtain the relation connecting the weighted tractions to the displacements of the boundary C0 . It may be of use as a check of the approximation of the numerical solution: the weighted traction value would have to be zero on each boundary element of Cþ 1. ð01Þ D00
3.3. Macro-elements like C The macro-element is delimited by a free boundary C2 and by an interface boundary C0 with neighbouring macro-elements. The previous Eq. (8) is reformulated after rearrangement the row blocks and the column blocks according to the boundary type. One has 2 32 3 2 3 2 3 c Au2u2 Au2u0 Hu2f 2 Au2f 0 0 Wþ F2 W 2X 2 7 6 7 6 7 6 7 þ6 W0 6 Au0u2 Au0u0 Au0f 2 Hu0f 0 76 F0 7 6 c W 7 607 ð30aÞ 76 7 þ 6 0X 7 ¼ 6 7 þ 6 b 2X 5 4 0 5 P2 4 Ef 2u2 Af 2u0 Af 2f 2 Af 2f 0 54 U2 5 4 P b 0X Af 0u2 Ef 0u0 Af 0f 2 Af 0f 0 0 Pþ U0 P 0 where, in addition to the position (20d,e), one sets Ef 2u2 ¼ Af 2u2 Cf 2u2 ;
Hu2f 2 ¼ Au2f 2 þ Cu2f 2
ð30b; cÞ
and where among the components of the vector which defines the boundary quantities the force vector F2 is assigned on the C2 . The latter equations may be rewritten in explicit form as follows: c Wþ 2 Au2u2 F2 þ Au2u0 F0 þ Hu2f 2 ðU2 Þ þ Au2f 0 ðU0 Þ þ W 2X ¼ 0;
ð31aÞ
c Wþ 0 Au0u2 F2 þ Au0u0 F0 þ Au0f 2 ðU2 Þ þ Hu0f 0 ðU0 Þ þ W 0X ¼ 0;
ð31bÞ
b Pþ 2 Ef 2u2 F2 þ Af 2u0 F0 þ Af 2f 2 ðU2 Þ þ Af 2f 0 ðU0 Þ þ P 2X ¼ 0;
ð31cÞ
b Pþ 0 Af 0u2 F2 þ Ef 0u0 F0 þ Af 0f 2 ðU2 Þ þ Af 0f 0 ðU0 Þ þ P 0X ¼ 0:
ð31dÞ
The use of the procedure adopted for the macro-element like C, and therefore the writing of Eq. (31d) on C0 instead of Cþ 0 and the arrangement of the block matrices and of the block vectors in Eq. (30a) permit us to reformulate the problem as follows: ð32Þ with the following equation: c Wþ 2 Au2u0 F0 þ Hu2f 2 ðU2 Þ þ Au2f 0 ðU0 Þ þ Au2u2 ðF2 Þ þ W 2X ¼ 0:
ð33Þ
The matrix in Eq. (32) is symmetric. The block equation may be written in the following compact form: 0 ¼ BX þ B0 ðU0 Þ þ c Wf þ c WX;
ð34aÞ
bf þ P bX P0 ¼ BT0 X þ B00 ðU0 Þ þ P
ð34bÞ
3356
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
where the vector XT ¼ FT0 UT2 has been introduced and where the remaining terms take on an obvious meaning. Obtaining the vector X by using Eq. (34a) and substituting it in Eq. (34b) the following equation can be attained: ð02Þ bf þ L bX P0 ¼ D00 U0 þ L
ð35Þ
where ð02Þ
D00 ¼ BT0 B1 B0 B00 ;
ð36aÞ
b f BT B1 c bf ¼ P Wf ; L 0
ð36bÞ
b X BT B1 c bX ¼ P WX: L 0
ð36cÞ
In this case too the final expression (35) connects the weighted traction vector on the discretized boundary elements of C0 to the node displacement vector, whose components are defined on the same boundary, to the domain external actions and to the force vector F2 related to the free boundary, both latter vectors also being computed on C0 . ð02Þ The matrix D00 is symmetric and singular: indeed, no portion of the boundary of the macro-element is constrained on the ground. As a consequence, when the vector U0 is chosen in such a way as to define a rigid motion, we must have ð02Þ
D00 U0 ¼ 0:
ð37Þ
3.4. Macro-elements like D All the boundary types are present, macro-elements. The block form of Eq. 2 Au1u1 Au1u2 Au1u0 Hu1f 1 Wþ 1 þ6 W2 6 Au2u1 Au2u2 Au2u0 Au2f 1 6 6 Au0u2 Au0u0 Au0f 1 Wþ 0 6 Au0u1 þ 6 P1 6 Ef 1u1 Af 1u2 Af 1u0 Af 1f 1 6 4 Af 2u1 Ef 2u2 Af 2u0 Af 2f 1 Pþ 2 þ P0 Af 0u1 Af 0u2 Ef 0u0 Af 0f 1
i.e. C1 constrained, C2 free, C0 interfacial among neighbouring (8) becomes the following: 32 3 2c 3 2 3 Au1f 2 Au1f 0 0 F1 W 1X 7 6 7 6 7 7 6 c 7 6 Hu2f 2 Au2f 0 76 F2 7 W 2X 7 6 0 7 76 7 6 7 607 c 7 Au0f 2 Hu0f 0 76 F0 7 6 W 0X 7 6 7 76 7þ6 ð38Þ 7¼6 6 6 7 7 6 b 1X 7 6 0 7 Af 1f 2 Af 1f 0 76 U1 7 6 P 7 7 6 6 7 76 7 b 2X 7 Af 2f 2 Af 2f 0 54 U2 5 6 5 405 4P b 0X Af 0f 2 Af 0f 0 0 U0 P
where the positions (20b,c), (20d,e) and (30b,c) are valid. The evaluation of the vector P0 on C0 instead of þ the vector Pþ 0 on C0 and an appropriate arrangement produce the following matricial equation:
ð39Þ
with the equations c Wþ 2 0 ¼ Au2u1 F1 þ Au2u0 F0 þ Hu2f 2 ðU2 Þ þ Au2f 0 ðU0 Þ þ Au2u2 F2 þ Au2f 1 ðU1 Þ þ W 2X ;
ð40aÞ
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
b Pþ 1 0 ¼ Ef 1u1 F1 þ Af 1u0 F0 þ Af 1f 2 ðU2 Þ þ Af 1f 0 ðU0 Þ þ Af 1u2 F2 þ Af 1f 1 ðU1 Þ þ P 1X :
3357
ð40bÞ
Eq. (39) is rewritten in compact form with an obvious meaning of the symbols 0 ¼ BX þ B0 ðU0 Þ þ c Wf þ c Wu þ c WX;
ð41aÞ
bf þ P bu þ P bX P0 ¼ BT0 X þ B00 ðU0 Þ þ P
ð41bÞ
and therefore ð012Þ bb þ L bX P0 ¼ D00 U0 þ L
ð42Þ
where ð012Þ
D00
¼ BT0 B1 B0 B00 ;
ð43aÞ
bf þ L b u ¼ ðP b f BT B1 c b u BT B1 c bb ¼ L Wf Þ þ ðP W u Þ; L 0 0
ð43bÞ
b X BT B1 c bX ¼ P WX: L 0
ð43cÞ
ð012Þ D00
The matrix is symmetric and nonsingular: indeed, a portion of the boundary is constrained on the ground. A rigid motion involves the following check: b u ¼ 0: D00 U0 þ L ð012Þ
ð44Þ
4. Solving system In order to obtain the elastic response of a structure subjected to known external actions, a two-dimensional continuum body serving as an explanation is considered in Fig. 4. As a consequence of the different physical ðE; mÞ and geometrical (s) characteristics it is subdivided into three macro-zones specified by the capital Latin letters A, B and C. The macro-elements A, B and C belong to macro-elements of type D, C and C respectively, and each of them is characterized by a relationship between the weighted tractions and the node displacements, both quantities defined on C0 , as can be seen in Eqs. (17), (27), (35) and (42) in the previous section. The macro-element interfaces are called by the small Latin letters a, b and c and for simplicity the double indices 00 characterizing the interface and the apices characterizing the macro-element type are suppressed. In each stiffness matrix D employed a rearrangement of the block rows and the block columns is made in order to assemble the quantities associated with every boundary dividing the vector P in function of the connection points among macro-elements. For the element A we will have
Fig. 4. Body subdivided into three macro-zones A, B and C.
3358
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
2
32 A Daa PAa 6 A 76 A 4 Pc 54 Dca PA0 DA0a
DAac DAcc DA0c
32 A 3 2 A 3 b P Ua DAa0 a 7 6 7 6 bA 7 DAc0 54 UAc 5 þ 4 P c 5 b DA00 UA0 PA
ð45Þ
0
or in compact form b A: P A ¼ DA UA þ P
ð46Þ A
UAa
UAc
The displacement vector U collects the node displacements of the boundary a, those of the boundary c, and UA0 those of the common node with both the boundary elements a and b. An equation like (46) can be written for the macro-element B and C. In this way one obtains a global relation connecting the weighted tractions to the displacements of the nodes concerning the same interface boundary b P ¼ DU þ P
ð47Þ
where the following definitions are valid: bT ¼ P b AT PT ¼ PAT PBT PCT ; P UT ¼ UAT
UBT
UCT ;
b CT ; P
b BT P
D ¼ diag DA
DB
ð48a; bÞ
DC :
ð48c; dÞ
The system compatibility equation is obtained by introducing the displacement vector of the interface nodes and that of the nodes which are common vertices with the macro-elements, i.e. uT ¼ uTa uTb uTc uT0 : ð49Þ As a consequence U ¼ Zu
ð50Þ
where Z is the topological matrix. The regularity condition of the weighted tractions takes on the form ZT P ¼ 0:
ð51Þ
The use of Eqs. (47), (50), and (51) gives rise to the following solving system: Ku þ ^f ¼ 0
ð52Þ
where for the structure being examined
ð53aÞ
and h ^f T ¼ ðZT P b ÞT ¼ ð P bA þ P b B ÞT a a
bB þ P b C ÞT ðP b b
bA þ P b C ÞT ðP c c
i bA þ P bB þ P b C ÞT : ðP 0 0 0
ð53bÞ
Here too the rigid motion technique can be employed in order to check the correctness of the assembly performed, but it is to be noted that this technique must involve both the vector u and the equivalent force vector ^f through the displacements imposed on the constraints present.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3359
Using Eq. (52) gives the displacement vector u and in succession we obtain all the boundary quantities of each macro-element and therefore through the SI the elastic response in the domains.
5. Numerical applications The first example to be shown includes two types of check, one to prove the convergence, the other to verify the goodness of the solution when the substructuring changes its geometry. In both the checks the solution chosen for comparison is obtained by considering the body as a single macro-element and performing a dense discretization on the boundary or, where it is possible, by making reference to close form solutions. The check is made on both the displacement and the stress at meaningful points in the body. The convergence check is obtained by • subdividing the body into macro-elements utilizing a sparse discretization of the boundary, and • preserving the geometry of the macro-elements, but increasing the number of the boundary elements on the sides in the same ratio ð2; 3; 4; . . .Þ in each macro-element. The check on the goodness of the solution contemplating the change in the substructure geometry of the macro-elements is obtained by • subdividing the body into macro-elements whose shapes show several unevennesses and performing a pre-established boundary discretization, and • preserving the number of boundary elements at the sides, but changing the geometry of the macro-elements inside the domain towards more regular shapes. In the first two examples the physical and geometrical characteristics of each macro-element are assumed equal, and this is done to permit the checks, but the approach employed and the computing programme make it possible to have different values for the Young modulus E, the Poisson coefficient m and for the thickness s of each macro-element. The modelling of the displacements and of the tractions on the boundary elements is assumed linear. 5.1. First example: cantilever beam Let us consider a beam (l h ¼ 10 4) constrained at the first extreme against rigid motion and subjected to a unit shear force distributed by a parabolic law at the tip (Fig. 5). Let us assume E ¼ 1, m ¼ 0:3, and the thickness s ¼ 1. The vertical displacement w of the point in the corner of the tip (x ¼ 10, y ¼ 2) and the axial stress rx at the upper extreme of the constrained element 54–55 (x ¼ 0, y ¼ 1) are chosen as quantities to be checked. At first the beam is considered as a single domain whose boundary is discretized by introducing a denser discretization made up by 56 boundary elements and characterized by 47 free nodes. The values obtained at the points chosen by the classic symmetric BE analysis using the discretization mentioned are taken as comparison quantities in the computations below. 5.1.1. Convergence check Let us assume a domain subdivided by six macro-elements, as in Fig. 6a, characterized by nodes, these being the vertices of the elements. This domain, thus subdivided, is assumed unchanged in the test. The
3360
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 5. Cantilever beam subjected to a unit shear force at the tip.
Fig. 6. Cantilever beam subdivided into six MEs: (a) with a single boundary element for every side, (b) with four boundary elements for every side.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3361
constrained and the tip boundaries will have an unchanged discretization, each made up by eight elements, whereas each side of the macro-element boundaries will be subdivided into 2, 3 and 4 elements (see Fig. 6b). The analysis is performed using the displacement method through the interface node displacements, therefore the progressive subdivision of the boundary sides of each macro-element leads to the following numbers of the free nodes: ME ME ME ME ME ME
1: 2: 3: 4: 5: 6:
13–18–23–28. 4–8–12–16. 7–14–21–28. 6–12–18–24. 9–18–27–36. 13–18–23–28.
When the macro-element boundary discretization increases, the number of the interfacial nodes involved in the equation system increases in the following way: 15–29–43–57. It is to be noted that the analysis of the beam subdivided into six macro-elements using few interfacial variables permits us to obtain displacements and tractions which are very close to the values obtained using the classic symmetric approach in the beam thought of as a single element (see Table 1). The reference closed form values in Table 1 are obtained by considering as a hypothesis a tangential stress having a parabolic form in the constrained portion. The true tangential stress distribution shows a different shape and as a consequence the response of both the displacements and the stresses is different from the theoretical one. A variable response is obtained for the stress: indeed, the value is a floating one around the reference solution. 5.1.2. Check of the goodness of the solution The substructuring is made up by six macro-elements each having very irregular geometry and the boundary sides each constituting a single boundary element, with the exception of the constraint and tip boundaries where a dense discretization (eight elements) is made. A subsequent modification of the macro-elements towards more regular geometries is performed and for every change it is to be noted that the solution in terms of displacements becomes close to the values obtained by the classic SBEM without substructuring chosen as a reference. In this case too the stress value is a floating one around the reference solution (Fig. 7, Table 2).
Table 1 Comparison between the displacement and stress values obtained by substructuring SBE analysis and the reference values
SBEM with subdomains
Classic SBEM Closed form
Number of interfacial nodes
wð10; 2Þ
Error % w respect to the classic SBEM
rx (0,1)
Error % rx respect to the classic SBEM
15 29 43 57
157:603 165:543 166:776 167:070
5.79 1.04 0.31 0.13
4:300 4:125 4:311 4:271
2.21 1.95 2.47 1.52
wð10; 2Þ ¼ 167:291 wð10; 2Þ ¼ 166:400
rx ð0; 1Þ ¼ 4:207 rx ð0; 1Þ ¼ 5:00
3362
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 7. Cantilever beam subdivided into six MEs: with (a) irregular, (b) more regular boundary shape.
Table 2 Comparison of the displacement and stress values obtained by substructuring SBE analysis considering more regular geometries of the macro-elements
SBEM with subdomains
Regularity
wð10; 2Þ
Error % w respect to the classic SBEM
rx (0, 1)
Error % rx respect to the classic SBEM
Lower
Better
157.603 159.850 162.051 163.834
5.79 4.45 3.13 2.07
4:300 4:293 4.175 4:068
2.21 2.04 0.76 3.30
Classic SBEM
wð10; 2Þ ¼ 167:291
rx ð0; 1Þ ¼ 4:207
5.2. Second example: frame The frame of Fig. 8a (E ¼ 1, m ¼ 0:3 and s ¼ 1) is subjected to a linearly variable force distribution having a unit maximum value. 5.2.1. Convergence check Let us subdivide the structure into two macro-elements having regular geometry and perform a dense boundary discretization of the two domains. In these conditions the displacement uðT Þ at the higher corner and the vertical stress ry ðSÞ at the left extreme of the constrained element RS are computed. These quantities are chosen as reference values.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3363
Fig. 8. Frame subdivided into: (a) two regular MEs, (b) 11 irregular MEs.
Table 3 Comparison of the displacement and stress values obtained by substructuring SBE analysis considering the substructurings of Fig. 8(a) and (b) the first being chosen as reference system
SBEM with subdomains
Number of interfacial nodes
uðT Þ
Error % u respect to the classic SBEM
ry ðSÞ of BEs RS
Error % ry respect to the classic SBEM
22 33 44 55 66
284.135 619.923 765.700 821.344 843.015
66.90 27.79 10.81 4.33 1.80
13:15 20:94 25:23 27:24 28:23
54.50 27.54 12.70 5.74 2.32
Classic SBEM with 2 MEs
uðT Þ ¼ 858:502
ry ðSÞ ¼ 28:90
The previous subdivision is necessary because the frame is a structural system having a boundary enclosed inside another boundary. This is the case of not simply connected bodies in which the symmetric formulation is unable to produce a unique solution due to the singularity of the matrix of the solving system. The solution could become single by giving rise to nonsymmetric matrices. In order to prove the convergence check the frame is subdivided into 11 macro-elements as in Fig. 8(b) and the boundary of each side of the subdomains at first is discretized by sparse way and successively by denser and denser way through a subdivision of the sides of the boundary elements by 2, 3 and 4. Some results of this process related in Table 3 show that, when the discretization of the macro-element boundary is denser and therefore the interfacial displacement unknowns increase, the displacements converge to the displacement values obtained by using a very dense discretization made up on the boundary of two macro-elements.
3364
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 9. Symmetric frame subjected to symmetric mechanical loads.
5.3. Other examples dealt with by displacement method inside the BE formulation • Geometrical (l ¼ 105, h ¼ 40, s ¼ 1), physical (ME 1: E ¼ 100,000, m ¼ 0:3; ME 2: E ¼ 100,000, m ¼ 0:3; ME 3: E ¼ 100,000, m ¼ 0:3) and mechanical (p ¼ 10; q ¼ 10) data. • Boundary elements: ME 1 ¼ 52; ME 2 ¼ 12; ME 3 ¼ 12. • Total interfacial nodes: n ¼ 4. • Displacement uðT Þ ¼ 0:129. • Displacement wðT Þ ¼ 1:222. • Normal stress ry ðSÞ ¼ þ436:527 (Fig. 9). • Geometrical (l ¼ 105, h ¼ 40, s ¼ 1, w1 ¼ 0:25, w2 ¼ 0:5) and physical (ME 1: E ¼ 10,000, m ¼ 0:3; ME 2: E ¼ 10,000, m ¼ 0:3; ME 3: E ¼ 10,000, m ¼ 0:3) data. • Boundary elements: ME 1 ¼ 52; ME 2 ¼ 12; ME 3 ¼ 12. • Total interfacial nodes: n ¼ 4. • Displacement uðT Þ ¼ 0:667. • Displacement wðTÞ ¼ 0:424. • Normal stress ry ðSÞ ¼ þ49:175 (Fig. 10). • Geometrical (l ¼ 4, h ¼ 4, s ¼ 1), physical (ME 1: E ¼ 100,000, m ¼ 0:2; ME 2: E ¼ 100,000, m ¼ 0:2; ME 3: E ¼ 200,000, m ¼ 0:3; ME 4: E ¼ 200,000, m ¼ 0:3) and mechanical (p ¼ 1) data. • Boundary elements: ME 1 ¼ 20, ME 2 ¼ 20, ME 3 ¼ 24, ME 4 ¼ 40. • Total interfacial nodes: n ¼ 25. • Displacement uðT Þ ¼ 4E 6. • Displacement wðT Þ ¼ 24E 6. • Displacement wðA1Þ ¼ 1E 6. • Displacement wðA2Þ ¼ 14E 6. • Normal stress ry ðSÞ ¼ 1:03 (Fig. 11). • Geometrical (l ¼ 8:25, h ¼ 7:5, s ¼ 1, w1 ¼ 1, w2 ¼ 0:5), physical (ME 1: E ¼ 8,750, m ¼ 0:3; ME 2: E ¼ 6,250, m ¼ 0:22; ME 3: E ¼ 7,500, m ¼ 0:25) and mechanical (q1 ¼ 0:5, q2 ¼ 1) data. • Boundary elements: ME 1 ¼ 54, ME 2 ¼ 24, ME 3 ¼ 43.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 10. Symmetric frame subjected to constraint subsiding of the right pillar.
Fig. 11. Symmetric plate with crack made up by two types of material subjected to symmetric mechanical loads.
• • • •
Total interfacial nodes: n ¼ 18. Displacement uðT Þ ¼ 2:040. Displacement wðT Þ ¼ 1:008. Normal stress ry ðSÞ ¼ 1590:495 (Fig. 12).
3365
3366
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
Fig. 12. Frame subjected to mechanical and kinematical actions on the free and constraint boundaries, respectively.
6. Conclusions In the ambit of the boundary element theory applied to in-plane bi-dimensional systems, the symmetric formulation is taking on a bigger and bigger role in the analysis of continuum bodies. Inside this formulation a substructuring process permits us to obtain solutions for bodies made up by zone-wise different physical and/or geometrical characteristics. This goal is reached with the help of some researches permitting us to give stronger visibility to the symmetric formulation in the analysis processes. The approach proposed in the present paper employs macro-elements having large dimensions inside which the equilibrium and the compatibility are satisfied point by point by using the fundamental solutions. For each macro-element a relationship is written on the interfaces with other macro-elements connecting the weighed tractions at the boundary sides to the node displacements of the same boundary and to the external actions, after eliminating in the analysis process the unknowns defined in the remaining portions of the constrained and free boundaries. The regularity equations written on the interface sides among macro-elements give rise to a reduced number of unknowns constituted, in this case, by interface node displacements, thus making it possible to improve the response through a dense discretization. Through the use of a computing programme developed by the present authors the solution changes very little when the substructuring geometry changes, the number of boundary sides being equal, and in addition the solution converges rapidly when the interface boundary discretization is denser.
T. Panzeca et al. / Comput. Methods Appl. Mech. Engrg. 191 (2002) 3347–3367
3367
This approach shows characteristics which are not far from the continuum theory and the presence of a computing programme makes it utilizable in order to study continuum problems as interface problems (crack, contact), composite material structures, coupling of SBEM-FEM, and, in any case, all the problems defined in large spaces for which a reduction of problem variables is absolutely necessary.
References [1] M. Bonnet, G. Maier, C. Polizzotto, Symmetric Galerkin boundary element method, Appl. Mech. Rev. 51 (1998) 669–704. [2] G. Maier, M. Diligenti, A. Carini, A variational approach to boundary element elastodynamic analysis and extension to multidomain problems, Comp. Meth. Appl. Mech. Engrg. 92 (1991) 192–213. [3] O. Steinbach, W.L. Wendland, Domain decomposition and preconditioning techniques in boundary element methods, in: W.L. Wendland (Ed.), Boundary Element Topics, Proceeding of the Final Conference of the Priority Research Programme Boundary Element Methods 1989–1995 of the German Research Foundation, Stuttgart, October 2–4, 1995, Springer, Berlin, 1997, pp. 471– 490. [4] L.J. Gray, G.H. Paulino, Symmetric Galerkin boundary integral formulation for interface and multi-zone problems, Int. J. Num. Meth. Engrg. 40 (1997) 3085–3101. [5] J.B. Layton, S. Ganguly, C. Balakrishna, J.H. Kane, A symmetric Galerkin multi-zone boundary element formulation, Int. J. Num. Meth. Engrg. 40 (1997) 2913–2931. [6] S. Ganguly, J.B. Layton, C. Balakrishna, J.H. Kane, A fully symmetric multi-zone Galerkin boundary element method, Int. J. Num. Meth. Engrg. 44 (1999) 991–1099. [7] T. Panzeca, M. Salerno, Macro-elements in the mixed boundary value problems, Comp. Mech. 26 (2000) 437–446. [8] T. Panzeca, M. Salerno, S. Terravecchia, Domain decomposition in the symmetric boundary element analysis, IABEM 2000 Simposium, Brescia, July 4–7, 2000. Special issue of the Comp. Mech. 25 (2002), in press.