Fourier series method for plane elastic problems of polygonal domain

Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

www.elsevier.com/locate/cma

Fourier series method for plane elastic problems of polygonal domain Jiann-Gang Deng, Fu-Ping Cheng * Department of Civil Engineering, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu, Taiwan, ROC Received 1 February 1999; received in revised form 18 April 2000

Abstract This study employs the Fourier series method based on the edge function approach to solve the plane elastic problem of polygonal domain described by Navier equations. The analytical solutions serve as a set of fundamental solutions for each edge. Superposing the solution function and matching the prescribed boundary conditions in each edge allows the solving of the unknown variables and the analysis of the problem. An extra corner function and regularization technique is utilized to enhance convergence rate. Only one element is required to analyze which polygon domain is convex, however, by dividing a non-convex shape into several convex shapes, the proposed method can be extended to irregular geometrical shape. Numerical examples demonstrate the merits of the developed scheme, as well as its eciency and accuracy. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction Linear elastic isotropic behavior is commonly assumed for stress analysis of solids in civil and mechanical engineering. Additionally, most solution schemes were developed for plane elastic problems and later extended to other ®elds of engineering. However, realistic engineering problems are solved only through numerical methods with the most popular methods being FDM, FEM and BEM. Although these methods have been very successful, their common drawback is their ineciency in computation or problem modeling. These ®ndings provide motivation for developing an ecient numerical method. The Fourier series method based on the edge functions has already been successfully applied to solve various boundary value problems. Bird and Steele [1,2] applied this method to analyze the bending problem of circular plate. Kwok surveyed the mode III fracture mechanics analysis [3], and later extended the problem to cracks in the polygonal shell of revolution [4]. Furthermore, Kang [5] developed a procedure to analyze the polygonal plate-bending problem. Also, Wu [6,7] applied the method in question to the stability of shallow shell. Fu and Steele [8] applied it to the Dugdale crack in a cylindrical shell, and Fu [9] extended this method to the in-plane and bending problems of the shallow shell. Gorman [10] presents vibration analysis of plates by this approach. However, most of these researches were aimed to second- and fourthorder di€erential equations governed by the Laplace's equation. This study employs the Fourier series method to solve a boundary value problem described by two coupled second-order partial di€erential operators, the Navier equations governing the plane elastic problem of the polygonal domain. The basic steps in the method are as follows:

*

Corresponding author.

0045-7825/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 0 0 ) 0 0 3 3 3 - 9

4570

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

1. Find the fundamental solutions that satisfy the plane elastic governing equation. By solving a half-plane elastic problem with periodic boundary inputs, the fundamental solutions based on the Fourier series expansion are derived. These solutions will serve as the edge function for each side of the polygonal. 2. Calculate the interaction coecients in each edge. The solution functions derived in step 1 can be used to calculate the Fourier harmonic coecients for all edges in response to the contribution from each of these edge functions. 3. Superposing the contributions from all edges and matching the prescribed boundary conditions. The summation of the edge functions described by Fourier harmonic involved unknown coecients, which must be matched to the prescribed boundary condition. Therefore, a system of equations is constructed by assembling the interaction coecients and the unknown variables in a matrix form. Although the exponential decaying term of analytical solutions served as a solution function, leading to the matrix being diagonally dominant, the present method is restricted to the domain must be convex. However, by partitioning the domain into several convex sub-domains, an e€ective scheme proposes element connection by satisfying appropriate continuity conditions accounting for why this method can be extended to the irregular domain. The main drawback of the proposed method is its low convergence rate due to the Gibb's phenomenon. Two regularization techniques including arithmetic mean sequences and general character r factor are employed to improve the convergence rate [11,12]. Furthermore, an extra corner function is used to promote the convergence rate in the corner boundary condition. Thus, the proposed method can be utilized more eciently. Similar to the boundary integral method, this approach requires no mesh generation, and the number of elements and degrees of freedom are signi®cantly smaller than the FEM, FDM. The advantage of this approach over boundary elements is that the matrices are well-conditioned. By combining plate-bending problems, this approach can be extended to the plated structures.

2. Governing equation and solution function The governing equations of isotropic elastic continuum for plane elasticity with elastic support are illustrated as follows, where gravitation is absent [13]: orx osxy ‡ kd2 u ˆ 0; ox oy osyx ory ‡ kd2 v ˆ 0 in X; ox oy uN ˆ d1 …s†; vN ˆ d2 …s† on oX1 ;

rN

ˆ d3 …s†;

sN

ˆ d4 …s†

on oX2 ;

