The scaled boundary finite-element method – a primer: derivations

Computers and Structures 78 (2000) 191±210

www.elsevier.com/locate/compstruc

The scaled boundary ®nite-element method ± a primer: derivations John P. Wolf *, Chongmin Song 1 Department of Civil Engineering, Institute of Hydraulics and Energy, Swiss Federal Institute of Technology, Lausanne, CH-1015 Lausanne, Switzerland Received 28 May 1999; accepted 20 October 1999

Abstract The scaled boundary ®nite-element method is a semi-analytical fundamental-solution-less boundary-element method based solely on ®nite elements. Using the simplest wave propagation problem and discretizing the boundary with a twonode line ®nite element, which preserves all essential features, two derivations of the scaled boundary ®nite-element equations in displacement and dynamic sti€ness are presented. In the ®rst, the scaled-boundary-transformation-based derivation, the new local coordinate system consists of the distance measured from the so-called scaling centre and the circumferential directions de®ned on the surface ®nite element. The governing partial di€erential equations are transformed to ordinary di€erential equations by applying the weighted-residual technique. The boundary conditions are conveniently formulated in the local coordinates. In the second, the mechanically based derivation, a similar ®ctitious boundary is introduced. A ®nite-element cell is constructed between the two boundaries. Standard ®nite-element assemblage and similarity lead to the scaled boundary ®nite-element equations after performing the limit of the cell width towards zero analytically. Ó 2000 Civil-Comp Ltd. and Elsevier Science Ltd. All rights reserved. Keywords: Boundary element; Dynamics; Finite element; Radiation condition; Soil-structure interaction; Wave motion

1. Introduction The most widely used computational procedures in solid mechanics are the ®nite-element method and the boundary-element method with their own speci®c features, advantages and disadvantages. For the sake of illustration, the linear analysis in elastodynamics of a bounded medium and of an unbounded (in®nite) medium is addressed (Fig. 1). In the domain V, the equations of motion apply, which are partial di€erential equations in displacements u. Body loads p are present. On the boundary S, either the displacements u on Su , or the surface tractions t on St are prescribed. For an un-

* Corresponding author. Tel.: +4121-693-2405; fax: +4121693-2264. 1 Present address: Hilti Corporate Research; FL-9494 Schaan, Liechtenstein.

bounded medium (Fig. 1b), an additional boundary condition at in®nity (``radiation condition'') must be satis®ed. As initial conditions, the displacements and velocities vanish in the domain. In the ®nite-element method, the domain is spatially discretized, as shown for a bounded medium in Fig. 2a. In each ®nite element, shape functions in the form of polynomials interpolate the displacements. Standard numerical integration of regular functions leads to the static-sti€ness and mass matrices which are then assembled enforcing compatibility and equilibrium. A great ¯exibility exists in representing the geometry and the materials. In the boundary-element method [1,2], only the boundary is discretized (Fig. 2b), leading to a reduction of the spatial dimension by one. This diminishes the e€ort of data preparation and leads to fewer unknowns. However, a fundamental solution satisfying the governing equations in the domain must be available.

0045-7949/00/$ - see front matter Ó 2000 Civil-Comp Ltd. and Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 0 ) 0 0 0 9 9 - 7

192

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

Fig. 1. Problem de®nition in elastodynamics: (a) bounded medium and (b) unbounded medium.

Fig. 2. Spatial discretization: (a) ®nite-element method and (b) boundary-element method.

This analytical solution is often very complicated exhibiting singularites. In each boundary element shape functions in the form of polynomials interpolate the displacements and surface tractions. Special numerical integration of the polynomials with the fundamental solution involving singularities yields non-symmetric coecient matrices. For an unbounded medium (Fig. 1b), the ®niteelement method cannot represent the radiation condi-

tion at in®nity exactly. The spatial discretization of ®nite elements is terminated on an arti®cial boundary, where the dynamic behaviour of the truncated domain outside of this boundary is modelled only approximately. In contrast, the boundary-element method is well suited to model an unbounded medium as the fundamental solution satis®es the radiation condition at in®nity exactly. As already mentioned, only the boundary is discretized. In a frequency-domain analysis ®ctitious eigenfrequencies occur which requires a special treatment. The ®nite-element and boundary-element methods converge, in general, to the exact solution for decreasing element size. The convergence is slow close to the point of stress singularities as occurring around the crack tip in fracture mechanics, because both methods use polynomials to interpolate the displacements. Special techniques and a large number of elements are necessary to achieve high accuracy. A novel semi-analytical procedure called the scaled boundary ®nite-element method has been developed. After discretizing the boundary with surface ®nite elements, the governing partial di€erential equations in the domain are transformed to ordinary di€erential equations which can be solved analytically. A short characterization follows which will become apparent further on. The procedure is a fundamental-solution-less boundary-element method based solely on ®nite elements. The spatial dimension is reduced by one. No singular integrals must be evaluated. General anisotropic material can be analysed without increasing the computational e€ort. For an unbounded (semi-in®nite or in®nite) medium the radiation condition at in®nity is satis®ed exactly. No ®ctitious eigenfrequencies are present when modelling an unbounded medium. The scaled boundary ®nite-element method combines the advantages of the ®nite-element and boundary-element methods. In addition, the novel procedure exhibits advantages of its own. As will be demonstrated, part of the boundary Au , At does not have to be discretized, as shown in Fig. 3a for a bounded medium and in Fig. 3b for an unbounded medium. This feature allows the crack faces in a bounded medium and the free surface in an unbounded (semi-in®nite) medium such as a half-space to be represented without discretization. These boundary conditions on Au , At are satis®ed exactly. The same advantage also applies to certain interfaces between different materials. In the standard boundary-element method, the surfaces Au , At would have to be discretized when the fundamental solution of the full space is used. As another unique advantage of the scaled boundary ®nite-element method, the semi-analytical solution in the domain permits an ecient and highly accurate calculation of stress intensity factors. The derivation and solution procedures of the fundamental equations of the scaled boundary ®nite-

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

193

sponding truncated semi-in®nite wedge. In Section 3, this illustrative example is de®ned. Two derivations of the scaled boundary ®nite-element equations are presented. In Section 4, the so-called scaled boundary transformation and the weighted-residual technique lead to the scaled boundary ®nite-element equations in displacement and in dynamic sti€ness. In Section 5, the same equations are obtained by the mechanically based derivation using similarity and ®nite-element assemblage. In the accompanying paper [9] the analytical and numerical solution procedures of the scaled boundary ®nite-element equations in displacement and dynamic sti€ness are discussed using the illustrative example.

2. Concepts of scaled boundary transformation of geometry and of similarity Fig. 3. Part of boundary with no spatial discretization: (a) bounded medium and (b) unbounded medium.

element method were published in various journals as the development proceeded. In addition, applications in di€erent ®elds, such as wave propagation, di€usion and statics, were addressed. Key papers include a mechanically based derivation for elastodynamics of the scaled boundary ®nite-element equation in dynamic sti€ness where the procedure is called the consistent in®nitesimal ®nite-element cell method re¯ecting its derivation [3], a scaled-boundary-transformation-based derivation starting from the governing di€erential equations using a Galerkin weighted-residual technique resulting in the scaled boundary ®nite-element equation in displacement [4], a numerical solution of the scaled boundary ®niteelement equation in dynamic sti€ness [5], a formulation and a numerical solution of the displacement unitimpulse response [6], and an analytical solution in the frequency domain of the scaled boundary ®nite-element equation in displacement [7]. A book addressing the consistent in®nitesimal ®nite-element cell method with extensions to di€usion and bounded media is also available [8]. It is the goal of this paper and the accompanying paper [9] to streamline the derivations and solution procedures. In addition, the solution procedures are expanded to include body loads. In Section 2 of this paper, the concepts of the scaled boundary transformation of the geometry and of similarity are discussed, which permits the domain to be described by the boundary discretization with ®nite elements and the radial coordinates. To make the paper easy to follow, the simplest wave propagation problem is discussed which still contains all essential features: the scalar wave equation in a two-dimensional wedge and the corre-

To explain the concept of the scaled boundary ®niteelement method, a three-dimensional linear elastic bounded medium with a section shown in Fig. 1a is addressed. The so-called scaling centre O is chosen in a zone from which the total boundary must be visible (Fig. 4a). This can always be achieved by subdividing the total domain in subdomains, which can be regarded as superelements. For a bounded medium the scaling centre will, in general, be located inside the domain. Only the boundary is discretized with doublycurved surface ®nite elements with any arrangement of nodes (Fig. 4a). A typical ®nite element with the surface S e (superscript e for element) is shown in Fig. 5. Two local curvilinear coordinates g, f in the circumferential directions of the surface are de®ned. This is similar to the procedure used in an isoparametric surface ®nite

Fig. 4. Modelling of bounded medium with surface ®nite elements (section) and choice of scaling centre: (a) inside domain and (b) on boundary.

194

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

Fig. 5. Scaled boundary transformation of geometry of surface ®nite element forming pyramid with volume V e .

element. The radial coordinate n points from the scaling centre O to a point on the boundary (Figs. 4a and 5). n ˆ 0 in O and n ˆ 1 on the boundary are chosen. Connecting the edge of the surface ®nite element by straight lines to the scaling centre O forms the side-face Ae , which corresponds to either g ˆ 1 or f ˆ 7. The side-face Ae and the surface S e de®ne a pyramid with the apex O and the volume V e . The geometry of this pyramid is conveniently de®ned by the scaled boundary coordinates n, g, f with 0 6 n 6 1. A transformation of the geometry from the Cartesian coordinates x^, y^, ^z to the dimensionless coordinates n, g, f is performed. The transformation is unique due to the choice of the scaling centre. This transformation is similar to the representation of a spherical surface introducing spherical coordinates where on the surface the radial coordinate is constant. The dimensionless radial coordinate n can be interpreted as a scaling factor. As already mentioned, n ˆ 1 corresponds to the surface ®nite element on the boundary. A surface corresponding to a constant n smaller than 1 is similar to the surface ®nite element on the boundary. O is denoted as the scaling centre as this point is obtained for n ˆ 0. O can also be called the similarity centre, and n represents a (dimensionless) characteristic length. This transformation de®ned by the scaling factor and the local coordinates of the surface ®nite elements is called the scaled boundary transformation. Assembling all the pyramids by connecting their side-faces which corresponds to enforcing compatibility and equilibrium, results in the total medium with volume V and the closed boundary S. No side-faces Ae passing through the scaling centre remain. When during the assemblage process, some side-faces of the pyramids are not connected or some pyramids are missing, an additional boundary A passing through the scaling centre O is created. The case of a missing pyramid is illustrated in Fig. 4b with the displacements and surface tractions prescribed on Au and At , respectively.

