- Email: [email protected]
The exact thick arch finite element
PERGAMON
Computers and Structures 68 (1998) 369±379
The exact thick arch ®nite element Przemysl/ aw Litewka, Jerzy Rakowski Institute of Structural Engineering, Poznan University of Technology, ul. Piotrowo 5, 60-965 PoznanÂ, Poland Received 8 January 1997; received in revised form 1 January 1998
Abstract The exact stiness matrix is derived for a curved beam element with constant curvature. The plane two-node six degree-of-freedom element is considered in which eects of ¯exural, axial and shear deformations are taken into account. The analytical shape functions describing radial and tangential displacements as well as cross-section rotations are found in the algebraic±trigonometric form. They contain the coupled in¯uences of shear and membrane eects. Based on these shape functions, using the strain energy formula, the stiness matrix for shear ¯exible and compressible arch element is formulated. Obviously, this element is completely free of shear and membrane locking eects. The advantage of the elaborated element is its applicability to any combination of geometrical properties of the arch structure, e.g. the depth±length ratio of element. In presented numerical examples the shear and membrane in¯uences on the displacements for various cases of boundary conditions and loading are investigated. The results coincide exactly with the analytical ones obtained for continuous arches. # 1998 Elsevier Science Ltd. All rights reserved. Keywords: Thick arch; Finite element; Shape function; Stiness matrix
1. Introduction The question of the choice of the proper ®nite elements to the analysis of the curved structures is a subject of numerous papers and many authors have devoted their eorts to it [1±3]. The elaboration of the simple ®nite element which would give the correct results for the wide range of element geometrical properties still remains the main goal of their work. Also, the improvement of the existing elements appears to be the purpose of many papers. These improvements are intended, for example, to avoid parasitic eects as shear and membrane locking phenomena which can completely destroy the accuracy of the calculations. A broad overview of the literature dealing with this subject is given in ref. [4]. In most of the papers, curved beam elements with polynomial shape functions are used. In ref. [5] the cubic polynomials were adopted as the shape functions to the displacements of thin elements and the improvement of the results was obtained by the application of reduced integration and adoption of a least-square ®t. The use of the linear
shape functions for the separated ®elds of rotations, radial and tangential displacements was the basis for the preparation of the C0 curved beam element, where the reduced integration was adopted in ref. [6]. In ref. [4] the combination of linear, quadratic and cubic shape functions for the components of the displacements was applied. Good convergence of the results was obtained when the elements corrected by the mode decomposition were used. The numerical results for thin shells and deep arches were presented there. The shear and membrane locking eects were investigated for the mixed linear±cubic element and quadratic and cubic isoparametric elements in ref. [7]. The reduced integration improved the accuracy of the results in this case as well. The use of hybrid elements did not reveal the locking phenomena. In ref. [8] the thin curved beam elements described by the Winkler beam model were analysed. Three types of approximate trigonometric shape functions were presented and their accuracy was analysed. However, none of them has the analytical exact form. The modi®ed element with good convergence to the exact results was elaborated. The
0045-7949/98/$19.00 # 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 0 5 1 - 0
370
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
elements with other types of curvature are also the subject of some papers, e.g. in ref. [9] the stiness matrix is derived for a parabolic beam element in which the ¯exural, axial and shear deformation eects are taken into account. In the present paper the full set of the shape functions for the constant curvature beam element is derived. The coupled in¯uences of bending, shear and axial forces are included in them. Based on these shape functions, the exact stiness matrix for shear-¯exible and compressible element is given. The introduction of two parameters d and e enables the simple identi®cation of the shear and membrane in¯uences in the shape functions and in the element stiness matrix. The application of this matrix gives the exact results for arches loaded statically for any element mesh (for any number of ®nite elements used in the calculation). The elaborated element is totally free of the shear and membrane locking eects. Setting d or e to zero enables the analysis of the cases where the shear or membrane in¯uences on the deformations are neglected, respectively. In Section 2, the derivation of the exact shape functions in the algebraic±trigonometric form for the curved element is presented. The coupled ¯exural, shear and membrane eects are contained in them. The comparison with the exact shape functions for the straight beam elements [10] is done. The convergence of the shape functions for both elements is proved in the case when the radius of curvature tends to in®nity. Also, the in¯uence of d (shear eect) and e (membrane eect) on the form of these shape functions is analysed. In Section 3 the exact stiness matrix is presented. Its components are derived based on the exact shape functions. Despite the intricate form of the shape functions the stiness matrix is very simple. In Section 4, the formulation of the equilibrium conditions in the form of the dierence equations for the regular mesh structures is presented. Adopting the regular mesh and using the stiness matrices for two adjacent elements, the equilibrium equations of the nodes are transformed to the equivalent sixth-order dierential equation. The limit transformation (length of element tends to zero) gives the exact form of the dierential equation for a continuous arch. In Section 5 some numerical results are given. The shear and membrane in¯uences on the displacements for the structures with various loading and boundary conditions are examined. Section 6 contains concluding remarks.
Fig. 1 is considered. Let a be the angle of the element (a = 2a0) and a its length (a = Ra). The angular coordinate of element cross-section is x (ÿa0ExE a0) and its corresponding curvilinear coordinate is s. The vector of generalized displacements has a form: qT fu1 ; v1 ; f1 ; u2 ; v2 ; f2 gT ;
fi a ji
where ui are the nodal tangential displacements, vi are the radial displacements and ji are the total rotations of cross-sections (i = 1, 2). The following equilibrium condition for the element is ful®lled: K q F;
2
where K is the element stiness matrix, F is a vector of generalized nodal forces corresponding to (1): FT fH1 ; T1 ; m1 ; H2 ; T2 ; m2 gT ;
mi Mi =a;
3
where Hi, Ti, Mi are the axial forces, shear forces and bending moments acting at the nodal cross-sections, i = 1, 2. The constant cross-section of the considered element has the following characteristics: EI is the ¯exural stiness, GA is the shear stiness and EA is the axial stiness, where E and G are Young's and Kirchho's moduli and I, A is the moment of inertia and crosssection area, respectively. We introduce two additional dimensionless parameters which characterize the shear and membrane eects in the elements (see [10]): 2 EI 1 2 1 v i 2 EI 1 i d ; e 2 2 kGA a k a EA a a 4 where n is Poisson's ratio, k is the shear factor, i is the radius of gyration. The strain energy formula written for the curvilinear coordinate has a form:
2. Exact shape functions The plane two-node curved element with six degreesof-freedom and the radius of curvature R presented in
1
Fig. 1. The analysed arch element.
0
ÿ DA2 s0
as0
ac0 00
ac0 ÿ 2s0
1 ÿ
2 aD1
a s0
00
0
ÿ 12
h
0 1 ÿ c
n
2 D1
o 2 1 ÿ c ÿ as 1 1 a2
N(2) vj
s
0
D sÿ
a c0
2 aD1 00
1 2
h
as0
A0 aD1
ÿ D cA ac0 s0 i ÿ 00 A0 1 D sÿ 1 ÿ c A c0 ÿ D1 2 ac0 ÿ s0 1a c0 n o ÿ a12 1 D21 2 1 ÿ c ÿ as s0
00 00
cs0 ÿ 12 D0
2 D1 2
A0 ÿ 2D 1
ÿ D21 N(2) vv
N(1) vj
A 2D as0 ÿ 1a s0 2
cc0 12 D0
2 D1 2
A0 ÿ 2D 1 A0 D1
0
A ÿ 2D a2 s0 ÿ 12 D00 A00 s0 as ÿ D21 c0 s s0 2
0
ÿ DA1 ac0 a2 c0 12 D00 A00 c0 as ÿ D21 s0 s ÿ c0
A0 2D1 2 D1
N(2) vu
0
0
0
A ac0 ÿ 2s0 ÿ aD 1
h
i
ÿ 12 D0 s ÿ aD2 1 1 ÿ c A00 s0 n o A0 a12 1 D21 2 1 ÿ c ÿ as c0 2D as0 2
00
1 ÿ ÿ D cA as0 c0 i n o 00 1 1 ÿ c A s0 ÿ a2 1 D21 2 1 ÿ c ÿ as c0
A0 D2
A0 D2
0
ÿ DA2 s0
0
as0
ÿ DA2 ac0
1 ÿ c n o 1 ÿ a2 D21 2 1 ÿ c ÿ as 1
0 s
N(1) vv
0
a c0 ÿ 12 D0 A00 c0 as D21 s0 s a2 s0 D21 1 ÿ cs0 12 D0 c C1
2
ÿ D21 ÿ D21
It is worth noting that setting d = 0 means neglecting of the shear eect as well as setting e = 0 neglects the compressibility. In Figs. 2±5 some of the shape functions are presented. They are calculated for the arch element with the following data: a = 1 rad, R = 1.0 m (values of d and e are: d = 0.0156, e = 0.005, if not
C0
c0 cos a0 ; d1 a2 d; e1 a2 e D1 a a sin a 1 d1 e1 ÿ 22 1 ÿ cos a ad1 sin a; D2 a a sin a 1 d1 e1 ÿ 2a sin a 1 e1 ; 1 1 1 1 D0 ; D00 ÿ ; D1 D2 D1 D2 8 A0 1 d1 e1 ; A00 1 ÿ d1 e1 :
N(1) vu
s sin a; s0 sin a0 ; c cos a;
Table 1 The coecients Ck for the radial displacement shape functions v
7
ÿ D00 A00 ac0 h i ÿ 1 0 2 00 A0 1 2 D s ÿ aD1 1 ÿ c A c0 D1 2 ac0 ÿ s0 n o a12 1 D21 2 1 ÿ c ÿ as s0
w j; v; u;
i1
where W(x) indicates the displacement function (i) (i) j = j(x), v = v(x) or u = u(x) and N(i) wj, Nwv, Nwu are the corresponding element shape functions (6) related to the unit nodal displacements ji=1, vi=1, ui=1 with i = 1, 2. The coecients Ck included in the formulas for N(i) are presented in Tables 1±3. The following notation is adopted for simplicity:
as0
ac0
C4
A0 D1 A0 D1
According to the superposition rule the displacements of the cross-section with the coordinate x can be written in the form:
C3
N x C0 C1 x C2 C3 xsin x C4 C5 xcos x: 6
2 X i i i Nwj ji Nwv vi Nwu ui ;
A0 2 ÿ 2D a s0 12 D0 A00 s0 as D21 c0 s 1 0 A 2 2 1 0 2D2 a c0 D1 1 ÿ cc0 ÿ 2 D c ÿ
D00 A00 as0
where j = j(s), v = v(s), u = u(s) are the functions of total rotations, radial and tangential displacements, respectively. For the beginning the auxiliary problem is solved: the clamped±clamped arch subjected to the unit support displacements (ji=1, vi=1, ui=1, i = 1, 2) is considered. For that arch the internal forces are determined using the ¯exibility method. Then, using the principle of virtual work, the displacements functions: j, v, u along the element taking into account the in¯uences of bending, shear and axial forces are obtained. In this way the full set of 18 shape functions related to unit support displacements is derived. All of them are the algebraic±trigonometric functions of the angular coordinate x:
W x
A0 D2
C5
5
A0 ÿ 2D 1 A0 ÿ 2D 1
s
Z
kGA u 2 j ÿ v;s ÿ ds 2 R Zs EA v 2 ds; u;s ÿ 2 R s
EI 2 j ds 2 ;s
C2
Z U
1 D2
1 D2
ÿ D12 as
ÿ D12
1 D2
N(1) uv
N(1) uj
N(2) uu
N(2) uv
N(2) uj
1 s 2a
ÿ 12
s
h
1 ÿ c
as 1 D1 ÿ 2 2 1ÿc ÿ D1
ÿ a22
ÿ D21
2 D1
ÿ 2 1ÿc D1
2 as a2 D1
h
ÿ D21 1 ÿ c
ÿ D21 s
C1
a c0
2
a s0 s0
A0 2D2
D21 as0 1 ÿ cc0 c0 h i 2s A0 ÿ 12 D0 D00 c ÿ aD A00 c0 2D as0 ÿ sa0 1n 2 o ÿ c0 2 2 2 D1 c0 ÿ a s0 ÿ a2 D1 2 1 ÿ c ÿ as 1
ÿ
a sc0 ÿ
2
A0 2D2
D00 cA00 ac0
ÿ 12 D00 sA00 ac0
ÿ D21
1 0 2 D
D21 as0 1 ÿ cc0 h i 1 0 00 2s 00 A0 2 D c D ÿ aD1 A c0 2D2 as0 n o ÿ D21 2a s0 ÿ c0 ac02 D21 2 1 ÿ c ÿ as 1
A ÿ 12 D0 sA00 ac0 2D a2 c0 2
0
A D21 a sc0 ÿ 2D a2 s0 2
0
ÿ 12 D0 c D00 A00 ac0
C2
a N(1) ju a N(1) jv N(1) jj a N(2) ju a N(2) jv N(2) jj ÿ D12 a2 s 1 2 D2 a 1 c 1 1 D2 as 2 1 2 ÿ D2 a s ÿ D12 a2 1 c 1 1 D2 as 2
C0 ÿ D21 as ÿ D2n a 1 ÿ c o 1 ÿ 1a D21 2 1 ÿ c ÿ as 1 2 D1 as ÿnD21 a 1 ÿ c o 1 2 a D1 2 1 ÿ c ÿ as 1
C1
Table 3 The coecients Ck for the rotation shape functions j
a 1 c
1 s 2a
a 1 c
ÿ D12 as
N(1) uu
C0
Table 2 The coecients Ck for the tangential displacement shape functions u
2
a c0 a2 s0 2s0 ÿ ac0 ÿ D21 a2 c0 2 2 D1 a s0 ÿ D21 2s0 ÿ ac0
2 D1 2 D1 2 D1
C2
0
ac0
as0
ÿ DA2 s0
A0 D2
A0 D2
0
ÿ DA2 s0
0
as0
ÿ DA2 ac0
A D2
0
C3 as0
0 0 0 0 0 0
C3
A0 2D1
2 D1 s0 s
c0
2
a s0 ÿ D22 a2 c0 ÿ D22 as0 2 2 D2 a s0 2 2 D2 a c0 ÿ D22 as0
2 D2
C4
2
0 0 0 0 0 0
C5
a s0 h i 1 2 ÿ D1 1 ÿ cs0 D12 ac0 ÿ s0 h i 1 0 00 2s 00 s0 ÿ D22 s0 ÿ ca0 2 D D c ÿ aD1 A n o s0 A0 2D1 ac0 ÿ 2s0 a2 D21 2 1 ÿ c ÿ as 1
a c0
2
1 00 00 2 D sA as0
A0 ÿ 2D 1
ÿ 12 D0 D00 cA00 as0 D22 as0
h i ÿ2 D11 1 ÿ cs0 D12 ac0 h i 2s A00 s0 ÿ D22 s0 ÿ 12 D0 c D00 ÿ aD 1 n o 0 A ÿ 2D ac0 ÿ 2s0 ÿ as02 D21 2 1 ÿ c ÿ as 1 1
1 0 00 00 2 2 D c D A as0 D2 A0 2 2 2D1 a c0 ÿ D1 s0 s 1 0 00 A0 2 2 D sA as0 2D1 a s0
C4 0
0
A ÿ aD ac0 ÿ 2s0 1
0
ac0
ac0 ÿ 2s0
ÿ DA1 as0
A0 D1
A0 aD1
0
ÿ DA1 as0
ÿ DA1 ac0
C5
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
373
The shape functions caused by the unit displacements of the second node can be also obtained from those above using the following formulas: 2 1 x Nvv ÿx; Nvv 2 1 x ÿNvu ÿx; Nvu 2 Nuj x 2 Njv x 2 Nju x
1 Nuj ÿx; 1 ÿNjv ÿx; 1 Nju ÿx
2 1 Nvj x ÿNvj ÿx; 2 1 Nuv x ÿNuv ÿx; 2 1 Nuu x Nuu ÿx; 2 1 Njj x Njj ÿx;
and these results coincide with the shape functions N(2) related to coecients Ck given in Tables 1±3. In Fig. 3(a, b) the comparison between the shape functions N(1) vv for the curved and straight beam elements is presented. In Fig. 3(a) the shear eect is taken into account (d = 0.0156) and the membrane eect is neglected and in Fig. 3(b) both eects are neglected. The solid lines indicate the shape functions for the curved beam element and the dashed lines for the
Fig. 2. (a) Shape functions of radial displacement v; (b) shape functions of tangential displacement u; (c) shape functions of rotation j.
speci®ed), x is a dimensionless angular coordinate (x = x/a). In Fig. 2(a±c) the shape functions of radial displacement v, tangential displacement u and rotation j, respectively due to the unit displacements of the node 1 are presented. Solid lines indicate the shape functions related to the unit radial displacement v1, dashed lines, the unit tangential displacement u1 and dotted lines, the unit rotation j1.
Fig. 3. (a) Shape functions N(1) vv for arch and beam (d$0, e = 0); (b) shape functions N(1) vv for arch and beam (d = 0, e = 0).
374
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
1 Nvv
Fig. 4. Shape functions N(1) uu for arch and beam (d$ 0, e = 0).
straight beam element. It is worth to note the considerable dierence between both kinds of elements. In Fig. 4 the similar comparison is done for the shape function N(1) uu for the case when the shear eect is taken into account and the membrane eect is neglected. For the straight beam where the separate ®elds of transverse and longitudinal displacements are adopted the linear shape function for the longitudinal displacement is considered (dashed line). The shape function for the curved beam element despite its algebraic±trigonometric form is nearly linear (solid line). In Fig. 5(a) the in¯uence of the parameter d (shear eect) on the shape function N(1) vj is presented. The membrane eect is neglected (e = 0) and three dierent values of d are considered: solid line indicates the shape function N(1) vj for d = 0.10, dotted line for d = 0.01 and dashed line for d = 0. In Fig. 5(b) the in¯uence of the parameter e (membrane eect) on the same shape function N(1) vj is presented. The shear eect is neglected (d = 0) and three dierent values of e are considered: solid line indicates N(1) vj for e = 0.10, dotted line for e = 0.01 and dashed line for e = 0. Usually, in the statical analysis of the bar structures, the separated ®elds of longitudinal and transverse displacements in the ®nite elements are adopted. For the longitudinal ones, the linear shape functions are used, for the transverse displacements, the thirdorder polynomials, and for the cross-section rotations, the second-order polynomials. The exact shape functions for the shear ¯exible straight elements (Timoshenko beam) were presented in ref. [10]. It is worthwhile investigating the created exact shape functions for the curvilinear element with the assumption that the radius of curvature for the element in hand, tends to in®nity. Let us take the function of transverse displacement:
A0 2 2 a sin a0 1 ÿ cos asin a0 D1 2D1 1 D0 cos a D00 A00 a cos a0 sin x 2 0 A 2 2 a cos a0 1 ÿ cos acos a0 D1 2D1 1 D0 cos a D00 A00 a sin a0 cos x 2 0 2 A 1 ÿ cos a as0 x sin x ÿ D1 D1 0 A a cos a0 x cos x: ÿ D2 ÿ
9
We introduce to (9) the following notation: a = a/R, x = a(x ÿ 1/2)/R (now the origin of the dimensionless coordinate x coincides with node 1). Calculation of the limit for R 4 1 yields: 1 d 1= d 1, where 2x3 ÿ 3x2 ÿ dx limR41 Nvv d=12 d, which exactly coincides with the shape function for the beam element presented in ref. [10]. Making some more complicated operations the correct
Fig. 5. (a) Shape functions N(1) vj for dierent values of d (e = 0); (b) shape functions N(1) vj for dierent values of e (d = 0).
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
results for the other shape functions of straight beam element are obtained [10]: 1 lim Njv
R41
1 lim Nvj
R41
2 6D x ÿ x; a Da 3 2 2 dx; 2x ÿ 4 dx 2
k16 a3 D0 cos a D00 ÿ
3 2 lim Nvv Dÿ2x 3x2 dx;
6D ÿx2 x; a Da 2 ÿ dx; 2x3 ÿ 2 ÿ dx 2
k23 ÿa3 D00 sin a
2 lim Nvj
k24 ÿk15 ;
R41
k12 a4 D00 sin a;
2 2 ÿ 2 ÿ dx lim Njj D3x
k15 a4 D0 sin a; 2 2 a sin a; D1
2 2 a 1 ÿ cos a; D1
k25 ÿa4 D0 cos a ÿ D00 ;
k26 a3 D0 sin a ÿ
R41
2 2 a 1 ÿ cos a; D1
k33 a2 DD00 cos a 1
where 1 : 1 12d
k34 k16 ;
For the other shape functions one gets the following results: 1 2 1 ÿ x; lim Nuu x lim Nuu
R41
2 2 a sin a; D1
k22 a4 D0 ÿ D00 cos a;
2 lim Njv
R41
D
k11 a4 D0 D00 cos a;
k14 ÿa4 D0 cos a D0 ;
1 d 1; 2 ÿ 4 dx lim Njj D3x
R41
where
k13 ÿa3 D0 D00 cos a
R41
375
k35 ÿk26 ;
k36 ÿa2 D0 cos a D00 ÿ 1 ÿ k44 k11 ; k56 ÿk23 ;
R41
4 1 ÿ cos a ÿ a sin a; D1
k45 ÿk12 ;
4 1 ÿ cos a ÿ a sin a; D1
k46 k13 ;
k55 k22 ;
k66 k33 :
and i i i i lim Nuj lim Nju lim Nuv lim Nvu 0;
R41
R41
R41
R41
i 1; 2:
3. The exact element stiness matrix Having substituted the shape functions (7) given in Tables 1±3 to the strain energy formula (5) and using the minimum conditions with respect to the nodal displacements the components of the element stiness matrix are obtained. Here they are presented for the matrix K de®ned in Eq. (2). This matrix is exact, totally free of the shear and membrane locking eects. It contains the in¯uence of the shear forces (important for the cross-sections with the large depth) and the in¯uence of the axial forces (important for arches with small height). Setting d1 or e1 to zero neglects these in¯uences. Using the classical procedure of determining of the stiness matrix from the strain energy formula yields the following: EI K 3 kij ; a
kij kji
10
4. Dierence equation formulation The circular arch, with arbitrary boundary conditions and statical loading acting on it, is considered. We assume that the arch is divided into n identical elements (regular mesh). The equilibrium conditions formulated by FEM may be transformed to the set of three dierential equations. Knowing the element stiness matrix these three equilibrium equations for the node r, which is a common cross-section of two adjacent elements (r ÿ 1, r) and (r, r + 1), can be written in the following form [11]: A1 A2 D2 4A1 fr A3 E ÿ E ÿ1 vr a3 mr EI ÿ A3 E ÿ E ÿ1 fr A6 A7 D2 4A6 vr A4 A5 D2 4A4 u4
a3 Pr EI 2 A4 A5 D 4A4 fr ÿ A8 E ÿ E ÿ1 vr A8 E ÿ E ÿ1 ur
A9 A10 D2 4A9 ur
a3 Hr EI
11
376
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
where a2 1 ÿ cos a ; D2 4a sin a ÿ a2 1 cos a ÿ 4 1 ÿ cos a ÿ 1; D2 2a2 1 ÿ cos a ÿa3 1 ÿ cos a D0 a3 sin a; A4 ; D1 D2 a2 a 1 cos a ÿ 2 sin a a4 1 ÿ cos a ; A6 ; D1 D1 4 a 1 cos a ; A8 ÿD0 a4 sin a; D2 ÿa4 1 cos a ÿA4 =a; A10 ; D1 Mr =a; fr jr a:
A1 A2 A3 A5 A7 A9 mr
In Eq. (11) the nodal loads Mr, Pr and Hr are calculated as equivalent to the given loading according to FEM procedure and stand for concentrated bending moments, radial and tangential forces. Operator E i0E ir is the shifting operator [E i(fr) = fr + i ] and D2 is the second-order dierential operator de®ned as: D20E + E ÿ1ÿ2. Eliminating two unknown nodal displacements in Eq. (11), we can transform the set of equations to one sixth-order dierential equation. For example, after the elimination of jr and ur by means of elementary operations, one gets the simple dierential equation with one unknown (radial displacement of the node vr), which exactly describes the statics of discrete arch with the regular mesh (division into identical elements): D6 4 1 ÿ cos aD4 4 1ÿ cos a2 D2 vr 1 LM LP LH 12 EI where LM ÿ R2 1 ÿ cos aD2 2 1 ÿ cos a2 E ÿ E ÿ1 Mr ; 1 LP 4 fD1 1 ÿ cos a ÿ D2 1 cos a 4a 4 1 ÿ cos a2 D4 4D1 1 ÿ cos a 4 1 ÿ cos a2 D2 gPr ; 1 LH ÿ R4 f4 1 ÿ cos a a ÿ sin a ÿ D1 D2 sin aD2 4a 8a 1 ÿ cos a2 g E ÿ E ÿ1 Hr : In Eq. (12) the fourth-order dierential operator has the form: D40E2+E ÿ2ÿ4(E + E ÿ1) + 6 and the sixthorder operator: D6 0 E3 + E ÿ3 ÿ 6(E2 + E ÿ2) + 15(E + E ÿ1) ÿ 20. It is worth noting that the solution of the Eq. (12) which, of course, is an equivalent formulation for FEM equations, gives the exact values of
the nodal displacements independent of the mesh density. It is also interesting to do the limit transformation of Eq. (12). Dividing both sides of the equation by a6 and setting a 4 0 (the element length a 4 0 as well) one gets [10]: d6 v d4 v d2 v 1 Lm Lp Lh ; 2 dy6 dy4 dy2 EI
13
where y denotes the angular coordinate throughout the whole arch, d3 m dm Mr ÿ R3 ; m m y lim ; 3 a40 dy a dy d2 p EI d4 p EI d2 p Lp R4 2 ÿ R2 R2 ; 4 kGA dy EA dy2 dy Pr p p y lim ; a40 a dh EI d3 h EI d3 h R2 ; Lh ÿR4 R2 3 dy kGA dy EA dy3 Hr ; h h y lim a40 a
Lm ÿR3
which precisely describes the bending of the continuous arch loaded by distributed forces: radially by p(y), tangentially by h(y) and by distributed moments m(y). The following relations were used in the limit transformation: D6 d 4 1 ÿ cos a D4 d4 2 4; lim 6 6 ; lim 2 4 a40 a a40 a a dy dy 4 1 ÿ cos a2 D2 d2 lim 2: a40 a4 a2 dy
5. Numerical results The arch with radius of curvature R = 4 m, the angle o = 2p/3 (length l = 8p/3), the rectangular cross-section with depth h = 0.6 m and width b = 0.4 m, Young's modulus E = 30 GPa and Poisson's ratio v = 0.17 with dierent boundary conditions and loading is investigated. The arch is divided into four elements, hence d = 0.01920 and e = 0.00684. The results of the calculation of the normalized displacements: uC =uC/l, vC =vC/l and jC =jC/o of the arch midpoint C are presented in Table 4. It can be noted that taking into account the shear and membrane eects may be vital for the results of the calculations. In the case of clamped±clamped arch loaded by the vertical force the neglecting of these eects leads to the error of 55% in the vertical displacement vC . The case with the hydrostatic pressure act-
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
377
Table 4 The results of numerical calculations Arch
Displ.
d $0, e $0
d = 0, e$0
d$0, e = 0
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.0000 0.2488 0.0000
0.0000 0.2205 0.0000
0.0000 0.1430 0.0000
0.0000 0.1127 0.0000
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.1252 0.0000 ÿ0.3796
0.1136 0.0000 ÿ0.3683
0.0889 0.0000 ÿ0.4063
0.0773 0.0000 ÿ0.3952
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
ÿ0.0949 0.0000 1.0824
ÿ0.0921 0.0000 1.0375
ÿ0.1016 0.0000 1.0757
ÿ0.0988 0.0000 1.0304
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.0000 0.3047 0.0000
0.0000 0.2546 0.0000
0.0000 0.2010 0.0000
0.0000 0.1748 0.0000
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.2884 0.0000 ÿ0.8064
0.2770 0.0000 ÿ0.8122
0.2460 0.0000 ÿ0.8274
0.2348 0.0000 ÿ0.8332
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
ÿ0.2016 0.0000 1.3613
ÿ0.2030 0.0000 1.3383
ÿ0.2069 0.0000 1.3579
ÿ0.2083 0.0000 1.3350
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.0000 0.1180 0.0000
0.0000 0.1190 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
ing on the arch is interesting because with the membrane eect neglected it reveals zero displacements (circular arch is the neutral axis for this type of loading acting on the shear±¯exible structure). For the clamped±clamped arch loaded by the vertical unit force at its midpoint, the investigation of the d and e in¯uence on the calculated radial displacement of the midpoint vC is performed. All the data used remain the same as above except the fact that the cross-section depth is now varied from 0.1 to 1.0 m,
d = 0, e = 0
hence the values of correlated parameters d and e vary as well. The results are presented in Figs. 6±8. In Fig. 6 the in¯uence of possible neglecting of shear and membrane in¯uences is shown. The upper axis displays values of e corresponding to the given cross-section depth and the lower axis, values of d. The normalized displacement v' is the ratio between the calculated vertical displacement of the centre point vC and its exact value with both shear and membrane in¯uences taken into account. The solid line indicates the case when
378
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
Fig. 8. The in¯uence of setting the parameter d on the midpoint displacement vC for the clamped±clamped arch loaded by vertical force at midpoint. Fig. 6. The in¯uence of membrane and shear eects on the midpoint displacement vC for the clamped±clamped arch loaded by vertical force at midpoint.
shear eect is neglected (d = 0, e $ 0); the dashed line when membrane eect is neglected (d $ 0, e = 0) and the dotted line when both shear and membrane eects are neglected (d = 0, e = 0). With only the shear eect neglected, the calculated values of displacement vC are more accurate than in the case when only the membrane eect is neglected. When the cross-section depth increases, both shear and membrane eects are bigger (values of d and e increase) and the errors caused by the neglecting of these eects increase as well. In Fig. 7 the case when parameter d is taken strictly from the considered cross-section dimensions and e is set to a speci®ed value is presented. The relative displacement v0 is the ratio between the calculated displace-
Fig. 7. The in¯uence of setting the parameter e on the midpoint displacement vC for the clamped±clamped arch loaded by vertical force at midpoint.
ment vC and its value with neglected shear and membrane eects. Solid line indicates the case of e = 0.01, dashed line, e = 0.001 and dotted line, e = 0. In Fig. 8 the case when parameter e is taken strictly from the considered cross-section dimensions and d is set to a speci®ed value is presented. Solid line indicates the case of d = 0.025, dashed line, d = 0.0025 and dotted line, d = 0.
6. Concluding remarks In the present paper, the exact shape functions for the curvilinear beam element including the ¯exural, shear and axial eects are derived. The exact stiness matrix for the curvilinear element with constant curvature is obtained. This element is exact in the sense that for the arch with constant curvature and loading placed in nodes yields the exact results for any mesh density. However, it can be successfully applied for the analysis of any plane arch structures, e.g. the parabolic arch can be easily modelled as an assembly of dierent elements with constant curvature. The elaborated element, despite the fact it contains trigonometric functions, is simple and does not cause the parasitic shear and membrane eects. In the opinion of the authors, it can be successfully developed for the purposes of the shells analysis. The calculated shape functions can create the base for the derivation of mass matrix and geometrical stiness matrix for the analysis of dynamics and stability of the plane curvilinear structures. Such an exact±consistent element can be the starting point for the determination of the approximate element with polynomial shape functions. The results yielding from
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
its application would verify the approximated results coming from the polynomial element. Acknowledgements This work was supported by the grant DS 11-016/95 donated by the Committee of Scienti®c Research. References [1] Ashwell DG, Sabir AB. Further studies in the application of curved ®nite elements to circular arches. Int J Mech Sci 1971;13:507±17. [2] Ashwell DG, Sabir AB. Limitations of certain curved ®nite elements when applied to arches. Int J Mech Sci 1971;13:133±9. [3] Yamada Y, Ezawa Y. On curved ®nite elements for the analysis of curved beams. Int J Numer Meth Engng 1977;11:1635±51.
379
[4] Stolarski H, Chiang MYM. The mode-decomposition, C0 formulation of curved, two-dimensional structural elements. Int J Numer Meth Engng 1989;28:145±54. [5] Pandian N, Appa Rao TVSR, Chandra Satish. Studies on performance of curved beam ®nite elements for analysis of thin arches. Comput Struct 1989;31:997±1002. [6] Prathap G, Bhashyam GR. Reduced integration and the shear±¯exible beam element. Int J Numer Meth Engng 1982;18:195±210. [7] Stolarski H, Belytschko T. Shear and membrane locking in curved C0 elements. Comp Meth Appl Mech Engng 1983;41:279±96. [8] Guimaraes JEF, Heppler GR. On trigonometric basis functions for C1 curved beam elements. Comput Struct 1992;45:405±13. [9] Marquis JP, Wang TM. Stiness matrix of parabolic beam element. Comput Struct 1989;6:863±70. [10] Rakowski J. The interpretation of the shear locking in beam elements. Comput Struct 1990;37:769±76. [11] Rakowski J. A critical analysis of quadratic beam ®nite elements. Int J Numer Meth Engng 1991;341:949±66.
Computers and Structures 68 (1998) 369±379
The exact thick arch ®nite element Przemysl/ aw Litewka, Jerzy Rakowski Institute of Structural Engineering, Poznan University of Technology, ul. Piotrowo 5, 60-965 PoznanÂ, Poland Received 8 January 1997; received in revised form 1 January 1998
Abstract The exact stiness matrix is derived for a curved beam element with constant curvature. The plane two-node six degree-of-freedom element is considered in which eects of ¯exural, axial and shear deformations are taken into account. The analytical shape functions describing radial and tangential displacements as well as cross-section rotations are found in the algebraic±trigonometric form. They contain the coupled in¯uences of shear and membrane eects. Based on these shape functions, using the strain energy formula, the stiness matrix for shear ¯exible and compressible arch element is formulated. Obviously, this element is completely free of shear and membrane locking eects. The advantage of the elaborated element is its applicability to any combination of geometrical properties of the arch structure, e.g. the depth±length ratio of element. In presented numerical examples the shear and membrane in¯uences on the displacements for various cases of boundary conditions and loading are investigated. The results coincide exactly with the analytical ones obtained for continuous arches. # 1998 Elsevier Science Ltd. All rights reserved. Keywords: Thick arch; Finite element; Shape function; Stiness matrix
1. Introduction The question of the choice of the proper ®nite elements to the analysis of the curved structures is a subject of numerous papers and many authors have devoted their eorts to it [1±3]. The elaboration of the simple ®nite element which would give the correct results for the wide range of element geometrical properties still remains the main goal of their work. Also, the improvement of the existing elements appears to be the purpose of many papers. These improvements are intended, for example, to avoid parasitic eects as shear and membrane locking phenomena which can completely destroy the accuracy of the calculations. A broad overview of the literature dealing with this subject is given in ref. [4]. In most of the papers, curved beam elements with polynomial shape functions are used. In ref. [5] the cubic polynomials were adopted as the shape functions to the displacements of thin elements and the improvement of the results was obtained by the application of reduced integration and adoption of a least-square ®t. The use of the linear
shape functions for the separated ®elds of rotations, radial and tangential displacements was the basis for the preparation of the C0 curved beam element, where the reduced integration was adopted in ref. [6]. In ref. [4] the combination of linear, quadratic and cubic shape functions for the components of the displacements was applied. Good convergence of the results was obtained when the elements corrected by the mode decomposition were used. The numerical results for thin shells and deep arches were presented there. The shear and membrane locking eects were investigated for the mixed linear±cubic element and quadratic and cubic isoparametric elements in ref. [7]. The reduced integration improved the accuracy of the results in this case as well. The use of hybrid elements did not reveal the locking phenomena. In ref. [8] the thin curved beam elements described by the Winkler beam model were analysed. Three types of approximate trigonometric shape functions were presented and their accuracy was analysed. However, none of them has the analytical exact form. The modi®ed element with good convergence to the exact results was elaborated. The
0045-7949/98/$19.00 # 1998 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 0 5 1 - 0
370
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
elements with other types of curvature are also the subject of some papers, e.g. in ref. [9] the stiness matrix is derived for a parabolic beam element in which the ¯exural, axial and shear deformation eects are taken into account. In the present paper the full set of the shape functions for the constant curvature beam element is derived. The coupled in¯uences of bending, shear and axial forces are included in them. Based on these shape functions, the exact stiness matrix for shear-¯exible and compressible element is given. The introduction of two parameters d and e enables the simple identi®cation of the shear and membrane in¯uences in the shape functions and in the element stiness matrix. The application of this matrix gives the exact results for arches loaded statically for any element mesh (for any number of ®nite elements used in the calculation). The elaborated element is totally free of the shear and membrane locking eects. Setting d or e to zero enables the analysis of the cases where the shear or membrane in¯uences on the deformations are neglected, respectively. In Section 2, the derivation of the exact shape functions in the algebraic±trigonometric form for the curved element is presented. The coupled ¯exural, shear and membrane eects are contained in them. The comparison with the exact shape functions for the straight beam elements [10] is done. The convergence of the shape functions for both elements is proved in the case when the radius of curvature tends to in®nity. Also, the in¯uence of d (shear eect) and e (membrane eect) on the form of these shape functions is analysed. In Section 3 the exact stiness matrix is presented. Its components are derived based on the exact shape functions. Despite the intricate form of the shape functions the stiness matrix is very simple. In Section 4, the formulation of the equilibrium conditions in the form of the dierence equations for the regular mesh structures is presented. Adopting the regular mesh and using the stiness matrices for two adjacent elements, the equilibrium equations of the nodes are transformed to the equivalent sixth-order dierential equation. The limit transformation (length of element tends to zero) gives the exact form of the dierential equation for a continuous arch. In Section 5 some numerical results are given. The shear and membrane in¯uences on the displacements for the structures with various loading and boundary conditions are examined. Section 6 contains concluding remarks.
Fig. 1 is considered. Let a be the angle of the element (a = 2a0) and a its length (a = Ra). The angular coordinate of element cross-section is x (ÿa0ExE a0) and its corresponding curvilinear coordinate is s. The vector of generalized displacements has a form: qT fu1 ; v1 ; f1 ; u2 ; v2 ; f2 gT ;
fi a ji
where ui are the nodal tangential displacements, vi are the radial displacements and ji are the total rotations of cross-sections (i = 1, 2). The following equilibrium condition for the element is ful®lled: K q F;
2
where K is the element stiness matrix, F is a vector of generalized nodal forces corresponding to (1): FT fH1 ; T1 ; m1 ; H2 ; T2 ; m2 gT ;
mi Mi =a;
3
where Hi, Ti, Mi are the axial forces, shear forces and bending moments acting at the nodal cross-sections, i = 1, 2. The constant cross-section of the considered element has the following characteristics: EI is the ¯exural stiness, GA is the shear stiness and EA is the axial stiness, where E and G are Young's and Kirchho's moduli and I, A is the moment of inertia and crosssection area, respectively. We introduce two additional dimensionless parameters which characterize the shear and membrane eects in the elements (see [10]): 2 EI 1 2 1 v i 2 EI 1 i d ; e 2 2 kGA a k a EA a a 4 where n is Poisson's ratio, k is the shear factor, i is the radius of gyration. The strain energy formula written for the curvilinear coordinate has a form:
2. Exact shape functions The plane two-node curved element with six degreesof-freedom and the radius of curvature R presented in
1
Fig. 1. The analysed arch element.
0
ÿ DA2 s0
as0
ac0 00
ac0 ÿ 2s0
1 ÿ
2 aD1
a s0
00
0
ÿ 12
h
0 1 ÿ c
n
2 D1
o 2 1 ÿ c ÿ as 1 1 a2
N(2) vj
s
0
D sÿ
a c0
2 aD1 00
1 2
h
as0
A0 aD1
ÿ D cA ac0 s0 i ÿ 00 A0 1 D sÿ 1 ÿ c A c0 ÿ D1 2 ac0 ÿ s0 1a c0 n o ÿ a12 1 D21 2 1 ÿ c ÿ as s0
00 00
cs0 ÿ 12 D0
2 D1 2
A0 ÿ 2D 1
ÿ D21 N(2) vv
N(1) vj
A 2D as0 ÿ 1a s0 2
cc0 12 D0
2 D1 2
A0 ÿ 2D 1 A0 D1
0
A ÿ 2D a2 s0 ÿ 12 D00 A00 s0 as ÿ D21 c0 s s0 2
0
ÿ DA1 ac0 a2 c0 12 D00 A00 c0 as ÿ D21 s0 s ÿ c0
A0 2D1 2 D1
N(2) vu
0
0
0
A ac0 ÿ 2s0 ÿ aD 1
h
i
ÿ 12 D0 s ÿ aD2 1 1 ÿ c A00 s0 n o A0 a12 1 D21 2 1 ÿ c ÿ as c0 2D as0 2
00
1 ÿ ÿ D cA as0 c0 i n o 00 1 1 ÿ c A s0 ÿ a2 1 D21 2 1 ÿ c ÿ as c0
A0 D2
A0 D2
0
ÿ DA2 s0
0
as0
ÿ DA2 ac0
1 ÿ c n o 1 ÿ a2 D21 2 1 ÿ c ÿ as 1
0 s
N(1) vv
0
a c0 ÿ 12 D0 A00 c0 as D21 s0 s a2 s0 D21 1 ÿ cs0 12 D0 c C1
2
ÿ D21 ÿ D21
It is worth noting that setting d = 0 means neglecting of the shear eect as well as setting e = 0 neglects the compressibility. In Figs. 2±5 some of the shape functions are presented. They are calculated for the arch element with the following data: a = 1 rad, R = 1.0 m (values of d and e are: d = 0.0156, e = 0.005, if not
C0
c0 cos a0 ; d1 a2 d; e1 a2 e D1 a a sin a 1 d1 e1 ÿ 22 1 ÿ cos a ad1 sin a; D2 a a sin a 1 d1 e1 ÿ 2a sin a 1 e1 ; 1 1 1 1 D0 ; D00 ÿ ; D1 D2 D1 D2 8 A0 1 d1 e1 ; A00 1 ÿ d1 e1 :
N(1) vu
s sin a; s0 sin a0 ; c cos a;
Table 1 The coecients Ck for the radial displacement shape functions v
7
ÿ D00 A00 ac0 h i ÿ 1 0 2 00 A0 1 2 D s ÿ aD1 1 ÿ c A c0 D1 2 ac0 ÿ s0 n o a12 1 D21 2 1 ÿ c ÿ as s0
w j; v; u;
i1
where W(x) indicates the displacement function (i) (i) j = j(x), v = v(x) or u = u(x) and N(i) wj, Nwv, Nwu are the corresponding element shape functions (6) related to the unit nodal displacements ji=1, vi=1, ui=1 with i = 1, 2. The coecients Ck included in the formulas for N(i) are presented in Tables 1±3. The following notation is adopted for simplicity:
as0
ac0
C4
A0 D1 A0 D1
According to the superposition rule the displacements of the cross-section with the coordinate x can be written in the form:
C3
N x C0 C1 x C2 C3 xsin x C4 C5 xcos x: 6
2 X i i i Nwj ji Nwv vi Nwu ui ;
A0 2 ÿ 2D a s0 12 D0 A00 s0 as D21 c0 s 1 0 A 2 2 1 0 2D2 a c0 D1 1 ÿ cc0 ÿ 2 D c ÿ
D00 A00 as0
where j = j(s), v = v(s), u = u(s) are the functions of total rotations, radial and tangential displacements, respectively. For the beginning the auxiliary problem is solved: the clamped±clamped arch subjected to the unit support displacements (ji=1, vi=1, ui=1, i = 1, 2) is considered. For that arch the internal forces are determined using the ¯exibility method. Then, using the principle of virtual work, the displacements functions: j, v, u along the element taking into account the in¯uences of bending, shear and axial forces are obtained. In this way the full set of 18 shape functions related to unit support displacements is derived. All of them are the algebraic±trigonometric functions of the angular coordinate x:
W x
A0 D2
C5
5
A0 ÿ 2D 1 A0 ÿ 2D 1
s
Z
kGA u 2 j ÿ v;s ÿ ds 2 R Zs EA v 2 ds; u;s ÿ 2 R s
EI 2 j ds 2 ;s
C2
Z U
1 D2
1 D2
ÿ D12 as
ÿ D12
1 D2
N(1) uv
N(1) uj
N(2) uu
N(2) uv
N(2) uj
1 s 2a
ÿ 12
s
h
1 ÿ c
as 1 D1 ÿ 2 2 1ÿc ÿ D1
ÿ a22
ÿ D21
2 D1
ÿ 2 1ÿc D1
2 as a2 D1
h
ÿ D21 1 ÿ c
ÿ D21 s
C1
a c0
2
a s0 s0
A0 2D2
D21 as0 1 ÿ cc0 c0 h i 2s A0 ÿ 12 D0 D00 c ÿ aD A00 c0 2D as0 ÿ sa0 1n 2 o ÿ c0 2 2 2 D1 c0 ÿ a s0 ÿ a2 D1 2 1 ÿ c ÿ as 1
ÿ
a sc0 ÿ
2
A0 2D2
D00 cA00 ac0
ÿ 12 D00 sA00 ac0
ÿ D21
1 0 2 D
D21 as0 1 ÿ cc0 h i 1 0 00 2s 00 A0 2 D c D ÿ aD1 A c0 2D2 as0 n o ÿ D21 2a s0 ÿ c0 ac02 D21 2 1 ÿ c ÿ as 1
A ÿ 12 D0 sA00 ac0 2D a2 c0 2
0
A D21 a sc0 ÿ 2D a2 s0 2
0
ÿ 12 D0 c D00 A00 ac0
C2
a N(1) ju a N(1) jv N(1) jj a N(2) ju a N(2) jv N(2) jj ÿ D12 a2 s 1 2 D2 a 1 c 1 1 D2 as 2 1 2 ÿ D2 a s ÿ D12 a2 1 c 1 1 D2 as 2
C0 ÿ D21 as ÿ D2n a 1 ÿ c o 1 ÿ 1a D21 2 1 ÿ c ÿ as 1 2 D1 as ÿnD21 a 1 ÿ c o 1 2 a D1 2 1 ÿ c ÿ as 1
C1
Table 3 The coecients Ck for the rotation shape functions j
a 1 c
1 s 2a
a 1 c
ÿ D12 as
N(1) uu
C0
Table 2 The coecients Ck for the tangential displacement shape functions u
2
a c0 a2 s0 2s0 ÿ ac0 ÿ D21 a2 c0 2 2 D1 a s0 ÿ D21 2s0 ÿ ac0
2 D1 2 D1 2 D1
C2
0
ac0
as0
ÿ DA2 s0
A0 D2
A0 D2
0
ÿ DA2 s0
0
as0
ÿ DA2 ac0
A D2
0
C3 as0
0 0 0 0 0 0
C3
A0 2D1
2 D1 s0 s
c0
2
a s0 ÿ D22 a2 c0 ÿ D22 as0 2 2 D2 a s0 2 2 D2 a c0 ÿ D22 as0
2 D2
C4
2
0 0 0 0 0 0
C5
a s0 h i 1 2 ÿ D1 1 ÿ cs0 D12 ac0 ÿ s0 h i 1 0 00 2s 00 s0 ÿ D22 s0 ÿ ca0 2 D D c ÿ aD1 A n o s0 A0 2D1 ac0 ÿ 2s0 a2 D21 2 1 ÿ c ÿ as 1
a c0
2
1 00 00 2 D sA as0
A0 ÿ 2D 1
ÿ 12 D0 D00 cA00 as0 D22 as0
h i ÿ2 D11 1 ÿ cs0 D12 ac0 h i 2s A00 s0 ÿ D22 s0 ÿ 12 D0 c D00 ÿ aD 1 n o 0 A ÿ 2D ac0 ÿ 2s0 ÿ as02 D21 2 1 ÿ c ÿ as 1 1
1 0 00 00 2 2 D c D A as0 D2 A0 2 2 2D1 a c0 ÿ D1 s0 s 1 0 00 A0 2 2 D sA as0 2D1 a s0
C4 0
0
A ÿ aD ac0 ÿ 2s0 1
0
ac0
ac0 ÿ 2s0
ÿ DA1 as0
A0 D1
A0 aD1
0
ÿ DA1 as0
ÿ DA1 ac0
C5
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
373
The shape functions caused by the unit displacements of the second node can be also obtained from those above using the following formulas: 2 1 x Nvv ÿx; Nvv 2 1 x ÿNvu ÿx; Nvu 2 Nuj x 2 Njv x 2 Nju x
1 Nuj ÿx; 1 ÿNjv ÿx; 1 Nju ÿx
2 1 Nvj x ÿNvj ÿx; 2 1 Nuv x ÿNuv ÿx; 2 1 Nuu x Nuu ÿx; 2 1 Njj x Njj ÿx;
and these results coincide with the shape functions N(2) related to coecients Ck given in Tables 1±3. In Fig. 3(a, b) the comparison between the shape functions N(1) vv for the curved and straight beam elements is presented. In Fig. 3(a) the shear eect is taken into account (d = 0.0156) and the membrane eect is neglected and in Fig. 3(b) both eects are neglected. The solid lines indicate the shape functions for the curved beam element and the dashed lines for the
Fig. 2. (a) Shape functions of radial displacement v; (b) shape functions of tangential displacement u; (c) shape functions of rotation j.
speci®ed), x is a dimensionless angular coordinate (x = x/a). In Fig. 2(a±c) the shape functions of radial displacement v, tangential displacement u and rotation j, respectively due to the unit displacements of the node 1 are presented. Solid lines indicate the shape functions related to the unit radial displacement v1, dashed lines, the unit tangential displacement u1 and dotted lines, the unit rotation j1.
Fig. 3. (a) Shape functions N(1) vv for arch and beam (d$0, e = 0); (b) shape functions N(1) vv for arch and beam (d = 0, e = 0).
374
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
1 Nvv
Fig. 4. Shape functions N(1) uu for arch and beam (d$ 0, e = 0).
straight beam element. It is worth to note the considerable dierence between both kinds of elements. In Fig. 4 the similar comparison is done for the shape function N(1) uu for the case when the shear eect is taken into account and the membrane eect is neglected. For the straight beam where the separate ®elds of transverse and longitudinal displacements are adopted the linear shape function for the longitudinal displacement is considered (dashed line). The shape function for the curved beam element despite its algebraic±trigonometric form is nearly linear (solid line). In Fig. 5(a) the in¯uence of the parameter d (shear eect) on the shape function N(1) vj is presented. The membrane eect is neglected (e = 0) and three dierent values of d are considered: solid line indicates the shape function N(1) vj for d = 0.10, dotted line for d = 0.01 and dashed line for d = 0. In Fig. 5(b) the in¯uence of the parameter e (membrane eect) on the same shape function N(1) vj is presented. The shear eect is neglected (d = 0) and three dierent values of e are considered: solid line indicates N(1) vj for e = 0.10, dotted line for e = 0.01 and dashed line for e = 0. Usually, in the statical analysis of the bar structures, the separated ®elds of longitudinal and transverse displacements in the ®nite elements are adopted. For the longitudinal ones, the linear shape functions are used, for the transverse displacements, the thirdorder polynomials, and for the cross-section rotations, the second-order polynomials. The exact shape functions for the shear ¯exible straight elements (Timoshenko beam) were presented in ref. [10]. It is worthwhile investigating the created exact shape functions for the curvilinear element with the assumption that the radius of curvature for the element in hand, tends to in®nity. Let us take the function of transverse displacement:
A0 2 2 a sin a0 1 ÿ cos asin a0 D1 2D1 1 D0 cos a D00 A00 a cos a0 sin x 2 0 A 2 2 a cos a0 1 ÿ cos acos a0 D1 2D1 1 D0 cos a D00 A00 a sin a0 cos x 2 0 2 A 1 ÿ cos a as0 x sin x ÿ D1 D1 0 A a cos a0 x cos x: ÿ D2 ÿ
9
We introduce to (9) the following notation: a = a/R, x = a(x ÿ 1/2)/R (now the origin of the dimensionless coordinate x coincides with node 1). Calculation of the limit for R 4 1 yields: 1 d 1= d 1, where 2x3 ÿ 3x2 ÿ dx limR41 Nvv d=12 d, which exactly coincides with the shape function for the beam element presented in ref. [10]. Making some more complicated operations the correct
Fig. 5. (a) Shape functions N(1) vj for dierent values of d (e = 0); (b) shape functions N(1) vj for dierent values of e (d = 0).
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
results for the other shape functions of straight beam element are obtained [10]: 1 lim Njv
R41
1 lim Nvj
R41
2 6D x ÿ x; a Da 3 2 2 dx; 2x ÿ 4 dx 2
k16 a3 D0 cos a D00 ÿ
3 2 lim Nvv Dÿ2x 3x2 dx;
6D ÿx2 x; a Da 2 ÿ dx; 2x3 ÿ 2 ÿ dx 2
k23 ÿa3 D00 sin a
2 lim Nvj
k24 ÿk15 ;
R41
k12 a4 D00 sin a;
2 2 ÿ 2 ÿ dx lim Njj D3x
k15 a4 D0 sin a; 2 2 a sin a; D1
2 2 a 1 ÿ cos a; D1
k25 ÿa4 D0 cos a ÿ D00 ;
k26 a3 D0 sin a ÿ
R41
2 2 a 1 ÿ cos a; D1
k33 a2 DD00 cos a 1
where 1 : 1 12d
k34 k16 ;
For the other shape functions one gets the following results: 1 2 1 ÿ x; lim Nuu x lim Nuu
R41
2 2 a sin a; D1
k22 a4 D0 ÿ D00 cos a;
2 lim Njv
R41
D
k11 a4 D0 D00 cos a;
k14 ÿa4 D0 cos a D0 ;
1 d 1; 2 ÿ 4 dx lim Njj D3x
R41
where
k13 ÿa3 D0 D00 cos a
R41
375
k35 ÿk26 ;
k36 ÿa2 D0 cos a D00 ÿ 1 ÿ k44 k11 ; k56 ÿk23 ;
R41
4 1 ÿ cos a ÿ a sin a; D1
k45 ÿk12 ;
4 1 ÿ cos a ÿ a sin a; D1
k46 k13 ;
k55 k22 ;
k66 k33 :
and i i i i lim Nuj lim Nju lim Nuv lim Nvu 0;
R41
R41
R41
R41
i 1; 2:
3. The exact element stiness matrix Having substituted the shape functions (7) given in Tables 1±3 to the strain energy formula (5) and using the minimum conditions with respect to the nodal displacements the components of the element stiness matrix are obtained. Here they are presented for the matrix K de®ned in Eq. (2). This matrix is exact, totally free of the shear and membrane locking eects. It contains the in¯uence of the shear forces (important for the cross-sections with the large depth) and the in¯uence of the axial forces (important for arches with small height). Setting d1 or e1 to zero neglects these in¯uences. Using the classical procedure of determining of the stiness matrix from the strain energy formula yields the following: EI K 3 kij ; a
kij kji
10
4. Dierence equation formulation The circular arch, with arbitrary boundary conditions and statical loading acting on it, is considered. We assume that the arch is divided into n identical elements (regular mesh). The equilibrium conditions formulated by FEM may be transformed to the set of three dierential equations. Knowing the element stiness matrix these three equilibrium equations for the node r, which is a common cross-section of two adjacent elements (r ÿ 1, r) and (r, r + 1), can be written in the following form [11]: A1 A2 D2 4A1 fr A3 E ÿ E ÿ1 vr a3 mr EI ÿ A3 E ÿ E ÿ1 fr A6 A7 D2 4A6 vr A4 A5 D2 4A4 u4
a3 Pr EI 2 A4 A5 D 4A4 fr ÿ A8 E ÿ E ÿ1 vr A8 E ÿ E ÿ1 ur
A9 A10 D2 4A9 ur
a3 Hr EI
11
376
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
where a2 1 ÿ cos a ; D2 4a sin a ÿ a2 1 cos a ÿ 4 1 ÿ cos a ÿ 1; D2 2a2 1 ÿ cos a ÿa3 1 ÿ cos a D0 a3 sin a; A4 ; D1 D2 a2 a 1 cos a ÿ 2 sin a a4 1 ÿ cos a ; A6 ; D1 D1 4 a 1 cos a ; A8 ÿD0 a4 sin a; D2 ÿa4 1 cos a ÿA4 =a; A10 ; D1 Mr =a; fr jr a:
A1 A2 A3 A5 A7 A9 mr
In Eq. (11) the nodal loads Mr, Pr and Hr are calculated as equivalent to the given loading according to FEM procedure and stand for concentrated bending moments, radial and tangential forces. Operator E i0E ir is the shifting operator [E i(fr) = fr + i ] and D2 is the second-order dierential operator de®ned as: D20E + E ÿ1ÿ2. Eliminating two unknown nodal displacements in Eq. (11), we can transform the set of equations to one sixth-order dierential equation. For example, after the elimination of jr and ur by means of elementary operations, one gets the simple dierential equation with one unknown (radial displacement of the node vr), which exactly describes the statics of discrete arch with the regular mesh (division into identical elements): D6 4 1 ÿ cos aD4 4 1ÿ cos a2 D2 vr 1 LM LP LH 12 EI where LM ÿ R2 1 ÿ cos aD2 2 1 ÿ cos a2 E ÿ E ÿ1 Mr ; 1 LP 4 fD1 1 ÿ cos a ÿ D2 1 cos a 4a 4 1 ÿ cos a2 D4 4D1 1 ÿ cos a 4 1 ÿ cos a2 D2 gPr ; 1 LH ÿ R4 f4 1 ÿ cos a a ÿ sin a ÿ D1 D2 sin aD2 4a 8a 1 ÿ cos a2 g E ÿ E ÿ1 Hr : In Eq. (12) the fourth-order dierential operator has the form: D40E2+E ÿ2ÿ4(E + E ÿ1) + 6 and the sixthorder operator: D6 0 E3 + E ÿ3 ÿ 6(E2 + E ÿ2) + 15(E + E ÿ1) ÿ 20. It is worth noting that the solution of the Eq. (12) which, of course, is an equivalent formulation for FEM equations, gives the exact values of
the nodal displacements independent of the mesh density. It is also interesting to do the limit transformation of Eq. (12). Dividing both sides of the equation by a6 and setting a 4 0 (the element length a 4 0 as well) one gets [10]: d6 v d4 v d2 v 1 Lm Lp Lh ; 2 dy6 dy4 dy2 EI
13
where y denotes the angular coordinate throughout the whole arch, d3 m dm Mr ÿ R3 ; m m y lim ; 3 a40 dy a dy d2 p EI d4 p EI d2 p Lp R4 2 ÿ R2 R2 ; 4 kGA dy EA dy2 dy Pr p p y lim ; a40 a dh EI d3 h EI d3 h R2 ; Lh ÿR4 R2 3 dy kGA dy EA dy3 Hr ; h h y lim a40 a
Lm ÿR3
which precisely describes the bending of the continuous arch loaded by distributed forces: radially by p(y), tangentially by h(y) and by distributed moments m(y). The following relations were used in the limit transformation: D6 d 4 1 ÿ cos a D4 d4 2 4; lim 6 6 ; lim 2 4 a40 a a40 a a dy dy 4 1 ÿ cos a2 D2 d2 lim 2: a40 a4 a2 dy
5. Numerical results The arch with radius of curvature R = 4 m, the angle o = 2p/3 (length l = 8p/3), the rectangular cross-section with depth h = 0.6 m and width b = 0.4 m, Young's modulus E = 30 GPa and Poisson's ratio v = 0.17 with dierent boundary conditions and loading is investigated. The arch is divided into four elements, hence d = 0.01920 and e = 0.00684. The results of the calculation of the normalized displacements: uC =uC/l, vC =vC/l and jC =jC/o of the arch midpoint C are presented in Table 4. It can be noted that taking into account the shear and membrane eects may be vital for the results of the calculations. In the case of clamped±clamped arch loaded by the vertical force the neglecting of these eects leads to the error of 55% in the vertical displacement vC . The case with the hydrostatic pressure act-
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
377
Table 4 The results of numerical calculations Arch
Displ.
d $0, e $0
d = 0, e$0
d$0, e = 0
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.0000 0.2488 0.0000
0.0000 0.2205 0.0000
0.0000 0.1430 0.0000
0.0000 0.1127 0.0000
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.1252 0.0000 ÿ0.3796
0.1136 0.0000 ÿ0.3683
0.0889 0.0000 ÿ0.4063
0.0773 0.0000 ÿ0.3952
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
ÿ0.0949 0.0000 1.0824
ÿ0.0921 0.0000 1.0375
ÿ0.1016 0.0000 1.0757
ÿ0.0988 0.0000 1.0304
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.0000 0.3047 0.0000
0.0000 0.2546 0.0000
0.0000 0.2010 0.0000
0.0000 0.1748 0.0000
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.2884 0.0000 ÿ0.8064
0.2770 0.0000 ÿ0.8122
0.2460 0.0000 ÿ0.8274
0.2348 0.0000 ÿ0.8332
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
ÿ0.2016 0.0000 1.3613
ÿ0.2030 0.0000 1.3383
ÿ0.2069 0.0000 1.3579
ÿ0.2083 0.0000 1.3350
uC 10ÿ6 vC 10ÿ6 jC 10ÿ6
0.0000 0.1180 0.0000
0.0000 0.1190 0.0000
0.0000 0.0000 0.0000
0.0000 0.0000 0.0000
ing on the arch is interesting because with the membrane eect neglected it reveals zero displacements (circular arch is the neutral axis for this type of loading acting on the shear±¯exible structure). For the clamped±clamped arch loaded by the vertical unit force at its midpoint, the investigation of the d and e in¯uence on the calculated radial displacement of the midpoint vC is performed. All the data used remain the same as above except the fact that the cross-section depth is now varied from 0.1 to 1.0 m,
d = 0, e = 0
hence the values of correlated parameters d and e vary as well. The results are presented in Figs. 6±8. In Fig. 6 the in¯uence of possible neglecting of shear and membrane in¯uences is shown. The upper axis displays values of e corresponding to the given cross-section depth and the lower axis, values of d. The normalized displacement v' is the ratio between the calculated vertical displacement of the centre point vC and its exact value with both shear and membrane in¯uences taken into account. The solid line indicates the case when
378
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
Fig. 8. The in¯uence of setting the parameter d on the midpoint displacement vC for the clamped±clamped arch loaded by vertical force at midpoint. Fig. 6. The in¯uence of membrane and shear eects on the midpoint displacement vC for the clamped±clamped arch loaded by vertical force at midpoint.
shear eect is neglected (d = 0, e $ 0); the dashed line when membrane eect is neglected (d $ 0, e = 0) and the dotted line when both shear and membrane eects are neglected (d = 0, e = 0). With only the shear eect neglected, the calculated values of displacement vC are more accurate than in the case when only the membrane eect is neglected. When the cross-section depth increases, both shear and membrane eects are bigger (values of d and e increase) and the errors caused by the neglecting of these eects increase as well. In Fig. 7 the case when parameter d is taken strictly from the considered cross-section dimensions and e is set to a speci®ed value is presented. The relative displacement v0 is the ratio between the calculated displace-
Fig. 7. The in¯uence of setting the parameter e on the midpoint displacement vC for the clamped±clamped arch loaded by vertical force at midpoint.
ment vC and its value with neglected shear and membrane eects. Solid line indicates the case of e = 0.01, dashed line, e = 0.001 and dotted line, e = 0. In Fig. 8 the case when parameter e is taken strictly from the considered cross-section dimensions and d is set to a speci®ed value is presented. Solid line indicates the case of d = 0.025, dashed line, d = 0.0025 and dotted line, d = 0.
6. Concluding remarks In the present paper, the exact shape functions for the curvilinear beam element including the ¯exural, shear and axial eects are derived. The exact stiness matrix for the curvilinear element with constant curvature is obtained. This element is exact in the sense that for the arch with constant curvature and loading placed in nodes yields the exact results for any mesh density. However, it can be successfully applied for the analysis of any plane arch structures, e.g. the parabolic arch can be easily modelled as an assembly of dierent elements with constant curvature. The elaborated element, despite the fact it contains trigonometric functions, is simple and does not cause the parasitic shear and membrane eects. In the opinion of the authors, it can be successfully developed for the purposes of the shells analysis. The calculated shape functions can create the base for the derivation of mass matrix and geometrical stiness matrix for the analysis of dynamics and stability of the plane curvilinear structures. Such an exact±consistent element can be the starting point for the determination of the approximate element with polynomial shape functions. The results yielding from
P. Litewka, J. Rakowski / Computers and Structures 68 (1998) 369±379
its application would verify the approximated results coming from the polynomial element. Acknowledgements This work was supported by the grant DS 11-016/95 donated by the Committee of Scienti®c Research. References [1] Ashwell DG, Sabir AB. Further studies in the application of curved ®nite elements to circular arches. Int J Mech Sci 1971;13:507±17. [2] Ashwell DG, Sabir AB. Limitations of certain curved ®nite elements when applied to arches. Int J Mech Sci 1971;13:133±9. [3] Yamada Y, Ezawa Y. On curved ®nite elements for the analysis of curved beams. Int J Numer Meth Engng 1977;11:1635±51.
379
[4] Stolarski H, Chiang MYM. The mode-decomposition, C0 formulation of curved, two-dimensional structural elements. Int J Numer Meth Engng 1989;28:145±54. [5] Pandian N, Appa Rao TVSR, Chandra Satish. Studies on performance of curved beam ®nite elements for analysis of thin arches. Comput Struct 1989;31:997±1002. [6] Prathap G, Bhashyam GR. Reduced integration and the shear±¯exible beam element. Int J Numer Meth Engng 1982;18:195±210. [7] Stolarski H, Belytschko T. Shear and membrane locking in curved C0 elements. Comp Meth Appl Mech Engng 1983;41:279±96. [8] Guimaraes JEF, Heppler GR. On trigonometric basis functions for C1 curved beam elements. Comput Struct 1992;45:405±13. [9] Marquis JP, Wang TM. Stiness matrix of parabolic beam element. Comput Struct 1989;6:863±70. [10] Rakowski J. The interpretation of the shear locking in beam elements. Comput Struct 1990;37:769±76. [11] Rakowski J. A critical analysis of quadratic beam ®nite elements. Int J Numer Meth Engng 1991;341:949±66.