…1a† …1b† …1c† …1d†

where, kd2 ˆ ks =E…1 m2 †, ks denotes the spring constant of the elastic support. Regarding numerically treating diculty to the constant term of the solution of Eqs. (1a)±(1d) ks is introduced. An extremely small ks will not induce an accuracy problem in the absence of elastic support. According to Fig. 1, X represents a polygonal domain, and u; m; r; s, denote the displacements and stresses in X. The relations of displacements u; m and stresses r; s in terms of axes x; y are de®ned as: rx ˆ

E

…ux ‡ m  vy †;

…2a†

…vy ‡ m  ux †; m2 E …uy ‡ vx †; sxy ˆ 2…1 ‡ m†

…2b†

ry ˆ

1 1

E

m2

…2c†

where ux ; uy ; vx ; vy are derivatives of displacements u; v, and E; m denote the elastic modulus and Poisson's ratio, respectively.

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

4571

Fig. 1. Edge coordinate system of analysis plane.

Four kinds of possibly prescribed boundary values in the outward normal direction to the edge include displacements uN ; vN and stresses rN ; sN . Two of these four quantities are ranked to construct appropriate boundary conditions. The four types of boundary condition used are         uN rN vN rN BC1 ˆ ; BC2 ˆ ; BC3 ˆ ; BC4 ˆ : vN sN sN uN By applying the method of separation variables, the general solution of Eqs. (1a)±(1d) can be derived as follows, herein, the exponential decay part of the solution is used for the half-plane problem. uˆ

1 X

an e

pn y

sin…kn x† ‡ bn e

nˆ1

v ˆ a0 e

q0 x

‡ b0 e

p0 y

‡

1 X

qn y

an c n e

sin…kn x†; pn y

…3a†

cos…kn x† ‡ bn dn e

qn y

cos…kn x†;

…3b†

nˆ1

where p0 ˆ kd ;

q0 ˆ k d

cn ˆ pn =kn ;

p 2=…1 m†;

pn ˆ

q k2n ‡ kd2 ;

qn ˆ

q k2n ‡ kd2 2=…1 m†;

dn ˆ qn =kn ; kn ˆ np=li :

li represents the length of edge i, which the boundary inputs represented by Fourier cosine expansions along yi ˆ 0 as in Fig. 2. kn represents the mode number of li , and unknown coef®cients an ; bn are the amplitude of mode n.

Fig. 2. Coordinate transform on edge i and j.

4572

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

3. Construct Fourier series coecient matrix According to Fig. 1, each edge has de®ned a local coordinate system that represents a half-plane problem. Additionally, the displacement and stress on each edge must satisfy the solutions of governing equations. Therefore, the problem can be calculated by superposition of the edge solutions from all the half-plane problems. However, the boundary values on each edge are contributed by each local solution, and hence, the in¯uence coecients on each edge can be calculated from all local coordinate. As an example of edge i; j (Fig. 2), where displacements uj ; vj exist on edge j, the in¯uence of displacements and stresses on edge i can be calculated from the coordinate transformation as uij ˆ uj cos cij ‡ vj sin cij ;

…4a†

vij ˆ uj sin cij

…4b†



E

i

…rN †j ˆ

vj sin cij ;

m2

1

…1 …sN †ij ˆ

…uj †x …sin2 cij ‡ m  cos2 cij † ‡ …vj †y …cos2 cij ‡ m  sin2 cij †; m†……uj †y ‡ …vj †x † sin cij cos cij ;

 E …vj †x …sin2 cij 2…1 ‡ m†

cos2 cij † ‡ …uj †y …sin2 cij

…vj †y † sin 2cij ;

‡ ……uj †x

…4c† cos2 cij †; …4d†

where cij represents an angle between edge j and i; …uj †x ; …uj †y ; …vj †x ; …vj †y , are derivatives of uj ; vj . The series terms in Eqs. (3a) and (3b), e pn y cos…kn x†; e pn y sin…kn x†; e qn y cos…kn x†; e qn y sin…kn x† are introduced by the following expansions: e

pn y

cos…kn x† ˆ

1 X

…5a†

1 X i …Dpmnn †j cos…am s†;

…5b†

1 X qn i …Cmn †j cos…am s†;

…5c†

mˆ0

e

pn y

sin…kn x† ˆ

i

pn …Cmn †j cos…am s†;

mˆ0

e

qn y

cos…kn x† ˆ

mˆ0

e

qn y

sin…kn x† ˆ

1 X mˆ0

i

i

i

…Dqmnn †j cos…am s†; i