As the side-faces are obtained by scaling the edges of the surface ®nite elements on the boundary, the boundary Au , At is modelled without spatial discretization. In this case, the total boundary is decomposed into two parts: that part of the boundary passing through the scaling centre denoted as the side-face A and the remaining part called the interface S (Fig. 4b). The spatial dimension of the side-face A is reduced by one when viewed from O. The term interface is used only when a distinction with the word boundary is necessary. Only the interface is discretized with (doubly curved) surface ®nite elements. The scaling centre O is thus on the (total) boundary. Applications of this powerful extension exist in fracture mechanics, where the scaling centre is chosen at the crack tip. The adjacent sides of pyramids forming the crack faces are not connected. Note that the crack faces are not discretized. The same concept also applies to an unbounded medium shown in Fig. 1b. The scaling centre O is chosen outside the unbounded medium, from which the total boundary is visible (Fig. 6a) with a typical ®nite element shown in Fig. 5. The curvi-linear coordinates g, f on the surface of the boundary remain unchanged which still corresponds to n ˆ 1. A truncated semi-in®nite pyramid with its apex at the scaling centre O is de®ned by 1 6 n < 1 (Fig. 6a). Assembling all the pyramids corresponding to the surface ®nite elements yields the unbounded medium. When during the assemblage process some side-faces are not connected or some pyramids are

Fig. 6. Modelling of an unbounded medium with surface ®nite elements (section) and choice of scaling centre: (a) outside domain and (b) on extension of boundary.

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

missing, the remaining side-faces of the pyramids form an additional boundary (Au and At ) extending to in®nity (Fig. 6b). The extension of this side-face Au , At of the unbounded medium passes through the scaling centre O. Again, the side-face is not discretized. This feature allows the free surface of a half-space to be modelled without discretization.

nomenclature, the circum¯ex^is omitted when used in a subscript to indicate a direction. The stresses acting perpendicular to the plane are equal to  fsg ˆ

sx sy



  u;x^ ˆG u;y^

…2†

with the shear modulus G. Substituting Eq. (2) in Eq. (1) leads to the scalar wave equation,

3. Illustrative example The illustrative example used throughout the paper is shown for the homogeneous bounded medium and the unbounded medium in Fig. 7. The coordinates of a point in the domain are denoted as x^, y^ as x, y are reserved for the boundary. The out-of-plane (anti-plane) motion with the displacement u perpendicular to the plane is discussed. Formulating equilibrium in the Cartesian coordinate system yields sx;^x ‡ sy;^y ÿ q u‡p ˆ0

195

…1†

with the body load p per unit area acting perpendicular to the plane and the mass density q. To simplify the

u;x^x^ ‡u;y^y^ ÿ

1 p  u‡ ˆ0 c2s G

…3†

p with the shear-wave velocity, cs ˆ G=q. For the bounded medium, a wedge in the shape of a triangle is addressed (Fig. 7a). The three vertices are the origin O and Points 1 and 2 with the indicated coordinates. On the sides O±1 and O±2, the displacement u ˆ 0 and the surface traction sn ˆ fngT fsg ˆ 0 (unit outward normal fng), respectively are prescribed. On the side 1± 2, the prescribed displacement varies linearly from 0 at Point 1 to u0 at 2. The medium is initially at rest. For the corresponding unbounded medium, a truncated semi-in®nite wedge is examined (Fig. 7b). The same displacement boundary condition as for the bounded medium is prescribed on the side 1±2. The same also applies to the sides extending from Points 1 and 2 to in®nity.

4. Scaled-boundary-transformation-based derivation In the ®nite-element method, the equations of motion of a single element are derived ®rst which leads to its static-sti€ness and mass matrices. Assemblage of the ®nite elements, enforcing equilibrium and compatibility then yields the equations of motion of the global system with its corresponding property matrices. The derivation of the scaled boundary ®nite-element equations is analogous. 4.1. Scaled boundary transformation of geometry

Fig. 7. Out-of-plane motion of (a) wedge and (b) truncated semi-in®nite wedge.

A single line ®nite element 1±2 with its corresponding triangle O±1±2 is ®rst addressed for the illustrative example (Fig. 8a) (as the illustrative example is solved with a single line ®nite element in the accompanying paper [9]). The scaling centre O coincides with the origin of the Cartesian coordinate system. Moving from 1 to 2, O must be on the left. The coordinates of a point on the line element are denoted as x, y. Two nodes with the coordinates (x1 , y1 ) and (x2 , y2 ) are arranged as

196

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

This equation de®nes the transformation from the Cartesian coordinates x^, y^ to the scaled boundary coordinates n, g, which forms a right-handed system. Lines S n with a constant n are parallel to the line ®nite element S e (Fig. 8a). Lines S g with a constant g pass through the scaling centre O and the point on the line ®nite element with the same g. The spatial derivatives in the two coordinate systems are related as     o=on o=o^ x ˆ ‰J^…n; g†Š …10† o=og o=o^ y with the Jacobian matrix de®ned as   x^;n y^;n ‰J^…n; g†Š ˆ : x^;g y^;g

…11†

The partial derivatives of x^, y^ with respect to n, g are determined using Eqs. (9a) and (9b) Fig. 8. Two-node line ®nite element: (a) scaled boundary transformation and (b) parent element.



 x1 ; fxg ˆ x2   y1 fyg ˆ : y2

…4a† …4b†

The line element is generated from its parent element (Fig. 8b) with the mapping functions   ‰N Š ˆ 12 …1 ÿ g† 12 …1 ‡ g† ; …5† as x ˆ ‰N Šfxg ˆ x ‡ 12Dx g; y ˆ ‰N Šfyg ˆ y ‡ 12Dy g

…6a† …6b†

with the abbreviations Dx ˆ x2 ÿ x1 ;

…7a†

D y ˆ y2 ÿ y1 ;

…7b†

x ˆ 12…x1 ‡ x2 †; y ˆ 12…y1 ‡ y2 †:

…8b†

…8a†

The dimensionless radial coordinate n points from O to a point on the boundary (with n ˆ 1). The scaled boundary transformation relates any point in the domain with the coordinates denoted as x^, y^ to the corresponding point on the line element as (Fig. 8a)   x^ ˆ nx ˆ n x ‡ 12Dx g ; …9a†   y^ ˆ ny ˆ n y ‡ 12Dy g : …9b†

x^;n ˆ x ˆ x ‡ 12Dx g; y^;n ˆ y ˆ y ‡ 12Dy g; x^;g ˆ nx;g ˆ 12Dx n; y^;g ˆ ny;g ˆ 12Dy n: ‰J^…n; g†Š can be written as   1 ‰J^…n; g†Š ˆ ‰J …g†Š; n where



x ‰J …g†Š ˆ x;g

y y;g

…12a† …12b† …12c† …12d†

…13†

 …14†

and with (Eqs. (6a) and (6b)) x;g ˆ ‰N Š;g fxg ˆ 12Dx ; y;g ˆ ‰N Š;g fyg ˆ 12Dy :

…15a† …15b†

‰J …g†Š is a function of the geometry of the line element only. The argument g is omitted for conciseness. The derivatives with respect to x^, y^ follow from Eq. (10)     o=o^ x o=on ˆ ‰J^…n; g†Šÿ1 : …16† o=o^ y o=og Substituting Eqs. (13) and (14) into Eq. (16) yields   o 1 o o=o^ x ˆ fb1 g ‡ fb2 g ; …17† o=o^ y on n og where

   y;g Dy 1 ˆ ; x1 y2 ÿ x2 y1 ÿDx ÿx;g ( )   ÿ2 y ÿ Dy g 1 ÿy 1 2 fb g ˆ ˆ jJ j x x 1 y2 ÿ x 2 y1 2x ‡ Dx g

fb1 g ˆ

1 jJ j



…18a† …18b†

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

with the determinant (Eq. (14)) jJ j ˆ

1 …x y 2 1 2

ÿ x2 y1 †:

…19†

2

1

…fb gjJ j†;g ˆ ÿfb gjJ j

…20†

applies. The in®nitesimal area dV e of the domain is calculated as (Eq. (13)) dV e ˆ jJ^…n; g†jdn dg ˆ njJ jdn dg:

corresponds to either n ˆ 1 or g ˆ 1 (Fig. 9). The surface traction is calculated as T sn ˆ fng fsg

For later use, note that

…21†

The in®nitesimal length of a line S n with constant n (Fig. 8a) equals (Eqs. (12c) and (12d)) q q dS n ˆ …^ x;g †2 ‡ …^ y ;g †2 dg ˆ 12n D2x ‡ D2y dg: …22† To calculate the surface traction, the unit normal vectors fnn g and fng g of the lines S n and S g with constant n and g, respectively, are addressed (Fig. 8a)   Dy 1 ; …23a† fnn g ˆ q D2x ‡ D2y ÿDx   ÿy 1 fng g ˆ p : …23b† x x2 ‡ y 2

197

…27†

with fng denoting the unit outward normal of the side. A side is positive or negative, when the outward normal forms an angle in absolute value smaller than 90° with the positive or negative direction of the corresponding coordinate axis, respectively. For the bounded and unbounded cases, the sides 1±2 corresponding to the line elements are positive and negative, respectively. The sides g ˆ ‡1 and g ˆ ÿ1 are positive and negative, respectively, for both the bounded and unbounded cases. The normals calculated in Eqs. (23a) and (23b) correspond to the outward normals of the positive sides. For the outward normals of the negative sides sign changes occur.

Note that fnn g and fng g are proportional to fb1 g and fb2 g, respectively (Eqs. (18a), (18b), (23a), and (23b)) q q D2x ‡ D2y D2x ‡ D2y 1 n fb g ˆ fnn g; fn g ˆ …24a† x 1 y2 ÿ x 2 y1 2jJ j

fb2 g ˆ

p p 2 x2 ‡ y 2 g x2 ‡ y 2 g fn g ˆ fn g: x 1 y2 ÿ x 2 y1 jJ j

