close

Вход

Забыли?

вход по аккаунту

?

BEAM STIFFNESS MATRIX BASED ON THE ELASTICITY EQUATIONS

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL.
40, 89–109 (1997)
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN
BEM BY EXACT TRANSFORMATION FOR 2-D
ANISOTROPIC ELASTICITY
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
Department of Mechanical and Aerospace Engineering, Carleton University, 3135 Mackenzie Bldg,
1125 Colonel By Drive, Ottawa, Canada K1S 5B6
SUMMARY
In the Boundary Element Method (BEM) based on the direct formulation, body-force eects manifest themselves as an additional volume integral term in the Boundary Integral Equation (BIE). The numerical solution
of the integral equation with this term destroys the notion of the BEM as a truly boundary solution method.
This paper discusses the treatment of this volume integral for two-dimensional anisotropic elasticity with bodyforces present. The analytical basis for transforming this integral exactly into boundary ones is presented for
geometrically convex regions. This restores the application of the BEM to such problems as a truly boundary
solution technique. Numerical examples are presented to demonstrate the veracity of the transformation and
implementation.
KEY WORDS:
BEM; numerical fracture mechanics
1. INTRODUCTION
The Boundary Element Method (BEM), also known as the Boundary Integral Equation (BIE)
method, is an ecient numerical tool for engineering stress analysis. In elastostatics, the analytical
basis of the so-called direct BIE formulation, which the present work is focused on, is the existence
of the fundamental solution to the governing dierential equations, and the reciprocal work theorem.
From these, and following the usual limiting process, Somigliana’s identity and the BIE, which
relate the boundary displacements to the boundary tractions, may be obtained.
The fact that integrals present in the BIE and Somigliana’s identity are surface integrals establishes the BEM as a boundary solution technique as distinct from domain solution techniques,
such as the nite element method. The implication here is that only the surface of the domain
needs to be modelled or discretized, instead of the entire domain as in FEM. However, at the
basic level of formulation, this holds true only if body forces and temperature changes are absent.
When body forces and thermal eects are considered, the basic form of Somigliana’s identity and
the BIE will not only include surface integrals, but volume integrals as well. The direct numerical
evaluation of these integrals necessitates the discretization of the entire domain into ‘cells’, so that
integration over each of these ‘cells’ may be carried out. Although these integration ‘cells’ are not
‘nite elements’, since they do not contribute to the coecients corresponding to the unknowns
in the nal equation matrix, the notion of BEM as a boundary solution technique is destroyed in
such cases. Thus the solution to the ‘volume integral problem’ is of fundamental importance to
the development of BEM.
CCC 0029–5981/97/010089–21
? 1997 by John Wiley & Sons, Ltd.
Received 27 July 1995
Revised 8 February 1996
90
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
There are evidently two dierent general schemes to resolve the volume integral problem in the
BIE formulation. The rst is to attempt to perform the volume integration without any explicit
discretization of the solution domain, such as the Monte Carlo Method, developed by Gipson
and Camp,1 or the domain Fanning method, proposed by Camp and Gipson.2 Although a very
imaginative idea, the Monte Carlo method is not ecient because a large number of calculation
points must be randomly chosen, and the accuracy of the solution is not always assured. The
domain Fanning method takes advantage of the properties of the fundamental solution and improves
the eciency of evaluating the volume integral through an implicit domain discretization scheme.
However, for domains with certain boundary shapes, diculties will arise.
The second and more attractive procedure is to further transform the volume integrals to surface or boundary integrals. The exact transformation of volume integrals to boundary integrals,
by Rizzo and Shippy,3 the dual reciprocity method, by Nardini and Brebbia,4 the multiple reciprocity method, by Nowak,5 provide examples of this general scheme. These methods have met
with tremendous success in isotropic elasticity. However, the developments for anisotropic elasticity have, hitherto, not been that successful. This is primarily due to the signicantly more
complicated form of the fundamental solutions in anisotropic cases. A limited amount of work
in this regard has been reported recently, using the ‘particular integral approach’. It was rst
discussed by Watson,6 and later by Banerjee and Buttereld,7 for isotropic elasticity. This approach has recently been further extended to anisotropic elasticity by Deb and Banerjee,8 and
Deb et al.9 Yet another approach, developed by Ao17 very recently, is to reformulate the elasticity equations into a BIE which is in terms of unknown scalar potentials to be solved for. The
method, as presented in Reference 17 and Ao and Shippy,18 applies, however, to only orthotropic
materials.
The particular integral approach mentioned above has been successful in solving problems where
the stress distribution is smooth in anisotropic media. However, it lacks generality to treat problems
involving stress singularities such as those encountered in cracks. The motivation of the present
work is to overcome these diculties, and an attempt to develop a method which transforms
exactly the volume integrals to boundary ones in 2-D anisotropic elasticity is presented here.
Although the approach discussed is restricted to simply connected convex regions, it nevertheless
represents, to the authors’ knowledge, the rst reported successful eort in the use of this scheme
for treating volume integrals in the BIE method with body force eects for generally anisotropic
domains.
2. THE BIE FOR PLANE ANISOTROPIC ELASTICITY
For a plane stress state of a 2-D homogeneous anisotropic elastic region, the generalized Hooke’s
law can be written as (see Reference 12)
 