…5d†

i

pn qn where …Cmn †j ; …Dpmnn †j ; …Cmn †j ; …Dqmnn †j are derived in Appendix A. Concerning an arbitrary boundary value dj …x† existing on edge j, this value can be replaced by the Fourier cosine expansion along axis xj .

dj …x† ˆ

1 X …dm †j cos…am xj †;

…6†

mˆ0

where …dm †j represents the mth coef®cient of the function of dj …x†. Herein, the periodically extended function is served as the boundary condition of a half-plane problem. Combining Eqs.(3a) and (3b), Eqs. (4a)±(4d), (5a)±(5d), the displacements and stresses of the mth term of in¯uence functions from edge j to edge i are expressed by Fourier series as follows, where the nhj harmonic terms are employed to analyze the problem.

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

i

……uN †m †j ˆ ……vN †m †ij ˆ i

……rN †m †j ˆ i

……sN †m †j ˆ

nhj X1

i

j

……u1mn †j  ain ‡ …u2mn †j  bin †;

…7a†

……v1mn †ij  ain ‡ …v2mn †ij  bin †;

…7b†

nˆ0 nhj X1 nˆ0

nhj X1

i

i

……s1mn †j  ain ‡ …s2mn †j  bin †;

nˆ0 nhj X1 nˆ0

4573

i

…7c†

i

……t1mn †j  ain ‡ …t2mn †j  bin †:

…7d†

De®nitions of …u1mn †ij ; …u2mn †ij ; …v1mn †ij ; …v2mn †ij ; …s1mn †ij ; …s2mn †ij ; …t1mn †ij ; …t2mn †ij , are derived in Appendix B. For each edge, the coecients of any Fourier series term accumulated from all edges must have the same value as the prescribed boundaries. From Eqs. (6) and Eqs. (7a)±(7d), by simply equating the boundary values, the system of simultaneous equation can be constructed as 2

1

…DMmn †1

6 6 …DM †2 6 mn 1 6 6 3 6 …DMmn †1 6 6. 6. 6. 4 N …DMmn †1

…DMmn †2

1











…DMmn †2

2











…DMmn †32











.. .

.. N

…DMmn †2

…DMmn †ij

.











32 3 2 3 1 …dm †1 …Cn † 76 7 6 7 2 2 6 …dm †2 7 …DMmn †M 7 76 …Cn † 7 7 7 6 6 76 7 6 7 7 3 7 6 …d † …DMmn †3M 76 m 3 7 …C n† 7 ˆ 6 7; 6 76 7 6 7 7 7 7 6 6 . .. 76 .. .. 7 7 6 7 . . 5 5 4 4 5 N N …dm †M …Cn † …DMmn †M 1

…DMmn †M

…8†

where …Cn †i ˆ f ain bin gT , is the unknown coecient of edge i, and …dm †j denotes the prescribed boundary value on edge j. Regarding Dirichlet problems, boundary condition is ……vN †m †j gT

fdm gj ˆ f ……uN †m †j

…9a†

and the corresponding in¯uence coecients matrix is: " i …DMmn †j

ˆ

…u1mn †j

i

…u2mn †j

i

…v1mn †ij

…v2mn †ij

# …9b†

for Neumann problems, boundary condition is: fdm gj ˆ



……rN †m †j

……sN †m †j

T

…9c†

and the corresponding in¯uence coecients matrix is: " …DMmn †ij

ˆ

…s1mn †ij

…s2mn †ij

i

…t2mn †j

…t1mn †j

i

# …9d†

For the mixed boundary conditions, substituting an exact u; v; s; t to calculate the coecients of sub-matrix i …DMmn †j is convenient. A system of simultaneous equations has been developed from Eqs. (9a)±(9d) and can easily be solved by any linear equation solver. The displacements and stresses at any point of domain can be calculated under global coordinates, de®ned as

4574

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

uˆ

N X jˆ1

vˆ

N X jˆ1

rx ˆ

uj cos cj

vj sin cj ;

…10a†

uj sin cj ‡ vj cos cj ;

…10b†

E 1

m2

N h X …uj †x …cos2 cj ‡ m  sin2 cj † ‡ …vj †y … sin2 cj ‡ m  cos2 cj † jˆ1

…1 ry ˆ

E 1

m2

…10c†

N h X …uj †x …sin2 cj ‡ m  cos2 cj † ‡ …vj †y …cos2 cj ‡ m  sin2 cj † jˆ1

‡ …1 sxy ˆ

i m†……vj †x ‡ …uj †y † sin cj cos cj ;