…24b†

4.2. Governing equations in scaled boundary coordinates The scaled boundary transformation is applied to the geometry of the domain only. The shear stress sx , sy are still de®ned in the original Cartesian coordinate system x^, y^. The equilibrium equation (Eq. (1)) is reformulated in the scaled boundary coordinates as follows. With the shear stresses and Eq. (17) 1 fb1 gT fsg;n ‡ fb2 gT fsg;g ÿq u‡p ˆ0 n

…25†

results. The stress±displacement relationship (Eq. (2)) is transformed to   1 fsg ˆ G fb1 gu;n ‡ fb2 gu;g : …26† n The boundary conditions are straightforwardly formulated in the scaled boundary coordinate system, as each side (line element and two side-faces) of the triangle

Fig. 9. (a) Bounded medium with positive side 1±2 (n ˆ 1), positive side-face O±2 (g ˆ ‡1) and negative side-face O±1 (g ˆ ÿ1) and (b) unbounded medium with negative side 1±2 (n ˆ 1), positive side-face 2±1 (g ˆ ‡1) and negative side-face 1±1 (g ˆ ÿ1).

198

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

4.3. Weighted-residual technique The triangle de®ned by scaling a single line element on the boundary is addressed (Fig. 8a). The superscript e is dropped for conciseness. To derive a ®nite-element approximation, the weighted-residual technique is applied to the equilibrium equation in the scaled boundary coordinates (Eq. (25)) (strong form). Multiplication by a weighting function w ˆ w…n; g† and integration over the domain of the triangle V yields Z Z 1 wfb1 gT fsg;n dV ‡ w fb2 gT fsg;g dV n V V Z Z ÿ wq u dV ‡ wp dV ˆ 0 …28† V

V

with dV speci®ed in Eq. (21). The second term of Eq. (28) denoted as I is examined: Z 1 Z ‡1 Iˆ wfb2 gT fsg;g jJ jdg dn …29† 0

ÿ1