11 12 16  11 
 ”11 
”22 =  12 22 26  22
(1)




12
16 26 66
12
where (”11 ; ”22 ; 12 ) and (11 ; 22 ; 12 ) are strains and stresses, respectively, and constants ij are
the elastic compliances of the materials. In the case of plane strain conditions, the coecients ij
in equation (1) are replaced by ij , where
ij = ij −
i3 j3
33
(2)
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
91
According to Lekhnitskii, the governing dierential equation for a two-dimensional problem of
anisotropic elasticity can be obtained in terms of Airy’s stress function F as follows:
22
@4 F
@4 F
@4 F
@4 F
@4 F
−
2
+
(2
+
)
−
2
+
=0
26
12
66
16
11
@41
@21 @22
@42
@31 @2
@1 @32
(3)
The corresponding characteristic equation of (3) is
22 − 226 + (212 + 66 )2 − 216 3 + 11 4 = 0
(4)
It can be shown that for anisotropic materials, equation (4) has only complex roots, and they
are distinct (see Reference 12). In addition, because the coecients of equation (4) are all real
constants, the roots are in complex conjugate pairs. The four roots are denoted by
1 ; 1 ; 2 ; 2
In terms of generalized complex variables, i.e.
z1 = (1 − y1 ) + 1 (2 − y2 )
z2 = (1 − y1 ) + 2 (2 − y2 )
(5)
where (1 ; 2 ) is the eld point, and (y1 ; y2 ) is the source point, the fundamental solution for an
anisotropic 2-D elastic body can be obtained as (see Reference 11)
Uji (; y) = 2 Re[ri1 Aj1 log z1 + ri2 Aj2 log z2 ]
(6)
where Re [·] is an operator which takes the real part of a complex function. The constants rij
are
r11 = 11 12 + 12 − 16 1
r12 = 11 22 + 12 − 16 2
r21 = 12 1 +
22
− 26
1
r22 = 12 2 +
22
− 26
2
(7)
The explicit equations for constants Aij can be found in Tan et al.10 The corresponding traction
fundamental solutions of equation (6) can be obtained as
Tj1 = 2n1 Re[12 Aj1 =z1 + 22 Aj2 =z2 ]
− 2n2 Re [1 Aj1 =z1 + 2 Aj2 =z2 ]
Tj2 = − 2n1 Re [1 Aj1 =z1 + 2 Aj2 =z2 ]
(8)
+ 2n2 Re [Aj1 =z1 + Aj2 =z2 ]
where nj are unit outward normal components at a boundary point (x1 ; x2 ).
The direct formulation of BEM for anisotropic materials may be developed by following the
same steps as in the isotropic cases. The resulting boundary integral equation provides a direct
92
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
integral relation between boundary tractions, ti , and boundary displacements, ui , and is given as
(see Reference 11)
I
Cji ui (y) +
ui (x)Tji (x; y) dSx
@
I
=
ti (x)Uji (x; y) dSx
Z@
+
Xi ()Uji (; y) d
;
x; y ∈ @
∈ (9)
where Xi are body forces, and @
and represent the boundary and the domain, respectively. In
elastostatics, the body forces can usually be expressed as the gradient of a scalar potential , i.e.
Xi = ; i
; ii = k0
(10)
where k0 is a constant.
The distinct advantage of the BEM as a numerical tool over domain solution methods is that
only boundary discretization of the problem domain is needed. However, with the volume integral
(VI), i.e.
Z
(VI)j =
Xi ()Uji (; y) d
; y ∈ @
; ∈ (11)
due to the body forces Xi , appearing in the boundary integral equation, the notion of the method
being a boundary solution technique is now violated. Because closed-form solutions to the above
VI cannot be obtained generally, numerical evaluation is usually required. But the direct numerical evaluation of equation (11) needs discretization of the problem domain. Here, an attempt to
transform the VI into surface or boundary integrals is made so as to restore equation (9) as a
boundary-only integral equation.
3. EXACT TRANSFORMATION FOR CONVEX REGIONS
This paper discusses the method for transforming the VI to boundary integrals for convex regions.
The term convex region applies to a region whose tangent at every boundary point lies outside the
region. For example, regions shown in Figure 1 are convex regions, while those shown in Figure
2 are not.
With reference to the co-ordinate system shown in Figure 1, one can always nd a highest
point A (with maximum 2 ), and a lowest point C (with minimum 2 ) along the boundary of a
convex region. In some cases, these points may lie along straight lines of constant 2 , as shown
in Example II, Figure 1. The points A and C divide the entire boundary of a convex region into
˙
˙
two parts, namely the arcs ABC and CDA . For a case as in Example II, Figure 1, one can always
specify any point along line A as the point A, and any point along line C as the point C.
In dealing with the VI in equation (11), one has to deal with the logarithmic function of a
generalized complex variable, i.e.
z = (1 − y1 ) + (2 − y2 )
(12)
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
93
Figure 1. Typical convex regions
where = a + ib represents the characteristic roots, and b never vanishes for anisotropic materials
since equation (4) has only complex roots. A two-step mapping of equation (12) follows. First,
the (1 ; 2 ) plane is mapped to an intermediate (1 0 ; 2 0 ) plane through
01 = 1 − y1
02 = 2 − y2
(13)
Then, the intermediate (01 ; 02 ) plane is mapped to the nal (1 ; 2 ) plane through
1 = 01 + a02
2 = b02
(14)
In the (1 ; 2 ) plane, the variable z can be written as standard complex variable,
z = 1 + i2
(15)
Through equations (13) and (14), the original domains in Figure 1 are mapped into the new
domain , as shown in Figure 3.
94
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
Figure 2. Typical non-convex regions
From equation (13), one may note that the origin of the (01 ; 02 ) system is always on the
boundary and is coincident with the source point. Furthermore, it is evident that for source points
˙
along arc ABC , the positive 01 axis cuts across the domain, while for source points along arc
˙
CDA , the negative 01 axis cuts across the domain. The mapping depicted by equation (14) has
the distinctive feature of keeping the points along the 01 axis unchanged. It can be easily seen
that for points 02 = 0, one always has 1 = 01 . In addition, it can be noted that the image of
point A remains the highest point if b ¿ 0, and it becomes the lowest point if b ¡ 0. Conversely,
the image of point C in the (1 ; 2 ) plane remains the lowest point if b ¿ 0, and it becomes
the highest point if b ¡ 0. Hence, the images of points A and C remain as the dividing points
for the boundary @
of the mapped region in the (1 ; 2 ) plane. Also, the image of the arc
˙
˙
ABC remains on the left side of the region, and the image of the arc CDA remains on the right
side of the region.
Thus, one can conclude that if from a source point, the negative 1 0 axis cuts across the region
, then the negative 1 axis will cut across the corresponding region ; and if from a source
point, the positive 01 axis cuts across the region , then the positive 1 axis will cut across the
corresponding region . Similarly, if the 01 axis does not cut across the region , then the 1
˙
axis does not cut across the corresponding region . Thus, for any point along the arc ABC , the
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
95
Figure 3. The images of original regions in the 1 2 plane
˙
positive 1 axis cuts across the region ; and for any point along arc CDA , the negative 1 axis
cuts across the region .
By denition, the complex function log(z) is
log (z) = log | z | + i arg (z)
(16)
In computing the fundamental solution, the principal value of arg (z) is employed, i.e.
− ¡ arg (z)6
(17)
It is clear from equation (17) that along the negative 1 axis, arg (z) is not continuous. Thus,
˙
for any source point along ABC , the function log (z) is analytic in the problem domain; but for
˙
any source point along CDA , the function log (z) is not analytic in the problem domain because
the negative 1 axis cuts across the region. However, this diculty can be overcome by simply
˙
redening the principal value of arg(z) for the source points along CDA of a convex region as
0 ¡ arg (z)62
(18)
96
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
In equation (18), arg(z) is now discontinuous along the positive axis, but continuous along
the negative 1 axis, and log (z) becomes analytical in the problem region for source points along
˙
CDA .
With the understanding that log (z) is made to be analytic in the problem domain for every source
point along the boundary by redening the limits of its argument, a method for transforming the
volume integral in equation (11) to boundary integrals may now be formulated.
First, the following notation is introduced
1 1
[km ] =
(19)
1 2
so that with the application of the summation convention, equation (5) can be simplied as
zm = km (k − yk );
m; k = 1; 2
(20)
Next, two functions G1 (; y) and G2 (; y) are introduced such that they satisfy the following
dierential equations respectively,
Gm; jj (; y) =
=
1
zm
1
km (k − yk )
;
m; k = 1; 2
(21)
where dierentiation is taken to be with respect to the eld point . To nd a solution to equation
(21), it is important to note that the second derivative of a real function log () with respect to
variable is 1=. This suggests that the following are particular solutions of equation (21),
G1 = C1 z1 log (z1 )
G2 = C2 z2 log (z2 )
(22)
where Cm (m = 1; 2) are constants to be determined. Applying the Laplacian operator to equation
(22), and noting that log (z) is now analytic in the problem domain, gives
G1; jj = C1 (k1 kj log (z1 ) + k1 kj ); j
1
= C1 j1 j1 ·
z1
1
G2; jj = C2 j 2 j 2 ·
z2
(23)
where ij is the Kronecker delta function. By substituting equation (23) into equation (21), the
constants Cm can be obtained as
C1 =
1
;
j1 j1
C2 =
1
j 2 j 2
(24)
Thus, the functions Gm ; m = 1; 2, can be nally determined as
G1 =
1
1
z1 log (z1 ) = 2
z log (z1 )
2 1
j1 j1
11 + 21
1
1
G2 =
z2 log (z2 ) = 2
z log (z2 )
2 2
j2 j2
12 + 22
(25)
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
97
From functions Gm and by using a similar structure for the fundamental solution, two important
functions can be constructed, i.e.
Hj (; y) = 2 Re[ri1 i1 Aj1 G1 (; y) + ri2 i2 Aj 2 G2 (; y)];
j = 1; 2
(26)
Applying the Laplacian operator to functions Hj , and noting equation (21), one can obtain
Hj; kk = 2 Re[ri1 i1 Aj1 G1; kk (; y) + ri2 i2 Aj 2 G2; kk (; y)]
1
1
= 2 Re ri1 i1 Aj1 + ri2 i2 Aj 2
z1
z2
(27)
The relation between the functions Hj ; j = 1; 2, and the fundamental solution Uji (x; y); i; j =
1; 2, may now be established. Applying the divergence operator to the fundamental solution gives
Uji; i (; y) = 2 Re[ri1 Aj1 log (z1 ) + ri2 Aj 2 log (z2 )]; i
1
1
= 2 Re ri1 Aj1 k1 ki + ri2 Aj 2 k2 ki
z1
z2
1
1
= 2 Re r i1 Aj1 i1 + ri2 Aj 2 i2
z1
z2
(28)
Comparing equation (27) to equation (28), it is evident that
Hj; ii (; y) = Uji; i (; y)
(29)
Two vector functions E1 = (E11 ; E21 )T , and E2 = (E12 ; E22 )T are introduced next. They are
related to functions Gm (m = 1; 2) by
Ei1; i = G1 (; y) =
1
z1 log (z1 )
j1 j1
Ei2; i = G2 (; y) =
1
z2 log (z2 )
j 2 j 2
(30)
That is, the function Gm (m = 1; 2) is the divergence of vector Em . To nd a particular solution
to the above dierential equations, it is important to note that the derivative of a real function
2 log () with respect to is 2 log () + . This suggests that the vectors Em may have the
following structure,
Ei1 = Di1 z 21 log (z1 ) + Fi1 z 21
Ei2 = Di2 z22 log (z2 ) + Fi2 z22
(31)
where constants Dim and Fim ; i; m = 1; 2, are to be determined. By applying the divergence
operator to the rst of equation (31), i.e.
Ei1; i = Di1 (z12 log (z1 )); i + Fi1 (z12 ); i
(32)
(z12 ); i = 2z1 k1 ki = 2z1 i1
(33)
and noting that
(z12 log (z1 )); i = 2z1 k1 ki log (z1 ) + z1 k1 ki
= 2i1 z1 log (z1 ) + i1 z1
(34)
98
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
one can reduce equation (32) to
Ei1; i = 2Di1 i1 z1 log (z1 ) + Di1 i1 z1 + 2Fi1 z1 i1
By comparing equation (35) to equation (30), one obtains
1
2Di1 i1 =
j1 j1
(35)
(36)
Di1 i1 + 2Fi1 i1 = 0
(37)
From equations (36) and (37), the constants Di1 and Fi1 can be chosen as
1
Di1 =
4i1 j1 j1
1
1
Fi1 = − Di1 = −
2
8i1 j1 j1
(38)
(39)
Similarly, Di2 and Fi2 are determined to be
1
4i2 j 2 j 2
1
=−
8i2 j 2 j 2
Di2 =
(40)
Fi2
(41)
Thus the vectors E1 (; y) and E2 (; y) are found as
1
Ei1 (; y) =
z 2 (log z1 − 12 )
4i1 j1 j1 1
1
z 2 (log z2 − 12 )
Ei2 (; y) =
4i2 j 2 j 2 2
(42)
(43)
By employing the vector functions E1 (; y) and E2 (; y), two other vector functions W1 =
(W11 ; W12 ) and W2 = (W21 ; W22 ) can be constructed by using a structure similar to that of the
fundamental solution, i.e.
Wji = 2 Re[rk1 Aj1 k1 Ei1 + rk2 Aj 2 k2 Ei 2 ];
i; j; k = 1; 2
(44)
Applying the divergence operator to the vectors Wj yields
Wji; i = 2 Re[rk1 Aj1 k1 Ei1; i + rk2 Aj 2 k2 Ei2;i ]
= 2 Re[rk1 Aj1 k1 G1 + rk2 Aj 2 k2 G2 ]
(45)
Comparing equation (26) to equation (45), one can immediately recognize that
Wji; i (; y) = Hj (; y)
(46)
With the help of functions Hj ; j = 1; 2 and vectors Wi ; i = 1; 2, it is now possible to transform
the volume integral equation (11) to boundary ones. Substituting equation (10) into equation (11)
and applying Green’s theorem gives
Z
(VI)j =
I
; i Uji (; y) d
=
@
Z
Uji (x; y)ni dSx −
Uji; i (; y) d
;
∈ ; y; x ∈ @
(47)
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
99
where n = (n1 ; n2 ) is the outward normal of the boundary at point (x1 ; x2 ). From equation (29),
the divergence of the fundamental solution Uji can be expressed in terms of the functions Hj .
Thus, by applying Green’s theorem consecutively, one has
I
Z
(VI)j =
Uji (x; y)ni dSx −
Hj; ii (; y) d
I @
I
Z
=
Uji (x; y)ni dSx −
Hj; i (x; y)ni dSx +
; i Hj; i (; y) d
I @
I @
=
Uji (x; y)ni dSx −
Hj; i (x; y)ni dSx
@
@
I
Z
+
; i Hj (x; y)ni dSx −
; ii Hj (; y) d
(48)
@
From equation (10) and (46), the last volume integral in the above equation can be further reduced
to
Z
Z
; ii Hj (; y) d
=
k0 Wji; i (; y) d
I
=
k0 Wji (x; y) ni dSx
(49)
@
Thus, the volume integral in equation (11) nally can be expressed as boundary integrals, i.e.
I
I
(VI)j =
Uji (x; y)ni dSx −
Hj; i (x; y)ni dSx
@
@
I
I
+
; i Hj (x; y)ni dSx −
k0 Wji (x; y)ni dSx
(50)
@
@
Equation (50) successfully transforms the volume integral associated with body forces in the
BIE method into boundary ones. Therefore, the basic integral equation of the BEM for the 2-D
anisotropic elastostatic problem can be written as
I
Cji ui (y) +
ui (x)Tji (x; y) dSx
I @
=
ti (x)Uji (x; y) dSx
@
I
I
+
Uji (x; y)ni dSx −
Hj; i (x; y)ni dSx
@
@
I
I
+
; i Hj (x; y)ni dSx −
(51)
k0 Wji (x; y)ni dSx
@
@
Equation (51) is now a boundary-only integral equation, which, at least for convex regions, overcomes the volume integral problem of the BEM for 2-D anisotropic elastostatics. The numerical
solution of equation (51) may be obtained using, for example, the conventional quadratic isoparametric element formulation, the details for which may be found in e.g. Reference 6. The evaluation
of the transformed surface integrals containing functions H and W due to the body force term
presents no diculties as they are not singular.
The method developed above can also be applied to the case of some non-convex regions. For
instance, in the case of Example III, Figure 2, the principal value of the argument of a complex
variable can be redened to be within the limits
3
(52)
− ¡ arg (z)6
2
2
100
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
for all the source points along the inner semi-circular boundary. In a similar manner, log (z) can
be made analytic in the problem domain for all source points along the whole boundary. It should
perhaps also be remarked that, in general, a non-convex region can be judiciously divided into
several convex sub-regions, or sub-regions like the case of Example III, Figure 2, on each of
which the method described herein may be applied. This concept of sub-regioning is, of course,
quite common in BEM analysis.
4. NUMERICAL EXAMPLES
Three test problems are treated to demonstrate the validity of the exact transformation procedure
developed in this paper. In all these examples, the properties of materials are described in terms
of the engineering elastic constants, which are related to the elastic compliances ij as follows
(see, Reference 12),
1
E1
12
21
=−
= −
E1
E2
1
=
E2
12;1
1;12
=
=
E1
G12
1
=
G12
12; 2
2; 12
=
=
E2
G12
11 =
12
22
16
66
26
(53)
where Ek is the Young’s modulus in the k direction, G12 is the shear modulus in the (1 ; 2 )
plane, and jk is the Poisson’s ratio which is dened as the compressive strain in the k direction
due to unit extensional strain in the j direction. Constants ij; k and i; jk are referred to as the
coecients of mutual inuence of the rst and second kind, respectively.
The quadratic isoparametric element formulation was employed in the present study. In the
rst test problem, the volume integral associated with the body force was also evaluated for
each boundary node by means of direct two-dimensional numerical integration and the values are
compared with those obtained from evaluating the transformed boundary integrals. The second
test problem involves the determination of the stress distribution in a rotating, plane, circular,
orthotropic disc for which the exact analytical solution exists. Finally, a rotating, plane, circular,
anisotropic disc with a central crack is analysed. The 8-point Gauss quadrature was used throughout
for evaluating the boundary integrals in equation (51)
Problem 1
A square anisotropic region of 2 m × 2 m is considered with the BEM mesh discretization as
shown in Figure 4. The material constants are chosen to be as follows:
E1
(GPa)
E2
(GPa)
G12
(GPa)
12
12;1
12;2
345
516
173
0·131
0·0
0·0
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
101
Figure 4. Problem 1—A square domain and the BEM mesh
Figure 5. Body force potential for Problem 1
The body forces are taken to be
X1 = C [(e−2 − e2 ) cos(1 ) + (e−1 − e1 ) sin(2 )]
X2 = C [−(e−2 − e2 ) sin(1 ) − (e−1 − e1 ) cos (2 )]
(54)
with the corresponding potential
= C [(e−2 − e2 ) sin(1 ) − (e−1 − e1 ) sin(2 )]
(55)
which is shown graphically in Figure 5. C is a constant, and is chosen to be 106 in the numerical
calculations, assigning both X1 and X2 to have units of N (Newton). Although this may not be
a physically encountered situation, the choice of this potential is to demonstrate the mathematical
soundness of the approach developed.
In the analysis, the boundary of the square is divided into 8 elements with a total of 16 nodes
(Figure 4). Table I shows the results of the body force volume integral (VI) as obtained by the
Direct Integration Method (DIM) of evaluating equation (11), and by the ‘Exact Transformation
Method’ (ETM) of evaluating equation (50). In the former method, the solution domain was
102
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
Table I. Values of the volume integral for the square region—Problem 1
(VI)1 (mm)
Node
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
(VI)2 (mm)
ETM
DIM
ETM
DIM
0·624095E − 02
0·762926E − 02
− 0·149011E − 07
− 0·762925E − 02
− 0·624094E − 02
− 0·567059E − 02
− 0·383093E − 02
− 0·567054E − 02
− 0·624099E − 02
− 0·762930E − 02
0·163917E − 07
0·762930E − 02
0·624091E − 02
0·567060E − 02
0·383092E − 02
0·567057E − 02
0·624251E − 02
0·733150E − 02
0·174622E − 09
− 0·733150E − 02
− 0·624251E − 02
− 0·566292E − 02
− 0·383115E − 02
− 0·566292E − 02
− 0·624251E − 02
− 0·733150E − 02
0·000000E + 00
0·733150E − 02
0·624251E − 02
0·566292E − 02
0·383115E − 02
0·566292E − 02
0·533362E − 02
0·514822E − 02
0·402040E − 02
0·514822E − 02
0·533363E − 02
0·629162E − 02
0·167638E − 07
− 0·629167E − 02
− 0·533360E − 02
− 0·514820E − 02
− 0·402040E − 02
− 0·514823E − 02
− 0·533361E − 02
− 0·629162E − 02
0·223517E − 07
0·629164E − 02
0·533267E − 02
0·515314E − 02
0·402029E − 02
0·515314E − 02
0·533267E − 02
0·645052E − 02
− 0·116415E − 09
− 0·645052E − 02
− 0·533268E − 02
− 0·515314E − 02
− 0·402029E − 02
− 0·515314E − 02
− 0·533267E − 02
− 0·645052E − 02
0·174622E − 09
0·645052E − 02
Figure 6. A circular rotating disc and the co-ordinate system—Problem 2
divided into 64 ‘cells’ and a 10 × 10 point Gauss quadrature was used per cell. As can be seen
in the table, even with the small number of boundary elements, the results for the VI from ETM
and the DIM are in very good agreement indeed.
Problem 2
The stress distribution in a thin circular orthotropic disc which is rotating about its centroidal
axis x3 is analysed next (Figure 6). The results of this BEM analysis can be compared to
those obtained from the analytical closed-form solutions by Lekhnitskii.13 For such a solid disc
with radius R and rotating at a constant angular velocity !, the stress components depend only
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
103
Figure 7. BEM mesh—Problem 2
on the distance r from the centre of the rotation and, in a polar co-ordinate system, they are
given as13
r = 12 !2 R2 (1 − A)(1 − (r=R)2 )
= 12 !2 R2 [(1 − A)(1 − (r=R)2 ) + 2A(r=R)2 ]
(56)
r = 0
A =
11 + 212 + 22
311 + 212 + 66 + 322
where the constant A is an invariant with the transformation of the co-ordinate system. For the purpose of analysis, the material used in this example is a glass–epoxy composite with the following
elastic constants:
E1
(GPa)
E2
(GPa)
G12
(GPa)
12
12;1
12;2
48·26
17·24
6·89
0·29
0·0
0·0
The body forces due to the rotation are
X1 = x1 !2
(57)
X2 = x2 !
(58)
2
Therefore, the body force potential is
= 12 !2 (x12 + x22 )
(59)
Since constant A is an invariant with the transformation of the co-ordinate system, for simplicity,
the material principal axes are also taken as x1 and x2 here. Due to the symmetry with respect to
104
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
Figure 8. Variation of radial and tangential stresses in the rotating disc with radial distance—Problem 2
Table II. Percentage errors for radial and tangential stresses—Problem
r=R
0·0000
0·0625
0·1250
0·1875
0·2500
0·3750
0·5000
% error of r
% error of r=R
% error of r
% error of 0·2058
1·7651
1·1089
1·5382
0·4369
1·9147
− 0·0156
0·2229
0·4814
0·2103
0·7569
− 0·4029
0·6843
− 0·5445
0·6250
0·7500
0·8125
0·8750
0·9375
1·0000
2·6281
0·9407
2·9434
1·6480
3·5840
0·6703
0·4282
− 0·9390
1·9364
− 2·7796
0·2496
0·9766
x1 and x2 axes, only the rst quarter of the disc need be considered with the appropriate boundary
condition. The corresponding boundary discretization is shown in Figure 7, where 20 boundary
elements with a total of 40 nodes are employed.
The BEM results of radial and tangential stresses, normalized by !2 R2 , together with the exact
solution from Lekhnitskii13 are shown in Figure 8, where it can be seen that the agreement between
the two sets of results is excellent. The corresponding percentage errors are listed in Table II.
Problem 3
Next, a thin circular anisotropic rotating disc with a central crack as shown in Figure 9 is
investigated. It is assumed that the crack can have any arbitrary orientation with respect to the
material principal axis. The origin of the co-ordinate system is taken to be at the centre of the
disc with the x1 axis aligned along the plane of the crack. The stress intensity factors KI and KII
are to be determined.
As the closed-form solution for this problem is not readily available, Bueckner’s superposition
method14 was used to obtain an independent set of results for comparison with those obtained
with the present approach. Using this method, the stress intensity factor for the crack in the
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
105
Figure 9. Circular rotating disc with a central crack—Problem 3
Figure 10. Circular non-rotating disc with normal tractions along crack faces
rotating disc is equal to that for a crack in an identical, but non-rotating disc for which the
crack faces are subjected to normal tractions equal to the circumferential stress acting on the
prospective crack plane in the uncracked rotating disc (Figure 10). This circumferential stress
distribution is readily available in the literature.13 In this manner, the original crack problem
with body forces is converted to a simple crack problem without any body forces. This problem can now be solved using any standard BEM code, the solution of which may thus be
used for verication of the results from the exact volume to boundary integral transformation
approach.
Due to the material anisotropy, the problem is not symmetric about either the x1 or x2 axis
in general. Hence the boundary of the whole domain needs to be discretized. First, the domain
is divided into two sub-regions along the plane of the crack; a typical mesh used is shown in
106
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
Figure 11. Typical BEM mesh—Problem 3
Figure 11. Traction-singular quarter-point boundary elements, see e.g. Reference 15, were employed
adjacent to the crack tips and stress intensity factors were computed directly from the computed
nodal traction coecients there, as is well established in the literature for BEM fracture mechanics
analysis.15; 16
The material used in this example is the T300 / N5208 composite which has the following elastic
constants:8
E1
(GPa)
E2
(GPa)
G12
(GPa)
12
12;1
12; 2
131
13
6·4
0·038
0·0
0·0
Several material property orientations, as measured by the angle (see Figure 9), namely, =
0◦ , 15◦ , 30◦ , 45◦ , 60◦ , 75◦ and 90◦ , as well as normalized crack length, a=R = 0·125, 0·25, 0·375,
0·5, 0·625, 0·75, 0·875 were considered. Figure 12 shows the results for the non-dimensionalized
stress intensity factor KI =Ka at the crack tips as functions of crack length, where
√
Ka = 0 a
(60)
and 0 is the tangential stress at r = 0 in the absence of the crack. The percentage deviations
between the values obtained by Bueckner’s method and the present method are generally less than
1 per cent, the greatest discrepancy being only 1·2 per cent. The corresponding results for KII =Ka as
a function of crack length a=R are shown in Figure 13 for the dierent orientations of the material
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
107
√
Figure 12. Variation of normalized stress intensity factor KI =Ka with relative crack length a=R; Ka = 0 a
axes. The results of KII for = 0 and 90◦ of the material axes are zero for any crack length,
as is to be expected. It is evident again that there is excellent agreement between the results
obtained using the BEM formulation presented in this paper and Bueckner’s superposition approach. The relative dierences for KII are however larger than the ones for KI . This is attributed to the fact that the KII are numerically much smaller in magnitude compared to KI ;
hence their values are more easily inuenced by the calculation round-o errors in numerical
analysis.
108
J. J. ZHANG, C. L. TAN AND F. F. AFAGH
√
Figure 13. Variation of normalized stress intensity factor KII =Ka with relative crack length a=R; Ka = 0 a
5. CONCLUSION
In the direct BEM for elastostatics, the presence of body force eects results in a volume integral
term in the boundary integral equations. In isotropic elasticity, this volume integral can be transformed exactly into boundary ones. This enables the BEM to maintain itself as a truly boundary
solution computational technique. Hitherto, however, the same cannot be said of the BEM for
general anisotropic elasticity involving body forces.
In this study, the exact transformation of the volume integral to boundary integrals, for two
dimensional anisotropic regions, has been presented. The validity of the analysis has been demon-
TREATMENT OF BODY-FORCE VOLUME INTEGRALS IN BEM
109
strated by three examples. Although this exact transformation has been successfully implemented
in the BEM for convex geometric regions, the study nevertheless suggests the feasibility of the
approach to treat body forces and, in principle, thermoelastic eects. Research to remove this
limitation on the geometry of the anisotropic region analysed and to extend the procedures to the
thermoelastic problem is currently in progress.
REFERENCES
1. G. S. Gipson and C. V. Camp, ‘Eective use of Monte Carlo quadrature for body force integrals occurring in integral
form of elastostatics’, Proc. 7th Int. Conf. on Boundary Elements, 1985, pp. 17–26.
2. C. V. Camp and G. S. Gipson, Boundary Element Analysis of Nonhomogeneous Biharmonic Phenomena, Springer,
Berlin, 1992.
3. F. J. Rizzo and D. J. Shippy, ‘An advanced boundary integral equation method for three dimensional thermoelasticity’,
Int. j. numer. methods eng., 11, 1753–1768 (1977).
4. D. Nardini and C. A. Brebbia, ‘A new approach to free vibration using boundary elements’, in C. A. Brebbia (ed.),
Boundary Element Methods in Engineering, springer, Berlin, 1982, pp. 312–326.
5. A. J. Nowak, ‘Temperature elds in domains with heat sources using boundary-only formulation’, Proc. 10th BEM
Conf., Springer, Berlin 1988, pp. 233–247.
6. J. O. Watson, ‘Advanced implementation of the boundary element method for two and three-dimensional elastostatics’,
in P. K. Banerjee and R. Buttereld (eds.), Developments in BEM, Applied Science Publisher, 1979.
7. P. K. Banerjee and R. Buttereld, Boundary Element Methods in Engineering Science, McGraw-Hill, London, 1981.
8. A. Deb and P. K. Banerjee, ‘BEM for general anisotropic 2D elasticity using particular integrals’, Comm. appl. numer.
methods, 6, 111–119 (1990).
9. A. Deb, D. P. Henry and R. B. Wilson, ‘Alternate BEM formulations for 2- and 3-D anisotropic thermoelasticity’, Int.
J. Solids Struct., 27, 1721–1738 (1991).
10. C. L. Tan, Y. L. Gao and F. F. Afagh, ‘Anisotropic stress analysis of inclusion problems using the boundary integral
equation method’, J. Strain Anal., 27, 67–76 (1992).
11. T. A. Cruse, Boundary Element Analysis in Computational Fracture Mechanics, Kluwer Academic, Dordrecht, 1988.
12. S. G. Lekhnitskii, Theory of Elasticity of an Anisotropic Body, Mir, Moscow, 1981.
13. S. G. Lekhnitskii, Anisotropic Plates, Gordon and Breach, London, 1968.
14. H. F. Bueckner and N. Y. Schenectady, ‘The propagation of cracks and the energy of elastic deformation’, Trans.
ASME, 80E, 1225–1230 (1958).
15. J. Martinez and J. Dominguez, ‘On the use of quarter-point boundary elements for stress intensity factor computations’,
Int. j. numer. methods eng., 20, 1941–1950 (1984).
16. C. L. Tan and Y. L. Gao, ‘Boundary element analysis of plane anisotropic bodies with stress concentration and cracks’,
Composite Struct., 20, 17–28 (1992).
17. Q. Ao, ‘Development of boundary methods of solution of anisotropic problems’, Ph.D. Dissertation, University of
Kentucky, 1994.
18. Q. Ao and D. J. Shippy, ‘BEM treatment of body forces in plane orthotropic elastostatics’, J. Eng. Mech., 121,
1130–1135 (1995).
Документ
Категория
Без категории
Просмотров
4
Размер файла
222 Кб
Теги
base, elasticity, matrix, beam, stiffness, equations
1/--страниц
Пожаловаться на содержимое документа