E 2…1 ‡ m†

N h X jˆ1

i m†……vj †x ‡ …uj †y † sin cj cos cj ;

…uj †y …cos2 cj

sin2 cj † ‡ …vj †x …cos2 cj

…10d† sin2 cj † ‡ ……uj †x

i …vj †y † sin 2cj :

…10e†

uj ; vj is de®ned as Eqs. (4a)±(4d), which represents the displacement of point j, but coordinates x; y in Eqs. (4a)±(4d) should be replaced by xj ; yj ; where xj ; yj ; represent the relative coordinate of the unknown point on the jth local coordinate system. Meanwhile, cj indicates the angle between the global coordinate system and jth local coordinate system, and N denotes the total number of de®ned edges. To verify the feasibility of the proposed method, an example with a simple geometry was tested (Fig. 3), comprising a symmetrical plane with two opposite free edges, ®xed in the bottom face, and supporting a 1000 uniform distributed load on the top. Four boundary conditions are assigned as d ˆ T ‰ BC11 BC22 BC23 BC24 Š . The prescribed values in each of the edges are zero except for fBC23 g ˆ T ‰ 1000 0


Š . This example was analyzed by the approach with 33 harmonics and ANSYS [14] with 400 elements. The two results are then compared. According to Fig. 4, the deformed shape is attributed to the H ;  CD;  and the interior of plane EF; G 33 harmonics that are given. Compared with ANSYS along edge AD; (Figs. 5±8), the solutions of these locations are very close for both displacement and stress. In some distance from boundary, both methods yield more accurate results (Figs. 7 and 8). To compare the matrix sizes, 33 harmonics required 272 degrees of freedom, while ANSYS with 400 elements required 882 degrees of

Fig. 3. Plane under uniform axial stress.

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

4575

Fig. 4. Plot of deformed shape.

Fig. 5. Comparison of displacements along edge BC, CD.

freedom, therefore, considerable memory space and CPU time are saved. However, the stress on the corner of ®xed edge AB in this approach di€ers markedly from the ANSYS results. Additionally, the stresses along the ®xed edge oscillate with evident amplitude (Fig. 6). Owing to the non-compatibility of the boundary condition, the Gibb's phenomenon arises in response to the stress discontinuity at the corner of the ®xed edge. Consequently, the low convergence rate of stress will reduce the accuracy and eciency of the proposed method even the more harmonic terms are used. However, the Gibb's phenomenon can be eliminated by regularization techniques. This work introduces the Cesaro procedure and the general character of r factor to smooth the stress oscillation caused by the Gibb's phenomenon. Comparing Figs. 6, 9 and 10, the stress distribution attribute to 33 harmonics by the Cesaro procedure and the r factor has demonstrated that the stress oscillation has been eliminated and the results match the ANSYS solution except near corner A; B.

4576

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

Fig. 6. Comparison of stress along AB.

Fig. 7. Comparison of displacements and stresses along edge EF .

4. Incorporating corner degree of freedom Although the Gibb's phenomenon can be eliminated by the regularization technique, the corner stress will lose it accuracy while reducing the amplitude of high frequency terms. However, this shortcomings can be remedied by incorporating the corner singular term. Based on the analytical solution, an extra corner function is de®ned below: i

…Ccpc †j ˆ cc e …Dpcc †ij ˆ e i

pc yji

pc yji

…Cnpn †j ˆ cn e

cos…kc xij †;

sin…kc xij †;

pn yji

cos…kn xij †;

…11a† …11b† …11c†

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

4577

Fig. 8. Comparison of displacements and stresses along GH .

Fig. 9. Comparison of stress along AB smooth by Lanczos method.

i

…Dpnn †j ˆ e …Cnqn †ij

pn yji

ˆ cn e

i

…Dqnn †j ˆ e

sin…kn xij †;

qn yji

qn yji

…11d†

cos…kn xij †;

…11e†

sin…kn xij †; i

…11f† i

where kc ˆ p=2lj ; …Ccpc †j ; …Dpcc †j represents the in¯uence coecients of corner j to corner i i i i i; …Cnpn †j ; …Dpnn †j ; …Cnqn †j ; …Dqnn †j the in¯uence coef®cients of corner j to edge i, and xij ; yji are the relative coordinates of the starting point of edge i on the jth local coordinate system. From Eqs. (11a)±(11f), the displacement of corner i can be de®ned as: uic

ˆ

N X jˆ1

" i …Dpcc †j



ajc

‡

nhj X1 nˆ0

i ……Dpnn †j

# 

ajn