holds (Fig. 8a). Applying integration by parts (Green's theorem) leads to Z 1   gˆ‡1 Iˆ wfb2 gT jJ jfsg gˆÿ1 0  Z ‡1   ÿ wfb2 gT jJ j ;g fsg dg dn: …30† ÿ1

Substituting Eqs. (24b) and (20) results in Z 1  p  Iˆ w x2 ‡ y 2 fng gfsg gˆ‡1 0  p  g ÿ w x2 ‡ y 2 fn gfsg gˆÿ1  Z ‡1   1 T ÿ wfb g ‡ w;g fb2 gT fsgjJ jdg dn: ÿ

After substituting Eq. (21), Eq. (28) is satis®ed by setting the integrand of the integral over n equal to zero. This corresponds to enforcing the equation of equilibrium (Eq. (25)) exactly in the n-direction. Using Eq. (31) Z ‡1 q 1 T 2 2 n wfb g fsg;n jJ jdg ‡ x2 ‡ y2 …wsn † Z ÿ gˆÿ1

‡1 ÿ1

Z  ‡ w;g fb2 gT fsgjJ jdg ÿ n Z ‡n

‡1

ÿ1

‡1

ÿ1

wpjJ jdg ˆ 0



u…n; g† ˆ ‰N…g†Šfu…n†g with the displacement (Fig. 10)   u1 …n† : fu…n†g ˆ u2 …n†

…33†

…34†

with (Eqs. (18a),(18b) and (5)) …31†

q 2 2 ‡ x1 ‡ y1 …wsn †

The displacement of the line element on the boundary n ˆ 1 is interpolated using the mapping functions (Eq. (5)) as the shape functions ‰N …g†Š. It is postulated that the same shape functions apply for all lines with a constant n

The shear stresses are calculated substituting Eq. (33) into Eq. (26)   1 fsg ˆ G ‰B1 Šfu…n†g;n ‡ ‰B2 Šfu…n†g …35† n

ÿ1

ÿ1

Fig. 10. Displacement in domain.

gˆ‡1

ÿ wfb1 gT

wq ujJ jdg

‰B1 Š ˆ fb1 g‰N Š ˆ

1 2…x1 y2 ÿ x2 y1 †   Dy …1 ÿ g† Dy …1 ‡ g†  ; ÿDx …1 ÿ g† ÿDx …1 ‡ g† …36a†

‰B2 Š ˆ fb2 g‰N Š;g ˆ

1 2…x1 y2 ÿ x2 y1 † " # 2 y ‡ Dy g ÿ2 y ÿ Dy g  ; ÿ2x ÿ Dx g 2x ‡ Dx g

…36b†

‰B1 Š and ‰B2 Š are independent of n. The weighting function is discretized in the same way as the displacement is …32†

is obtained. Note that no integrals over the domain are present in Eq. (32).

w…n; g† ˆ ‰N …g†Šfw…n†g:

…37†

Substituting Eq. (37) into Eq. (32) and using Eqs. (36a) and (36b) yields

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

 Z fw…n†gT n 

‡1

ÿ1

‰B1 ŠT fsg;n jJ jdg ÿ

‡ ‰B2 ŠT fsgjJ jdg ÿ n  ‡ fF t g ˆ 0 with

Z

b

‡1

fF g ˆ

Z

‡1

ÿ1

Z

‡1



ÿ1

yields

ÿ ‰B1 ŠT

  T ‰E0 Šn2 fu…n†g;nn ‡ ‰E0 Š ÿ ‰E1 Š ‡ ‰E1 Š nfu…n†g;n

‰N ŠT q ujJ jdg ‡ nfF b g

ÿ ‰E2 Šfu…n†g ÿ ‰M 0 Šn2 f u…n†g ‡ fF …n†g ˆ 0 …38†

T

…39† …40†

gˆ‡1

fF b g, which can be a function of n, represents the nodal loads resulting from the body load. fF t g, which can again be a function of n, corresponds to the nodal forces resulting from the surface tractions acting on the sidefaces A. As fw…n†g is arbitrary, its coecient in Eq. (38) vanishes. Substituting the stresses (Eq. (35)) and the displacement (Eq. (33)) leads to the di€erential equation expressed in the displacement fu…n†g  Z ‡1 1 n ‰B1 ŠT G ‰B1 Šfu…n†g;nn ‡ ‰B2 Šfu…n†g;n n ÿ1  Z ‡1  1 2 ÿ ‰B1 ŠT ÿ 2 ‰B Šfu…n†g jJ jdg ÿ n ÿ1    1 ‡ ‰B2 ŠT G ‰B1 Šfu…n†g;n ‡ ‰B2 Šfu…n†g jJ jdg n Z ‡1 ÿn ‰N ŠT q‰N Šfu…n†gjJ jdg ‡ nfF b g ‡ fF t g ˆ 0: ÿ1

…41† Multiplying Eq. (41) by n and introducing the coecient matrices, Z

Z

D2x ‡ D2y ‰B1 ŠT G‰B1 ŠjJ jdg ˆ G 6…x1 y2 ÿ x2 y1 † ÿ1   2 1  ; 1 2 ‡1

…42a†

D2x ‡ D2y ‰B2 ŠT G‰B1 ŠjJ jdg ˆ G 12…x 1 y2 ÿ x2 y1 † ÿ1     ÿ1 ÿ1 ÿ1 1 xDx ‡ yDy ;  ÿG 2…x1 y2 ÿ x2 y1 † 1 1 1 ÿ1   Z ‡1 D2x ‡ D2y ‡ 12 x2 ‡ y2 ‰E2 Š ˆ ‰B2 ŠT G‰B2 ŠjJ jdg ˆ G 12…x1 y2 ÿ x2 y1 † ÿ1   1 ÿ1  ÿ1 1

‰E1 Š ˆ

and ‰M 0 Š ˆ

Z

‡1 ÿ1

‰N ŠT q‰N ŠjJ jdg ˆ q



x1 y2 ÿ x2 y1 2 1 1 2 6

…42b†

…42c†

…43†

…46†

Eq. (46) represents the scaled boundary ®nite-element equation in displacement formulated in the time domain. The corresponding equation in the frequency domain equals   ‰E0 Šn2 fu…n†g;nn ‡ ‰E0 Š ÿ ‰E1 Š ‡ ‰E1 ŠT  nfu…n†g;n ÿ‰E2 Šfu…n†g ‡ x2 ‰M 0 Šn2 fu…n†g ‡ fF …n†g ˆ 0



…45†

Note that the coecient matrices ‰E0 Š, ‰E1 Š, ‰E2 Š, ‰M 0 Š are independent of n. Integrations over the line element (at n ˆ 1) only are involved. The expressions for ‰E0 Š, ‰E1 Š, ‰E2 Š are similar to that for the static-sti€ness matrix of a standard line ®nite element, and ‰M 0 Š resembles the mass matrix. Integrations of polynomials only in g occur. ‰E0 Š and ‰M 0 Š are positive de®nite, ‰E2 Š is semi-positive de®nite, and ‰E1 Š is non-symmetric. The derivation of Eq. (44) is performed for the bounded medium (Fig. 9a). Eq. (44) also applies to the unbounded medium (Fig. 9b) with the same de®nition of the coecient matrices ‰E0 Š, ‰E1 Š, ‰E2 Š, ‰M 0 Š (Eqs. (42a)± (42c) and (43)). The only modi®cation in the derivation consists of replacing the domain of the triangle V (0 6 n 6 1) in Eq. (28) by that of the corresponding ``trapezoid'' (1 6 n < 1). Eq. (44) applies to the domain of the triangular (or trapezoid) corresponding to one line ®nite element on the boundary. To model the total domain, an assemblage as in the conventional ®nite-element method is performed and the displacement boundary conditions on the side-face Au , if present, are enforced (Fig. 4b). To simplify the nomenclature, the same symbols are used for the assembled coecient matrices, the assembled displacements and the assembled nodal loads in the following. In the assemblage process, the nodal forces fF t g in Eq. (45) will cancel with the exception of those on the side-face At of the total boundary with prescribed surface tractions (Fig. 4b). The assemblage process yields   ‰E0 Šn2 fu…n†g;nn ‡ ‰E0 Š ÿ ‰E1 Š ‡ ‰E1 ŠT nfu…n†g;n ÿ ‰E2 Šfu…n†g ÿ ‰M 0 Šn2 f u…n†g ‡ fF …n†g ˆ 0:

‡1

…44†

with the contribution of the body load and the surface traction fF …n†g ˆ n2 fF b g ‡ nfF t g:

‰N Š pjJ jdg; 8ÿ1p 9 > > < x21 ‡ y12 sn = gˆÿ1 fF t g ˆ p > x2 ‡ y 2 s > : ; 2 2 n

‰E0 Š ˆ

199

…47†

with the amplitudes of the displacements fu…n†g and of the nodal loads fF …n†g. Eq. (47) is a system of linear

200

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

non-homogeneous second-order ordinary di€erential equations for the displacement amplitudes fu…n†g with the dimensionless radial coordinate n as the independent variable. As already pointed out, the coecient matrices ‰E0 Š, ‰E1 Š, ‰E2 Š, ‰M 0 Š are independent of n.

To determine the nodal force±displacement relationship as in standard ®nite elements, the internal nodal forces fQ…n†g, which are equal to the resultants of the stresses, on a line with a constant n are addressed. Applying the principle of virtual work yields Z fw…n†gT fQ…n†g ˆ wfnn gT fsg dS n : …48† Sn

Substituting Eqs. (22), (24a) and (37) leads for an arbitrary fw…n†g to Z ‡1 fQ…n†g ˆ n ‰N ŠT fb1 gT fsgjJ jdg: …49† ÿ1

…50†

Comparing with Eqs. (42a) and (42b) yields fQ…n†g ˆ ‰E0 Šnfu…n†g;n ‡‰E1 ŠT fu…n†g:

…51†

For a bounded medium (Fig. 9a), side 1±2 is positive, i.e. fnn g is the unit outward normal. The nodal forces fR…n†g are equal to the internal nodal forces fQ…n†g. For the unbounded medium (Fig. 9b), the opposite applies. The two cases can be uni®ed as fR…n†g ˆ fQ…n†g

…52†

with the upper and lower signs corresponding to the bounded and unbounded media, respectively. The same Eqs. (51) and (52) apply to the assembled system. In the frequency domain, the amplitudes of the displacements fu…n†g are related to those of the nodal forces fR…n†g as fR…n†g ˆ ‰S…x; n†Šfu…n†g ÿ fRF …n†g:

…53†

‰S…x; n†Š denotes the dynamic-sti€ness matrix on a line with a constant n and fRF …n†g represents the amplitudes of the nodal loads due to the body load and the surface traction. Substituting Eqs. (51) and (53) into Eq. (52) yields  ‰S…x; n†Šfu…n†g  RF …n† ˆ ‰E0 Šnfu…n†g;n ‡‰E1 ŠT fu…n†g: Di€erentiating Eq. (54) with respect to n

…55†

is obtained. Adding Eq. (47) and Eq. (55) multiplied by n results in

4.4. Dynamic-sti€ness matrix

Substituting Eqs. (35) and (36a) results in  Z ‡1 fQ…n†g ˆ n ‰B1 ŠT G ‰B1 Šfu…n†g;n ÿ1  1 ‡ ‰B2 Šfu…n†g jJ jdg: n

 ‰S…x; n†Š;n fu…n†g  ‰S…x; n†Šfu…n†g;n   RF …n† ;n ÿ‰E0 Šnfu…n†g;nn   ÿ ‰E0 Š ‡ ‰E1 ŠT fu…n†g;n ˆ 0

…54†

n‰S…x; n†Š;n fu…n†g ‡ …‰S…x; n†Š  ÿ ‰E1 Š†nfu…n†g;n n RF …n† ;n ÿ‰E2 Šfu…n†g ‡ x2 ‰M 0 Šn2 fu…n†g ‡ fF …n†g ˆ 0:

…56†

Solving Eq. (54) for nfu…n†g;n and substituting in Eq. (56) yields ÿ  
 ‰S…x; n†Š ÿ ‰E1 Š ‰E0 Šÿ1  ‰S…x; n†Š ÿ ‰E1 ŠT   n‰S…x; n†Š;n ÿ ‰E2 Š ‡ x2 n2 ‰M 0 Š fu…n†g  ÿ  n RF …n† ;n   ‰S…x; n†Š
 ÿ ‰E1 Š ‰E0 Šÿ1 RF …n† ‡ fF …n†g ˆ 0: …57† For arbitrary fu…n†g, its coecient matrix must vanish leading to   ÿ
 ‰S…x; n†Š ÿ ‰E1 Š ‰E0 Šÿ1  ‰S…x; n†Š ÿ ‰E1 ŠT  n‰S…x; n†Š;n ÿ‰E2 Š ‡ x2 n2 ‰M 0 Š ˆ 0:

…58†

which will permit ‰S…x; n†Š to be determined. In addition, the remaining part of Eq. (57) leads to   n RF …n† ;n ‡…‰S…x; n†Š ÿ ‰E1 Š†‰E0 Šÿ1 RF …n† ˆ fF …n†g

…59†

with ‰S…x; n†Š determined from Eq. (58). This is a system of linear ordinary di€erential equations of the ®rst order for the amplitudes of the nodal load fRF …n†g. The boundary condition is fRF …n ˆ 0†g ˆ 0 for a bounded medium, and fRF …n ! 1†g ˆ 0 for an unbounded medium. As n‰S…x; n†Š;n ˆ xn‰S…x; n†Š;xn applies and ‰E0 Š, ‰E1 Š, ‰E2 Š, ‰M 0 Š are independent of x and n, it follows from Eq. (58) that ‰S…x; n†Š is actually a function of the product xn only. This con®rms the well-known fact that ‰S…x; n†Š, is a function of the dimensionless frequency de®ned for the line with a constant n as aˆ

xnr0 cs

…60†

p with the shear-wave velocity cs ˆ G=q and the characteristic length r0 of the boundary (n ˆ 1). The derivative with respect to a can be interpreted either for varying n with ®xed x or for varying x with ®xed n. Thus, for the two-dimensional case ‰S…x; n†Š ˆ ‰S…a†Š satis®es

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

a‰S…a†Š;a ˆ n‰S…x; n†Š;n ˆ x‰S…x; n†Š;x :

201

…61†

This permits the spatial derivative to be replaced by the frequency derivative in Eq. (58) ÿ

 
 ‰S…x; n†Š ÿ ‰E1 Š ‰E0 Šÿ1  ‰S…x; n†Š ÿ ‰E1 ŠT  x‰S…x; n†Š;x ÿ‰E2 Š ‡ x2 n2 ‰M 0 Š ˆ 0:

…62†

On the boundary (n ˆ 1), the dynamic-sti€ness matrix for the bounded medium ‰S b …x†Š (superscript b for bounded) is expressed as ÿ

 
‰S b …x†Š ÿ ‰E1 Š ‰E0 Šÿ1 ‰S b …x†Š ÿ ‰E1 ŠT ‡ x‰S b …x†Š;x ÿ‰E2 Š ‡ x2 ‰M 0 Š ˆ 0

…63†

Fig. 11. Concept of mechanically based derivation for a bounded medium.

and for the unbounded medium ‰S 1 …x†Š, (superscript 1 for unbounded) as ÿ

 
‰S 1 …x†Š ‡ ‰E1 Š ‰E0 Šÿ1 ‰S 1 …x†Š ‡ ‰E1 ŠT ÿ x‰S b …x†Š;x ÿ‰E2 Š ‡ x2 ‰M 0 Š ˆ 0:

…64†

These represent the scaled boundary ®nite-element equations in dynamic sti€ness formulated in the frequency domain for two-dimensional bounded and unbounded media. Each of them is a system of non-linear ®rst-order ordinary di€erential equations with the frequency x as the independent variable. For the bounded medium, ‰S b …x†Š is usually not calculated. The static-sti€ness matrix ‰K b Š, and the mass matrix ‰M b Š are used, which follow from the lowfrequency expansion of ‰S b …x†Š ‰S b …x†Š ˆ ‰K b Š ÿ x2 ‰M b Š ‡ O…x4 †:

…65†

5. Mechanically based derivation The dynamic characteristics of a bounded or an unbounded medium can be described by the force± displacement relationship with respect to the degrees of freedom of the nodes on the boundary (interface). In the frequency domain, the amplitudes of the displacements fug are related to those of the nodal forces fRg (for vanishing body loads) as fRg ˆ ‰S…x†Šfug

…66†

with the dynamic-sti€ness matrix ‰S…x†Š. As an alternative to the scaled-boundary-transformation-based derivation (Section 4), the scaled boundary ®nite-element equations in dynamic sti€ness and

displacement can also be established using a mechanically based derivation, which is more familiar to engineers. The concept is explained using the bounded case of the illustrative example (Fig. 7a), the wedge. The side 1±2 represents the boundary (interface) which is discretized (Fig. 11a). O denotes the similarity centre and 1±2 corresponds to n ˆ 1 with the characteristic length re (e for exterior). A similar ®ctitious boundary with the scaling factor ni < 1 and the characteristic length ri (ri ˆ ni re ) is constructed (i for interior). The discretization of the similar ®ctitious boundary is selected to be similar to that of the boundary. The wedge (referred to re ) is divided by the ®ctitious boundary into two parts: the wedge (Fig. 11c) with a similar boundary (referred to ri ) and a trapezoid with its exterior boundary coinciding with the boundary and its interior boundary with the ®ctitious boundary (Fig. 11b). The discretizations of the exterior and interior boundaries are identical to those of the boundary and the ®ctitious boundary of the wedge. The ®nite-element discretization of the trapezoid is de®ned by connecting the corresponding nodes on its exterior and interior boundaries. The discretized trapezoid is called the ®nite-element cell. Adding the ®nite-element cell (Fig. 11b) to the wedge referred to ri (Fig. 11c) results in the wedge referred to re (Fig. 11a). The same applies to their dynamic-sti€ness matrices. This assemblage enforcing compatibility and equilibrium leads to a relationship linking the dynamic-sti€ness matrices at the two similar boundaries of the wedge via that of the ®nite-element cell. Another relationship of the dynamic-sti€ness matrices at the boundaries of the wedge follows from similarity. Thus, the dynamic-sti€ness matrices at the two boundaries of the wedge can be expressed as a function of that of the ®nite-element cell. After performing the limit of the cell width re ÿ ri ! 0

202

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

ÿ n1 ‡ n2 ÿ n3 ÿ 3n4 ˆ 0; n1 ‡ n3 ‡ n4 ˆ 0;

…69a† …69b†

ÿ 2n1 ÿ 2n3 ÿ n5 ˆ 0:

…69c†

The rank of the coecient matrix of the system of three equations (Eq. (69a)) with ®ve unknowns equals 3, which permits the two unknowns n1 , n5 to be chosen arbitrarily. For n1 ˆ 1, n5 ˆ 0, the other unknowns are n2 ˆ 0, n3 ˆ ÿ1, n4 ˆ 0, yielding the ®rst dimensionless variable ‰SŠGÿ1 . For n1 ˆ 0, n5 ˆ 1, the other unknowns are n2 ˆ 1, n3 ˆ ÿ0:5, n4 ˆ 0:5, resulting in the second dimensionless variable rGÿ0:5 q0:5 x, which is the dimensionless frequency, a0 ˆ Fig. 12. Concept of mechanically based derivation for an unbounded medium.

analytically, an equation for the dynamic-sti€ness matrix of the wedge is obtained. The same concept also applies to the truncated semiin®nite wedge (Fig. 12). Note that the boundary and the similar ®ctitious boundary corresponds to ri and re , respectively. The scaling factor ne is larger than 1. 5.1. Dynamic-sti€ness matrices at similar boundaries of a medium The dynamic-sti€ness matrices of similar boundaries with the characteristic length r ‰S…x; r†Š are addressed. A dimensional analysis is performed to identify the independent dimensionless variables of which ‰S…x; r†Š is a function. The characteristic length r, the shear modulus G, the mass density q and the frequency x are sucient to determine ‰S…x; r†Š for the out-of-plane motion in two dimensions. A bracket denotes the dimension; ‰LŠ, ‰MŠ and ‰TŠ are the dimensions of length, mass and time. The dimensions of the variables are ÿ1

ÿ2

‰‰SŠŠ ˆ ‰LŠ ‰MŠ‰TŠ ; ‰rŠ ˆ ‰LŠ; ÿ1

…67a† …67b†

ÿ2

‰GŠ ˆ ‰LŠ ‰MŠ‰TŠ ;

…67c†

ÿ3

…67d†

ÿ1

…67e†

‰qŠ ˆ ‰LŠ ‰MŠ; ‰xŠ ˆ ‰TŠ :

The product of all these variable raised to unknown powers ni (i ˆ 1; 2; . . . ; 5) must be dimensionless ‰‰SŠŠn1 ‰rŠn2 ‰GŠn3 ‰qŠn4 ‰xŠn5 ˆ ‰LŠÿn1 ‡n2 ÿn3 ÿ3n4  ‰MŠn1 ‡n3 ‡n4  ‰TŠÿ2n1 ÿ2n3 ÿn5 : This results in

…68†

xr cs

…70†

p with the shear-wave velocity cs ˆ G=q. Note that a0 varies with r. The ®rst dimensionless variable, ‰SŠGÿ1 will be a function ‰SŠ of the second dimensionless variable, a0 yielding ‰S…x; r†Š ˆ G‰S…a0 †Š

…71†

x and r do not appear explicitly in ‰S…a0 †Š, but only as a product in a0 (Eq. (70)). The variation of the dynamic-sti€ness matrix, as a function of the characteristic length r and of the excitation frequency x, is studied. Eq. (71) is addressed. ‰S…a0 †Š is a function of one independent variable, i.e. a0 . For constant x, ‰S…a0 †Š is a function of r, or just as valid, for constant r, ‰S…a0 †Š is a function of x. The same change in a0 can be achieved by varying the values of either r or x with the other ®xed. The derivative thus follows for a constant x varying r as ‰S…a0 †Š;a0 ˆ

cs ‰S…a0 †Š;r x

…72a†

and the same result is calculated for a constant r but varying x as ‰S…a0 †Š;a0 ˆ

cs ‰S…a0 †Š;x : r

…72b†

Equating the two right-hand sides of Eqs. (72a) and (72b) yields r‰S…a0 †Š;r ˆ x‰S…a0 †Š;x :

…73†

The partial derivative of ‰S…a0 †Š with respect to r can thus be replaced by that with respect to x for similar boundaries. Eqs. (71) and (73) lead to r‰S…r; x†Š;r ˆ x‰S…x; r†Š;x ; which is applicable to the two-dimensional case.

…74†

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

5.2. Coecient matrices The static-sti€ness and mass matrices of a single ®nite element of the cell are calculated, which determine the coecient matrices of the scaled boundary ®nite-element equations. A single 2-node line ®nite element 1±2 on the interior boundary i is addressed for the illustrative example (Fig. 13) (as the illustrative example is solved with a single line ®nite element in the accompanying paper). Moving from 1 to 2, O must be on the left. Starting from this line element a two-dimensional 4-node ®nite element of the cell with its exterior boundary e obtained by scaling the interior boundary with a factor larger than 1 is constructed. A single element in the radial direction is introduced. The coordinates of the two nodes of the line element are arranged as   x1 fxg ˆ ; …75a† x2   y1 fyg ˆ : …75b† y2 On the interior boundary ‰N Š ˆ ‰N …g†Š denotes the mapping functions of the line parent element (Fig. 13c) ‰N Š ˆ ‰ 12 …1 ÿ g†

1 …1 2

‡ g† Š:

…76†

203

The coordinates of a point on the line element are generated as x ˆ ‰N Šfxg ˆ x ‡ 12Dx g; y ˆ ‰N Šfyg ˆ y ‡ 12Dy g

…77a† …77b†

with the abbreviations Dx ˆ x2 ÿ x1 ; Dy ˆ y2 ÿ y1 ;

…78a† …78b†

x ˆ 12…x1 ‡ x2 †; y ˆ 12…y1 ‡ y2 †:

…79a† …79b†

The coordinates of the nodes on the exterior boundary are expressed as fxe g ˆ …1 ‡ w†fxg; fye g ˆ …1 ‡ w†fyg

…80a† …80b†

with the dimensionless cell width, wˆ

re ÿ ri : ri

…81†

The mapping functions of the two-dimensional parent ^ g†Š (Fig. 13b) are generated from element ‰N^ Š ˆ ‰N^ …n; those of the line element (Eq. (76)) with the local coor^ g. The symbol circum¯ex ^ denotes the twodinates n, dimensional ®nite element. The mapping functions are decomposed with respect to the interior and exterior boundaries:   ^ ˆ ‰N^i Š ‰N^e Š ; ‰NŠ …82† where ^ Š ‰N^j Š ˆ 12…1 ‡ n^j n†‰N

…j ˆ i; e†

…83†

with n^j ˆ ÿ1 for j ˆ i and n^j ˆ ‡1 for j ˆ e. The coordinates x^, y^ of a point in the two-dimensional ®nite element are expressed as   w ^ x; x^ ˆ ‰N^i Šfxg ‡ ‰N^e Šfxe g ˆ 1 ‡ …1 ‡ n† …84a† 2   w ^ y: y^ ˆ ‰N^i Šfyg ‡ ‰N^e Šfye g ˆ 1 ‡ …1 ‡ n† …84b† 2 The Jacobian matrix of the two-dimensional ®nite element equals   ^ g†Š ˆ x^;n^ y^;n^ ‰J^…n; x^;g y^;g w  ˆ 2 …85† ^ ‰J …g†Š; 1 ‡ w2 …1 ‡ n† where Fig. 13. (a) Line ®nite element 1±2 with adjacent two-dimensional ®nite element of cell; (b) two-dimensional parent element of ®nite element of cell and (c) line parent element.

‰J …g†Š ˆ



x x;g

y y;g

 …86†

204

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

and with

‰B1 Š ˆ fb1 g‰N Š ˆ

x;g ˆ ‰N Š;g fxg ˆ 12Dx ;

…87a†

y;g ˆ ‰N Š;g fyg ˆ 12Dy :

…87b†

‰J …g†Š is a function of the geometry of the line element only. The argument g is omitted for conciseness. The determinant of the Jacobian matrix equals   ^ g†j ˆ w 1 ‡ w …1 ‡ n† ^ jJ j jJ^…n; …88† 2 2 with jJ j ˆ xy;g ÿx;g y ˆ 12…xDy ÿ Dx y† ˆ 12…x1 y2 ÿ x2 y1 †:

…89†

The inverse of the Jacobian matrix equals & % 2 ÿ1 ÿ1 w ^ g†Š ˆ ‰J Š ‰J^…n; 1

…90†

^ 1‡w2 …1‡n†

with ‰J Šÿ1 ˆ

 1 y;g jJ j ÿx;g

 ÿy : x

…91†

The derivatives with respect to x^, y^ follow from     x o=on^ ˆ ‰J^…n; ^ g†Š o=o^ o=o^ y o=og as 

o=o^ x o=o^ y



  ^ g†Šÿ1 o=on^ ˆ ‰J^…n; o=og o 1 o ˆ fb1 g ‡ fb2 g ; ^ og on^ 1 ‡ w …1 ‡ n†

…92†

‰B2 Š ˆ fb2 g‰N Š;g ˆ





Substituting Eqs. (96) and (95) into Eq. (98) yields Z ‡1   n^j n^l 0 w ^ dn^ ‰Kjl Š ˆ ‰E Š 1 ‡ …1 ‡ n† 2 2w ÿ1 Z ‡1   ^ n^j n ‡ l ‰E1 Š 1 ‡ n^j n^ dn^ ‡ ‰E1 ŠT 4 4 ÿ1 Z ‡1   w  1 ‡ n^l n^ dn^ ‡ ‰E2 Š 8 ÿ1    Z ‡1 1 ‡ n^ n^ 1 ‡ n^ n^ j l  dn^ …99† ^ 1 ‡ w …1 ‡ n† ÿ1 2

…93†

with the coecient matrices Z ‡1 ‰E0 Š ˆ ‰B1 ŠT G‰B1 ŠjJ jdg ÿ1



y;g Dy 1 ˆ ; x1 y2 ÿ x2 y1 ÿDx ÿx;g ( )   ÿ2 y ÿ Dy g 1 ÿy 1 fb2 g ˆ ˆ : jJ j x x1 y2 ÿ x2 y1 2x ‡ Dx g fb1 g ˆ

1 jJ j

…94a† …94b†

The in®nitesimal area dV of the ®nite element is calculated as   ^ g†jdndg ^ ˆ w 1 ‡ w …1 ‡ n† ^ jJ jdn^dg: dV ˆ jJ^…n; …95† 2 2 The shape functions of the ®nite element are chosen to be identical to the mapping functions (Eq. (82)). Applying Eq. (93) to the shape functions (Eq. (83)) leads to strain-nodal displacement matrix, ‰Bj Š ˆ

n^j 1 1 ‡ n^j n^  ‰B2 Š ‰B Š ‡  w w ^ 2 1 ‡ …1 ‡ n† 2

with

…j ˆ i; e†

…97b†

V

ˆG 

1 2…x1 y2 ÿ x2 y1 † " # y ÿ Dy g 2 y ‡ Dy g ÿ2  : ÿ2x ÿ Dx g 2x ‡ Dx g

…97a†

The static-sti€ness matrix ‰KŠ is decomposed into submatrices with respect to the interior and exterior boundaries as (j ˆ i; e; l ˆ i; e) Z ‰Kjl Š ˆ ‰Bj ŠT G‰Bl Š dV : …98†

2

where

1 2…x1 y2 ÿ x2 y1 †   Dy …1 ‡ g† Dy …1 ÿ g†  ; ÿDx …1 ÿ g† ÿDx …1 ‡ g†

…96†

‰E1 Š ˆ

Z

  D2x ‡ D2y 2 1 ; 6…x1 y2 ÿ x2 y1 † 1 2 ‡1

ÿ1

…100a†

‰B2 ŠT G‰B1 ŠjJ jdg

  D2x ‡ D2y ÿ1 1 12…x1 y2 ÿ x2 y1 † 1 ÿ1   ÿ1 ÿ1 xDx ‡ yDy ÿG ; 2…x1 y2 ÿ x2 y1 † 1 1 Z ‡1 ‰E2 Š ˆ ‰B2 ŠT G‰B2 ŠjJ jdg ˆG

…100b†

ÿ1

ˆG

  D2x ‡ D2y ‡ 12…x2 ‡ y2 † 1 ÿ1 : 12…x1 y2 ÿ x2 y1 † ÿ1 1

…100c†

Performing the integration over n^ analytically and decomposing with respect to the dimensionless cell width w yields ‰Kjl Š ˆ where

1 0 ‰K Š ‡ ‰Kjl1 Š ‡ w‰Kjl2 Š ‡ O…w2 †; w jl

…101†

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

‰Kjl0 Š ˆ n^j n^l ‰E0 Š; n^j n^l 0 n^j n^ ‰Kjl1 Š ˆ ‰E Š ‡ l ‰E1 Š ‡ ‰E1 ŠT ; 2 2! 2 ^j n^l n 1 ‰Kjl2 Š ˆ 1‡ ‰E2 Š: 4 3 The mass matrix ‰MŠ follows from Z ‰Mjl Š ˆ ‰N^j ŠT q‰N^l ŠdV : V

Substituting Eqs. (83) and (95) results in Z ‡1    w 0 ‰Mjl Š ˆ q ‰M Š 1 ‡ n^j n^ 1 ‡ n^l n^ 8 ÿ1   w ^ dn^  1 ‡ …1 ‡ n† 2 with the coecient matrix Z ‡1 ‰M 0 Š ˆ ‰N ŠT q‰N ŠjJ jdg ÿ1   x 1 y2 ÿ x 2 y1 2 1 ˆq : 1 2 6

…102a† …102b† …102c†

…103†

…104†

…105†

…106†

To derive the scaled boundary ®nite-element equation, the following relationships will be used: 0 ‰Kii0 Š ˆ ÿ‰Kie0 Š ˆ ÿ‰Kei0 Š ˆ ‰Kee Š ˆ ‰E0 Š; ÿ
T 1 ‰Kie1 Š ‡ ‰Kee Š ˆ ÿ ‰Kii1 Š ‡ ‰Kie1 Š ˆ ‰E1 Š;

…107b†

2 ‰Kii2 Š ‡ ‰Kie2 Š ‡ ‰Kei2 Š ‡ ‰Kee Š ˆ ‰E2 Š;

…107c†

‰Mii2 Š ‡ ‰Mie2 Š ‡ ‰Mei2 Š ‡ ‰Mee2 Š ˆ ‰M 0 Š:

Substituting Eqs. (83) and (95) results in

‡1



ÿ1

  w ^ dnfF ^ bg 1 ‡ n^j n^ 1 ‡ …1 ‡ n† 2

…108†

…110†

with fF b g ˆ

Z

‡1 ÿ1

‰N ŠT pjJ jdg:

…111†

Integrating Eq. (110) over n^ leads to    w w2 1 fPjb g ˆ fF b g: ‡ 1 ‡ n^j 2 3 4

…112†

Analogously, the nodal load fP t g of the ®nite element due to the prescribed surface traction sn (superscript t for traction), if present, on the two sides S gˆÿ1 , S gˆ‡1 (Fig. 13a) is calculated. Again decomposing with respect to the interior and exterior boundaries yields (j ˆ i; e) Z Z fPjt g ˆ ‰N^j ŠT sn dS g ‡ ‰N^j ŠT sn dS g : …113† S gˆ‡1

The in®nitesimal length for a constant g follows from Eqs. (84a) and (84b) as q w p ^ dS g ˆ …d^ x2 ‡ y 2 dn: x†2 ‡ …d^ y †2 ˆ …114† 2 Substituting Eqs. (83) and (114) results in fPjt g ˆ

w 4

Z

‡1



 1 ‡ n^j n^ dn^

ÿ1    ‰N ŠT sn 

gˆÿ1

 ‡ ‰N ŠT sn

…107a†

Eqs. (107a)±(107c) follows from Eqs. (102a)±(102c) with n^i ˆ ÿ1, n^e ˆ ‡1 and Eq. (108) from Eq. (106). The right-hand sides of Eqs. (107a)±(107c) and (108) are de®ned in Eqs. (100a)±(100c) and (105). ‰E0 Š, ‰E1 Š, ‰E2 Š and ‰M 0 Š are the coecient matrices of one line ®nite element. The nodal load fP b g of the ®nite element due to the body load p (superscript b for body load) is examined. Decomposing in subvectors with respect to the interior and exterior boundaries (Fig. 13a) yields (j ˆ i; e) Z fPjb g ˆ ‰N^j ŠT p dV …109† V

Z

S gˆÿ1

Performing the integration over n^ analytically yields ‰Mjl Š ˆ w‰Mjl2 Š ‡ O…w2 † ! n^j n^l w 1‡ ˆ ‰M 0 Š ‡ O…w2 †: 4 3

w 4

fPjb g ˆ

205

q x21 ‡ y12

q w x22 ‡ y22 ˆ fF t g gˆ‡1 2

…115†

with ( p ) x21 ‡ y12 sn jgˆÿ1 fF g ˆ p : x22 ‡ y22 sn jgˆ‡1 t

…116†

As for the static-sti€ness and mass matrices of the ®nite-element cell, the coecient matrices ‰E0 Š, ‰E1 Š, ‰E2 Š and ‰M 0 Š are assembled to form those of the boundary (Fig. 12) enforcing the displacement boundary condition on Au (Fig. 6b). The assembled coecient matrices will thus be banded. To simplify the nomenclature, the same symbols are used for the assembled coecient, staticsti€ness and mass matrices in the following. For the assembled static-sti€ness and mass matrices equations (107a)±(107c) and (108) still hold. The same assemblage also applies to the nodal loads fP b g, fP t g (Eqs. (112) and (115)) and the corresponding vectors fF b g, fF t g (Eqs. (111) and (116)).

206

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

5.3. Assemblage of ®nite-element cell and medium For both the bounded and unbounded cases (Figs. 11 and 12), the assemblage enforcing compatibility and equilibrium formulated in the frequency domain links the dynamic-sti€ness matrices of the medium at the boundary and at the similar ®ctitious boundary. This relationship involves the dynamic-sti€ness matrix of the ®nite-element cell, which can be expressed by its staticsti€ness and mass matrices. In addition, the nodal load of the medium due to the body load and the surface traction can be expressed by that of the ®nite-element cell. The force±displacement relationship of the ®niteelement cell located between the interior and exterior boundaries (Figs. 11b and 12b) is written as fP g ˆ ‰S…x†Šfug ÿ fPg

…117†

with the amplitudes of the displacements fug and of the nodal forces fP g. The dynamic-sti€ness matrix equals ‰S…x†Š ˆ ‰KŠ ÿ x2 ‰MŠ

…118†

with the static-sti€ness matrix ‰KŠ and the mass matrix ‰MŠ of the ®nite-element cell. The amplitudes of the known nodal load of the ®nite-element cell due to the body load and prescribed surface traction are formulated as (Eqs. (112) and (115)) fPg ˆ fP b g ‡ fP t g:

…119†

Partitioning Eq. (117) into the nodes lying on the interior and exterior boundaries yields        fPi g ‰Sii …x†Š ‰Sie …x†Š fui g fPi g ˆ ÿ : fPe g ‰Sei …x†Š ‰See …x†Š fue g fPe g …120† The force±displacement relationship of the medium at the boundaries (interfaces) corresponding to the interior and the exterior boundaries of the ®nite-element cell is formulated as fRi g ˆ ‰Si …x†Šfui g ÿ fRFi g;

…121a†

fRe g ˆ ‰Se …x†Šfue g ÿ fRFe g

…121b†

with the amplitudes of the unknown nodal loads of the medium fRF g due to the body load and prescribed surface traction. For the bounded medium, the outward normals of the boundary of the medium (Fig. 11a) and of the exterior boundary of the ®nite-element cell (Fig. 11b) point in the same direction yielding fRe g ˆ fPe g, while the opposite applies for the similar ®ctitious boundary of the medium (Fig. 11c) and the interior boundary of the ®nite-element cell (Fig. 11b) which leads to fRi g ˆ ÿfPi g. For the unbounded medium (Fig. 12), the opposite applies, resulting in fRi g ˆ fPi g and fRe g ˆ ÿfPe g. The two cases can be uni®ed as

fRi g ˆ  fPi g; fRe g ˆ  fPe g

…122a† …122b†

with the upper and lower signs corresponding to the bounded and unbounded media, respectively. Eliminating fP g and fRg by substituting Eqs. (120), (121a), and (121b) into Eqs. (122a) and (122b) yields      ‰Sii …x†Š ‰Sie …x†Š fui g fPi g ÿ ‰Sei …x†Š ‰See …x†Š fue g fPe g    fui g ‰Si …x†Š ˆ ‰Se …x†Š fue g   F fRi g ÿ : …123† fRFe g Solving the second row of Eq. (123) for fue g and substituting in the ®rst row results in …‰Si …x†Š ÿ ‰Sii …x†Š ‡ ‰Sie …x†Š…‰Se …x†Š ‡ ‰See …x†Š†ÿ1 ‰Sei …x†Š†fui g ˆ fPi g  fRFi g ÿ ‰Sie …x†Š…‰Se …x†Š ‡ ‰See …x†Š†ÿ1 …fPe g  fRFe g†:

…124†

As fui g is arbitrary, its coecient matrix must vanish leading to ‰Si …x†Š ˆ ‰Sii …x†Š ÿ ‰Sie …x†Š…‰Se …x†Š ‡ ‰See …x†Š†ÿ1 ‰Sei …x†Š:

…125†

This is the so-called dynamic condensation equation which links the dynamic-sti€ness matrices of the medium at the boundary and the similar ®ctitious boundary. In addition, the right-hand side of Eq. (124) is equal to zero yielding a relationship for the amplitudes of the nodal loads at the same two boundaries of the medium   RFi ˆ ÿfPi g ‡ ‰Sie …x†Š…‰Se …x†Š   …126† ‡ ‰See …x†Š†ÿ1 fPe g  fRFe g : 5.4. Analytical limit of ®nite-element cell width Eqs. (125) and (126) are valid for an arbitrary cell width. Performing the analytical limit of the cell width towards zero eliminates the discretization error in the radial direction. The dynamic condensation equation (Eq. (125)) is reformulated as …‰Se …x†Š ‡ ‰See …x†Š†‰Sie …x†Šÿ1 …‰Si …x†Š ÿ ‰Sii …x†Š† ‡ ‰Sei …x†Š ˆ 0:

…127†

The dynamic-sti€ness submatrices of the cell are written as follows. For instance, ‰Sii …x†Š ˆ ‰Kii Š ÿ x2 ‰Mii Š

…128†

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

applies with ‰Kii Š speci®ed in Eq. (101) and ‰Mii Š in Eq. (106). This yields ‰Sii …x†Š ˆ

1 0 ‰K Š ‡ ‰Kii1 Š ‡ w‰Kii2 Š ÿ wx2 ‰Mii2 Š ‡ O…w2 †: w ii …129†

1 ‰Kii1 Š ‡ ‰Kie1 Š ‡ ‰Kei1 Š ‡ ‰Kee Š |‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚} 1

ÿ

ÿ1 ÿ 1 ÿ w  ‰Se …x†Š ‡ ‰Kie1 Š ‡ ‰Kee Š ‰E0 Š  ‰Si …x†Š ÿ ‰Kii1 Š ÿ ‰Kie1 Š |‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚} 2 0 1 2  @ ‰Se …x†Š ÿ ‰Si …x†Š A ‡ w…‰Kii2 Š ‡ ‰Kie2 Š ‡ ‰Kei2 Š ‡ ‰Kee Š |‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚} |‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚} 3

Analogously

The sum identi®ed by 1 in Eq. (134) vanishes according to Eq. (107b). Eq. (107b) is substituted in the term identi®ed by 2. The sums identi®ed by 4 and 5 are transformed using Eqs. (107c) and (108), respectively. Dividing Eq. (134) by w results in …‰Se …x†Š ÿ ‰E1 Š†‰E0 Šÿ1 …‰Si …x†Š ÿ ‰E1 ŠT †

‰Sie …x†Š ˆ ÿ

 …130b†

…130c†

‰Se …x†Š ÿ ‰Si …x†Š ÿ ‰E2 Š ‡ x2 ‰M 0 Š ˆ O…w†: w

…130d†

…131†

The coecient matrix of w is equal to the inverse of that of 1=w in Eq. (130b). ‰AŠ and ‰BŠ follow from (Eqs. (130b) and (131)): ‰IŠ ˆ ‰Sie …x†Š‰Sie …x†Šÿ1 ÿ   ˆ ‰IŠ ÿ w ‰Kie1 Š‰E0 Šÿ1 ‡ ‰E0 Š‰AŠ ÿ w2 ‰Kie2 Š 
ÿ x2 ‰Mie2 Š ‰E0 Šÿ1 ÿ ‰Kie1 Š‰AŠ ‡ ‰E0 Š‰BŠ ‡ O…w3 †: …132† Setting the coecient matrices of w and w2 equal to zero yields ‰AŠ ˆ ÿ‰E0 Šÿ1 ‰Kie1 Š‰E0 Šÿ1 ; …133a†  ÿ 2
 0 ÿ1 0 ÿ1 1 0 ÿ1 1 2 2 ‰BŠ ˆ ÿ‰E Š ‰Kie Š‰E Š ‰Kie Š ÿ ‰Kie Š ÿ x ‰Mie Š ‰E Š : …133b†