‡

i …Dqnn †j



bjn †

;

…12a†

4578

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

Fig. 10. Comparison of stress along AB by Lanczos method.

Fig. 11. Comparison of stress along AB by incorporating corner function.

vic

ˆ

N X jˆ1

" i …Ccpc †j



ajc

‡

nhj X1 nˆ0

i ……Cnpn †j

# 

ajn

‡

i …Cnqn †j



bjn †

:

…12b†

Herein, coordinate transformation was used by substituting uic ; vic and their derivatives into Eqs. (4a)±(4d). Therefore, the in¯uence coecients of displacement and stress at corner i can be calculated from corner j and edge j. By rearranging Eqs. (9a)±(9d), for the Dirichlet problems, …DMmn †ij and fdm gj is modi®ed as 2 3 8 9 i i i …u2mn †j ……a21†m †j …u1mn †j < ……uN †m †j = 6 i i i 7 i …v2mn †j ……a31†m †j 5; fdm gj ˆ ……vN †m †j ; …DMmn †j ˆ 4 …v1mn †j …13a† : …v † ; i i c j ……a12†n †j ……a13†n †j …a11†ij

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

4579

for the Neumann problems, …DMmn †ij ; fdm gj is modi®ed as 2

i

…s1mn †j

6 i i …DMmn †j ˆ 6 4 …t1mn †j i ……b12†n †j

i

…s2mn †j i

…t2mn †j i ……b13†n †j

3

9 8 ……rN †m †j > > = < i 7 ……s † † ; ; d ……b31†m †j 7 ˆ f g N m m j j 5 > > ; : …r † i c j …b11† i

……b21†m †j

…13b†

j

i i i i i i these coecients …a11†j ; ……a12†n †j ; ……a13†n †j ; …b11†j ; ……b12†n †j ; ……b13†n †j are de®ned in Appendix B. Addii i i i i i i i i i tionally, ……a21†m †j ˆ …u1cm †j ; ……a31†m †j ˆ …v1cm †j ; ……b21†m †j ˆ …s1cm †j ; ……b31†m †j ˆ …t1cm †j , and …u1cm †j ; …v1cm †j ; i i …s1cm †j ; …t1cm †j are also de®ned in Appendix B, but kn should be replaced by kc . By adding corner degree of

freedom, the unknown coecients vector is  T i …Cn † ˆ ain bin aic :

…13c†

After the corner degree of freedom is included, the previous example was reexamined, the displacement of cornerA; B of the ®xed edge was assigned to zero, and 33 harmonics were utilized for analysis, the Tukey window function [15] is selected to smooth the stress oscillation along the ®xed edge. The stresses of ®xed edge AB are illustrated in Fig. 11, and the corner stress is signi®cantly increased.

5. Element connection The function used herein to construct DM is the solution that holds only for the upper half-plane. Therefore, the problem with a non-convex domain where the divergent terms induced by exponential increase in some parts of the domain is probably found, and the numerical method may lose accuracy and yield unexpected results. However, a non-convex domain can be divided into several convex sub-domains by introducing certain internal boundaries and the continuity conditions on these internal boundaries can be used to connect the sub-domains. For instance, a simple non-convex domain in Fig. 12, is cut into two convex triangular elements. Since boundaries 1,2,4,5, have the prescribed boundary conditions, for the common boundary line 3, the displacement boundary is selected for an unknown variable. From the equilibrium and compatibility condition, the relationship between the common boundary and other prescribed boundary conditions can easily be established. From Eqs. (9a)±(9d), concerning the Dirichlet and Neumann problems, the coecients matrix and prescribed boundary vectors are ‰ DM Š  fC g ˆ fd g;

…14a†

‰ FM Š  fC g ˆ f f g:

…14b†

By eliminating the unknown term fCg of Eqs. (14a) and (14b), the relation between the boundary value fdg and ff g is expressed as: ‰KM Š  fd g ˆ f f g;

…15†

where ‰KMŠ ˆ ‰FMŠ  ‰DMŠ 1 .

Fig. 12. A non-convex domain divided into two convex sub-domains.

4580

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

Instead of using matrix form, the relationship of fdg; ff g of elements I and II; is written as: 2

K11

6 6 K21 4 K31 2 K33 6 6 K43 4 K53

K12 K22 K32 K34 K44 K54

3 I 8 9I 8 9I f1 > > > > > d1 > > = < > 7< = K23 7 d ˆ f ; 2 2 5> > > > > > ; ; : > : > K33 d3 f3 3II 8 9II 8 9 K35 > f3 > > > > d3 > > = < > 7 < = :II K45 7 d ˆ f 4 5 > 4> > > > > > > : ; : ; K55 d5 f5 K13