Substituting Eqs. (130a)±(130d) and (131) with Eqs. (133a) and (133b) into Eq. (127) leads to

…135†

The limit of the dimensionless cell width w ! 0 can now be performed. With w de®ned in Eq. (81) lim

w!0

‰Se …x†Š ÿ ‰Si …x†Š ‰Se …x†Š ÿ ‰Si …x†Š ˆ lim ri re !ri w re ÿ ri ˆ r‰S…x†Š;r

follow. The inverse ‰Sie …x†Šÿ1 appearing in Eq. (127) is expressed as a polynomial in w with the unknown coecient matrices ‰AŠ and ‰BŠ ‰Sie …x†Šÿ1 ˆ ÿw‰E0 Šÿ1 ‡ w2 ‰AŠ ‡ w3 ‰BŠ ‡ O…w4 †:

…134†

5

1 0 ‰E Š ‡ ‰Kii1 Š ‡ w…‰Kii2 Š ÿ x2 ‰Mii2 Š† ‡ O…w2 †: w …130a†

ÿ
1 0 ‰E Š ‡ ‰Kie1 Š ‡ w ‰Kie2 Š ÿ x2 ‰Mie2 Š w ‡ O…w2 †; ÿ
1 ‰Sei …x†Š ˆ ÿ ‰E0 Š ‡ ‰Kei1 Š ‡ w ‰Kei2 Š ÿ x2 ‰Mei2 Š w ‡ O…w2 † ÿ 2
1 1 ‰See …x†Š ˆ ‰E0 Š ‡ ‰Kee Š ‡ w ‰Kee Š ÿ x2 ‰Mee2 Š w ‡ O…w2 †

4

ÿ x2 …‰Mii2 Š ‡ ‰Mie2 Š ‡ ‰Mei2 Š ‡ ‰Mee2 Š †† ˆ O…w2 †: |‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚}

Using Eq. (107a) leads to ‰Sii …x†Š ˆ

207

…136†

results with ‰S…x†Š ˆ ‰Si …x†Š and r ˆ ri . With ‰Se …x†Š ˆ ‰S…x†Š ‡ O…w†, the limit of Eq. (135) equals   ÿ
 ‰S…x†Š ÿ ‰E1 Š ‰E0 Šÿ1  ‰S…x†Š ÿ ‰E1 ŠT  r‰S…x†Š;r ÿ‰E2 Š ‡ x2 ‰M 0 Š ˆ 0:

…137†

The derivative with respect to r can be replaced by that with respect to x using the relationship of the dynamicsti€ness matrices of similar boundaries (Eq. (74) for the two-dimensional case). This yields   ÿ
 ‰S…x†Š ‡ ‰E1 Š ‰E0 Šÿ1  ‰S…x†Š ‡ ‰E1 ŠT  x‰S…x†Š;x ÿ‰E2 Š ‡ x2 ‰M 0 Š ˆ 0:

…138†

In an actual application a speci®c boundary is addressed which ®xes r. The dynamic-sti€ness matrix thus becomes a function of x only. The partial derivative ‰S…r; x†Š;x in Eq. (74) is replaced by ‰S…x†Š;x . The dynamic-sti€ness matrix for the bounded medium ‰S b …x†Š (superscript b for bounded) is expressed as   ÿ b
‰S …x†Š ÿ ‰E1 Š ‰E0 Šÿ1 ‰S b …x†Š ÿ ‰E1 ŠT   ‡ x S b …x† ;x ÿ‰E2 Š ‡ x2 ‰M 0 Š ˆ 0 …139† and for the unbounded medium ‰S 1 …x†Š (superscript 1 for unbounded) as   ÿ 1
‰S …x†Š ‡ ‰E1 Š ‰E0 Šÿ1 ‰S 1 …x†Š ‡ ‰E1 ŠT ÿ x‰S 1 …x†Š;x ÿ‰E2 Š ‡ x2 ‰M 0 Š ˆ 0:

…140†

208

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

This represents the scaled boundary ®nite-element equation in dynamic sti€ness formulated in the frequency domain. It is a system of non-linear ordinary di€erential equations of ®rst order with the frequency x as the independent variable. Turning to the relationship for the nodal loads of the medium, Eq. (126) is reformulated as  fPe g  RFe ˆ …  ‰S…x†Š ‡ ‰See …x†Š†‰Sie …x†Šÿ1    …141†  fPi g  RFi : Substituting Eqs. (130d) and (131) with Eqs. (133a) and (133b) and using Eq. (107b) yields …  ‰S…x†Š ‡ ‰See …x†Š†‰Sie …x†Šÿ1 ÿ
ˆ ÿ‰IŠ ÿ w  ‰Se …x†Š ‡ ‰E1 Š ‰E0 Šÿ1 ‡ O…w2 †:

 n‰S…x; n†Š;n ÿ‰E2 Š ‡ x2 n2 ‰M 0 Š ˆ 0: …142†

Substituting Eq. (142) in Eq. (141) results in fPi g ‡ fPe g fRFe g ÿ fRFi g ÿ  ‡  ‰Se …x†Š w w  
‡ ‰E1 Š ‰E0 Šÿ1 fPi g  fRFi g ˆ O…w†:

…143†

…149†

Eq. (147) for nr is formulated as

 nfRF …n†g;n ‡…‰S…x; n†Š ÿ ‰E1 Š†‰E0 Šÿ1 RF …n† ˆ fF …n†g: …144†

results with fF g ˆ fF b g ‡ fF t g:

…148†