…16a†

…16b†

The equilibrium equations of the common boundary line 3 of elements I and II are I

I

I

‰K31 Š fd1 gI ‡ ‰K32 Š fd2 gI ‡ ‰K33 Š fd3 gI ˆ ff3 gI ; II

II

…17a†

II

‰K33 Š fd3 gII ‡ ‰K34 Š fd4 gII ‡ ‰K35 Š fd5 gII ˆ ff3 g;II

…17b†

where I



fd3 g ˆ

…uN †3 …vN †3

I

I

; ff3 g ˆ



…rN †3 …sN †3

I

II

; fd 3 g ˆ



…uN †3 …vN †3

II

II

; ff 3 g ˆ



…rN †3 …sN †3

II :

…18†

The vector of Fourier harmonics in the proposed approach, such as fdg; ff g, de®ned on any boundary are directional. For example, Fig. 12, illustrates that the direction of the common boundary line 3 in element I is in the opposite direction of element II. The continuity conditions in which all the Fourier harmonics must be transformed are de®ne in the same direction. For the cosine expansion, the transformation can be completed by simply changing the sign of the even harmonics of the vector. Herein, a transformation on any vector fAg is de®ned as n o n o m‡1 Am ˆ … 1†  Am ; …19† II

II

II

II

on the common boundary of element II, while the transformed vectors fuN g ; fvN g ; frN g ; fsN g are all de®ned in the positive direction of line 3. The continuity conditions can be written as: n oII ‡ … uN † 3 ˆ 0; oII  I n …vN †3 ‡ … vN † 3 ˆ 0; oII  I n …rN †3 ˆ … rN † 3 ; oII  I n …sN †3 ˆ … sN †3 : 

…uN †3

I

…20a† …20b† …20c† …20d†

For displacement boundary and force boundary, de®ned as Eq. (18), the continuity conditions can be written more concisely as: n oII fd3 gI ‡ d3 ˆ 0; …21a† n oII I …21b† ff3 g ˆ f3 ; from Eqs. (16a), (16b) and Eqs. (20a)±(20d), the common boundary displacement can be expressed as: ‰@Šfd3 g ˆ fRg;

…22†

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

4581

where

h i ~ ; ‰@Š ˆ ‰K33 ŠI ‡ K 33    h iII h iII ‰RŠ ˆ ‰K31 ŠI fd1 gI ‡ ‰K32 ŠI fd2 gI ‡ K34 fd4 gII ‡ K35 fd5 gII ;

 K~ represent the ``row transformation'' and ``column transformation'' of K as K; " # m … 1† Kmn  Kˆ ; m‡1 … 1† Kmn h i n n‡1 K~ ˆ … 1† Kmn … 1† Kmn :

Fig. 13. Connection of non-convex domain under axial load.

Fig. 14. Comparing of displacements along line EF , GH .

…23a† …23b†

…24a† …24b†

4582

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

Fig. 15. Comparison of stress along ®x edge AB.

Fig. 16. Comparison of stress along connection edge EF .

The example of the non-convex domain, divide into two sub-domains, with uniform distributed load and ®xed edge (Fig. 13). The 65 harmonics was used to analyze the problem, and results were compared to ANSYS with 10164 high-order elements. The results from two methods were compared with the displacement along the top edge and common boundary in Fig. 14. Nearly the same results are obtained in both the two methods. The stresses along the ®xed edge calculated by the proposed method matches the ANSYS solution but exhibit little oscillation near the corner (Fig. 15). Finally, comparing the stresses in the interior of the domain, the proposed method reaches almost the same solution from the ANSYS results (Fig. 16). 6. Conclusions This work presents a simple and straightforward method, using the Fourier series, applying the plane elasticity boundary value problem with arbitrary boundary conditions for the polygonal domain. The

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

4583

solution function of the governing equation with constant coecients is analytically derived. Meanwhile, the ecient Fast Fourier Transform is used to calculate the discretely prescribed boundary conditions. These conditions give the proposed method highly ecient. Additionally, less degrees of freedom are required than most available numerical methods. Finally, the method proposed herein has the following merits: 1. Little memory and computational e€ort is necessary, since no numerical integration is involved and less degrees of freedom is required. 2. Only one element is required for a convex domain, while a non-convex domain can be partitioned into several convex sub-domains. This fact suggests that the proposed method provide meshless input that can save considerable e€ort in problem modeling. 3. Since boundaries are approximated, errors arise only along boundaries. In the interior of the domain, the solution at each point is highly accurate. 4. The Gibb's phenomenon induced by the discrepant boundary conditions results in stress oscillation on the ®xed edge. However, the regularization technique can be improved the problem. 5. After incorporating the corner singularity, the stress convergence rate of corner is signi®cantly enhanced. The current approach has demonstrated the eciency and accuracy in the presented paper. By the incorporation of the plate-bending problem, the proposed approach analyze the plated structure is feasible. pn qn n n Appendix A. The calculation of Cmn , Dpmn , Cmn , Dqmn

The coecients of Eqs. (5a)±(5d) can be calculated by the following integration: For m 6ˆ 1;

pn Cmn

Dpmnn

Z 2 li ˆ e li 0 Z li 2 ˆ e li 0

pn y

cos…kn x† cos…am s† ds;

pn y

sin…kn x† cos…am s† ds:

The fact that, we always calculate the responses on line i due to the input on line j, accounts for why we can omit the superscripts i; j without losing the clarity. With respect to Fig. 2, the coordinate of any point …x; y† on line i can be written as x ˆ x0 ‡ s cos cij ; y ˆ y0 ‡ s sin cij ; where …x0 ; y0 † denotes the coordinate of the starting point of the line i; s represents the distance measured from …x0 ; y0 † along i, and cij is the angle between the linesj and i. Hence, we have pn Cmn ˆe

pn y 0

‰ cos…kn x†…ccm ‡ ccp†

Dpmnn

pn y 0

‰ sin…kn x†…ccm ‡ ccp† ‡ cos…kn x†…ssm ‡ ssp†Š;

ˆe

sin…kn x†…ssm ‡ ssp†Š;

and ccm; ccp; ssm; ssp represent as    1 1 hm e pn sin cij li sin…hm li † ‡ pn s 1 e pn sin cij li cos…hm li † ; 2 2 li …pn s† ‡ hm    1 1 hp e pn sin cij li sin…hp li † ‡ pn s 1 e pn sin cij li cos…hp li † ; ccp ˆ 2 2 li …pn s† ‡ hp

ccm ˆ

  1 1 hm  1 2 2 li …pn s† ‡ hm   1 1 ssp ˆ hp 1 e 2 2 li …pn s† ‡ hp

ssm ˆ

e

pn sin cij li

pn sin cij li

cos…hm li †

cos…hp li †





pn se

pn se

pn sin cij li

pn sin cij li

sin…hm li †

sin…hp li † ;



4584

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

where pn s ˆ pn sin cij ; hm ˆ kn cos cij for

m ˆ 1;

Z

2 li e li 0 Z 2 li ˆ e li 0

pn Cmn ˆ

Dpmnn

am ; hp ˆ kn cos cij ‡ am pn y

cos…kn x† ds;

pn y

sin…kn x† ds:

pn For a constant term, it more easily integrates into Cmn ; Dpmnn as

n ˆ 0;

Z 1 li e li 0 Z 1 li ˆ e li 0

pn Cmn ˆ

p0 x

cos…am s† ds;

Dpmnn

q0 y

cos…am s† ds;

qn to calculate Cmn ; Dqmnn , it only needs to replaces pn with qn in ccm; ccp; ssm; ssp.

Appendix B. The calculation of (a11), (a12), (a13), (b11), (b12), (b13) From Eqs. (3a), (3b), (4a)±(4d), (5a)±(5d), the values i; j are omitted according to the reason has mentioned earlier. Therefore, all of the coecients can be illustrated as pn u1mn ˆ Dpmnn  cos cij ‡ cn  Cmn  sin cij ; qn u2mn ˆ Dqmnn  cos cij ‡ dn  Cmn  sin cij ;

v1mn ˆ Dpmnn  sin cij

pn cn  Cmn  cos cij ;

qn v2mn ˆ Dqmnn  sin cij dn  Cmn  cos cij ;  pn  Et s1mn ˆ  kn …sin2 cij ‡ m  cos2 cij † pn cn …cos2 cij ‡ m  sin2 cij † Cmn 1 m2 ‡ …1 m†  …pn ‡ cn kn †Dpmnn sin cij cos cij ;  qn  Et s2mn ˆ  kn …sin2 cij ‡ m  cos2 cij † qn dn …cos2 cij ‡ m  sin2 cij † Cmn 1 m2 ‡ …1 m†  …qn ‡ dn kn †Dqmnn  sin cij cos cij ;  Et  pn t1mn ˆ …pn ‡ cn kn †…sin2 cij cos2 cij †Dpmnn ‡ …kn ‡ cn pn †Cmn sin 2cij ; 2…1 ‡ m†  Et  qn t2mn ˆ sin 2cij : …qn ‡ dn kn †…sin2 cij cos2 cij †Dqmnn ‡ …kn ‡ dn qn †Cmn 2…1 ‡ m†