Addressing the nodal loads, fF b g (Eq. (111)) is multiplied by n2 and fF t g by n (Eq. (116)), resulting in (Eq. (145)) fF …n†g ˆ n2 fF b g ‡ nfF t g:

From Eqs. (112) and (115) and using Eq. (119) fPi g ‡ fPe g ˆ fF g ‡ O…w† w

12b). To determine the response in the interior of the medium, the derivation can be performed for a similar boundary characterized by nr with its corresponding cell. This corresponds to multiplying the coordinates of the original boundary by the scaling factor n. All equations still apply, but formulated at the boundary characterized by nr. The determinant jJ j in Eq. (89) is multiplied by n2 . ‰B1 Š and ‰B2 Š in Eqs. (97a) and (97b) are divided by n. For this two-dimensional case, ‰E0 Š, ‰E1 Š and ‰E2 Š (Eqs. (100a)±(100c)) are independent of n, and ‰M 0 Š (Eq. (105)) is multiplied by n2 . Thus, Eq. (137) for nr is formulated as   ÿ
 ‰S…x; n†Š ÿ ‰E1 Š ‰E0 Šÿ1  ‰S…x; n†Š ÿ ‰E1 ŠT

…145†

The limit of w ! 0 can now be performed. With w de®ned in Eq. (81)  F  F  F  F Re ÿ Ri R ÿ Ri lim ˆ lim ri e re !ri w!0 w re ÿ ri  F ˆ r R ;r …146†  results with fRF g ˆ RFi and r ˆ ri . With ‰Se …x†Š ˆ ‰S…x†Š ‡ O…w† and fPi g ˆ O…w† (Eqs. (119), (112) and (115)), the limit of Eq. (143) yields   r RF ;r ‡…‰S…x†Š ÿ ‰E1 Š†‰E0 Šÿ1 RF ˆ fF g: …147† After determining ‰S…x†Š from Eq. (138), this is a system of linear ordinary di€erential equations of the ®rst order for the amplitudes of the nodal load of the medium fRF g. The boundary condition is fRF …r ˆ 0†g ˆ 0 for a bounded medium and fRF …r ! 1†g ˆ 0 for an unbounded medium. In a practical application, the dynamic-sti€ness matrix ‰S…x†Š on the boundary only is required, which can be calculated from Eq. (138). To derive this equation, it is sucient to introduce a ®nite-element cell adjacent to the boundary characterized by r only (Figs. 11b and

…150†

The scaled boundary ®nite-element equation in displacement is derived starting from that in dynamic sti€ness (Eq. (148)). The ®rst row of Eq. (123) formulated at the boundary characterized by nr is addressed ‰Sii …x; n†Šfui …n†g ‡ ‰Sie …x; n†Šfue …n†g ÿ fPi …n†g  ˆ ‰Si …x; n†Šfui …n†g  RFi …n† :

…151†

The dimensionless cell width w (Eq. (81)) is independent of n. The same applies for the two-dimensional case to ‰K 1 Š in Eqs. (130a)±(130d) as ‰E0 Š, ‰E1 Š are independent of n (Eq. (101)). From Eqs. (130a)±(130d) 1 0 ‰E Š ‡ ‰Kii1 Š ‡ O…w†; w 1 ‰Sie …x; n†Š ˆ ÿ ‰E0 Š ‡ ‰Kie1 Š ‡ O…w† w ‰Sii …x; n†Š ˆ

…152a† …152b†

follows. Substituting Eqs. (152a) and Eq. (152b) and fue …n†g ˆ fui …n†g ‡ w

fue …n†g ÿ fui …n†g w

…153†

in the ®rst two terms of Eq. (151) yields ‰Sii …x; n†Šfui …n†g ‡ ‰Sie …x; n†Šfue …n†g   1 ˆ …‰Kii1 Š ‡ ‰Kie1 Š†fui …n†g ‡ ÿ ‰E0 Š ‡ ‰Kie1 Š w w 

fue …n†g ÿ fui …n†g ‡ O…w†: w

…154†

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

For the limit of w ! 0 lim

w!0

homogeneous second-order ordinary di€erential equations for the displacement amplitude fu…n†g with the dimensionless radial coordinate n as the independent variable.

fue …n†g ÿ fui …n†g fue …n†g ÿ fui …n†g ˆ lim ri re !ri w re ÿ ri fue …n†g ÿ fui …n†g ˆ lim nri re !ri nre ÿ nri ˆ nrfu…n†g;nr ˆ nfu…n†g;n

…155†

applies. Performing the limit of w ! 0 of Eq. (151) with Eqs. (154) and (155) using Eq. (107b) leads to ‰E0 Šnfu…n†g;n ˆ



   ‰S…x; n†Š ÿ ‰E1 ŠT fu…n†g  RF …n† : …156†

Note that fPi …n†g vanishes, as it is proportional to w (Eqs. (119), (112) and (115)). ‰S…x; n†Š ˆ ‰Si …x; n†Š, fu…n†g ˆ fui …n†g and fRF …n†g ˆ fRFi …n†g apply. Eq. (148) is multiplied by fu…n†g and rearranged ÿ

 
 ‰S…x; n†Š ÿ ‰E1 Š ‰E0 Šÿ1  ‰S…x; n†Š ÿ ‰E1 ŠT   n‰S…x; n†Š;n ÿ ‰E2 Š ‡ x2 n2 ‰M 0 Š fu…n†g ˆ 0:

…157†

Substituting Eq. (156) in Eq. (157) yields ÿ

ÿ
 ‰S…x; n†Š ÿ ‰E1 Š ‰E0 Šÿ1 ‰E0 Šnfu…n†g;n  F
 R …n†  n‰S…x; n†Š;n fu…n†g ÿ ‰E2 Šfu…n†g ‡ x2 n2 ‰M 0 Šfu…n†g ˆ 0:

…158†

Substituting Eq. (150) leads to  ÿ‰E1 Šnfu…n†g;n n RF …n† ;n ‡fF …n†g  n…‰S…x; n†Šfu…n†g†;n ÿ‰E2 Šfu…n†g ‡ x2 n2 ‰M 0 Š  fu…n†g ˆ 0:

…159†

Di€erentiating Eq. (156) with respect to n yields  …‰S…x; n†Šfu…n†g†;n  RF …n† ;n   ˆ ‰E0 Šnfu…n†g;nn ‡ ‰E0 Š ‡ ‰E1 ŠT fu…n†g;n :

…160†

Substituting Eq. (160) in Eq. (159) results in the scaled boundary ®nite-element equation in displacement   ‰E0 Šn2 fu…n†g;nn ‡ ‰E0 Š ÿ ‰E1 Š ‡ ‰E1 ŠT nfu…n†g;n ÿ ‰E2 Šfu…n†g ‡ x2 n2 ‰M 0 Šfu…n†g ‡ fF …n†g ˆ 0:

209

…161†

Note that this equation applies to both bounded and unbounded media. Eq. (161) is a system of linear non-

6. Conclusions A novel boundary-element method based on ®nite elements to analyse bounded and unbounded media is developed. Two derivations for this scaled boundary ®nite-element method are presented. In the ®rst, the scaled-boundary-transformationbased derivation, a scaling centre is selected. The distance from this origin to a point is represented by the radial coordinate n. Discretizing the boundary with surface ®nite elements determines the other two coordinates g, f. Applying the scaled boundary transformation to the geometry, the governing partial di€erential equations in the Cartesian coordinates are formulated in the local coordinates n, g, f. The boundary conditions are conveniently formulated in n, g, f. Using the weighted-residual technique of ®nite elements in the two circumferential directions g, f parallel to the boundary results in the scaled boundary ®nite-element equation in displacement amplitudes, a system of linear secondorder ordinary di€erential equations with the radial coordinate n as the independent variable. The scaled boundary ®nite-element equation can also be formulated in dynamic sti€ness. In the second, the mechanically based derivation, a similarity centre which corresponds to the scaling centre is selected. A similar ®ctitious boundary is constructed by scaling the boundary. A ®nite-element cell is introduced between the two boundaries. Assembling the ®nite-element cell and the medium yields a relationship linking the dynamic-sti€ness matrices at the boundaries of the medium via the dynamic-sti€ness matrix of the ®nite-element cell. Another relationship of the dynamicsti€ness matrices at the boundaries follows from similarity. Performing the limit of the cell width towards zero analytically yields the scaled boundary ®niteelement equation in dynamic sti€ness, which can also be expressed in displacement. Only the boundary is discretized with curved surface ®nite elements. When modelling a semi-in®nite medium, the free surface extending to in®nity is not discretized. The same applies for straight crack faces in fracture mechanics. The analytical and numerical solution procedures for the scaled boundary ®nite-element equations in displacement and dynamic sti€ness for bounded and unbounded media in the frequency and time domains and in statics are discussed in the accompanying paper [9].

210

J.P. Wolf, C. Song / Computers and Structures 78 (2000) 191±210

References [1] Beskos DE. Boundary element methods in dynamic analysis. Appl Mech Rev ASME 1987;40:1±23. [2] Beskos DE. Boundary element methods in dynamic analysis: Part II (1986±1996). Appl Mech Rev ASME 1997; 50:149±97. [3] Song Ch, Wolf JP. Consistent in®nitesimal ®nite-element cell method: three-dimensional vector wave equation. Int J Numer Meth Engng 1996;39:2189±208. [4] Song Ch, Wolf JP. The scaled boundary ®nite-element method ± alias consistent in®nitesimal ®nite-element cell method ± for elastodynamics. Comput Meth Appl Mech Engng 1997;147:329±55.

[5] Wolf JP, Song Ch. Consistent in®nitesimal ®nite-element cell method in frequency domain. Earthquake Engng Struct Dynam 1996;25:1307±27. [6] Wolf JP, Song Ch. Unit-impulse response of unbounded medium by scaled boundary ®nite-element method. Comput Meth Appl Mech Engng 1998;159:355±67. [7] Song Ch, Wolf JP. The scaled boundary ®nite-element method: analytical solution in frequency domain. Comput Meth Appl Mech Engng 1998;164:249±64. [8] Wolf JP, Song Ch. Finite-Element Modelling of Unbounded Media. Wiley, Chichester, 1996. [9] Song Ch, Wolf JP. The scaled boundary ®nite-element method ± a primer: solution procedures. Comput Struct 2000;78:211±25.