For the corner function, it yields h i h i i i …a11† ˆ e pc yj sin…kc xij † sin cij cc e pc yj cos…kc xij † cos cij ; h i h i i i cn e pn yj cos…kn xij † cos cij ; …a12†n ˆ e pn yj sin…kn xij † sin cij h i h i i i dn e qn yj cos…kn xij † cos cij ; …a13†n ˆ e qn yj sin…kn xij † sin cij  Et  kc …sin2 cij ‡ m  cos2 cij † cc pc …cos2 cij ‡ m  sin2 cij † Ccpc 2 1 m ‡ …1 m†…pc ‡ cc kc † sin cij cos cij Dpcc ;  Et  …b12†n ˆ kn …sin2 cij ‡ m  cos2 cij † cn pn …cos2 cij ‡ m  sin2 cij † Cnpn 2 1 m o ‡ …1 m†…pn ‡ cn kn † sin cij cos cij Dpnn ; …b11† ˆ

J.-G. Deng, F.-.P. Cheng / Comput. Methods Appl. Mech. Engrg. 190 (2001) 4569±4585

…b13†n ˆ

4585

 Et  kn …sin2 cij ‡ m  cos2 cij † dn qn …cos2 cij ‡ m  sin2 cij † Cnqn 2 1 m ‡ …1 m†…qn ‡ dn kn † sin cij cos cij Dqnn ;

where Ccpc ˆ e

pc yji

cos…kc xij †

Dpcc ˆ e

pc yji

sin…kc xij †;

Cnpn ˆ e

pn yji

cos…kc xij † Dpnn ˆ e

pn yji

sin…kc xij †;

Cnqn ˆ e

qn yji

cos…kc xij †

qn yji

sin…kc xij †;

Dqnn ˆ e

where t is the thickness of plate. References [1] M.D. Bird, C.R. Steele, Separated solution procedure for bending of circular plate with circular holes, ASME. Appl. Mech. Rev. 44 (1991) 27±35. [2] M.D. Bird, C.R. Steele, A solution procedure for Laplace's equation on multiply connected circular domains, ASME J. Appl. Mech. 59 (1992) 398±404. [3] F.M.W. Kwok, L.C. Kang, C.R. Steele, Mode III fracture mechanics analysis with Fourier series method, ASME. Appl. Mech. Rev. 44 (1991) 166±170. [4] F.M.W. Kwok, Fourier Series Method for Cracks in Polygonal and Shells of Revolution, Ph.D. Thesis, Dissertation, Department of Mechanical Engineering, Stanford University, 1994. [5] L.C. Kang, A Fourier series method for polygonal domains; large element computation for plates, Ph.D. Thesis, Dissertation, Department of Mechanical Engineering, Stanford University, 1992. [6] C.H. Wu, A Fourier series method for stability of shallow shells with one large element computations for plates, Ph.D. Thesis, Dissertation, Department of Mechanical Engineering, Stanford University, 1993. [7] L.C. Kang, C.H. Wu, C.R. Steele, Fourier series for polygonal plate bending; a very large plate element, Appl. Math. Comput. 67 (1995) 157±225. [8] Y.C. Fu, C.R. Steele, Fourier series solution for the Dugdale crack in a cylindrical shell, ASME J. Appl. Mech. 62 (1995) 533±535. [9] Y.C. Fu, A large element computation for polygonal-shaped plates, shallow shells and cylindrical shells with the Fourier series and state-variable method, Ph.D. Thesis, Dissertation, Department of Mechanical Engineering, Stanford University, 1996. [10] D.J. Gorman, Vibration Analysis of Plates by the Superposition Method, World Scienti®c, Singapore, 1999. [11] G.H. Hardy, Divergent Series, Oxford University Press, London, 1949. [12] C. Lanczos, Appl. Analysis, Prentice-Hall, Englewood Cli€s, NJ, 1956. [13] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, McGraw-Hill, New York, 1951. [14] ANSYS version 5.1, User's Manual, Swanson analysis systems, Inc. Houston, PA, 1994. [15] D.G. Childers, A.E. Durling, Digital Filtering and Signal Processing, West Publishing Company, New York, 1975.