- Email: [email protected]
Buckling analysis of shear deformable shallow shells by the boundary element method
Engineering Analysis with Boundary Elements 31 (2007) 361–372 www.elsevier.com/locate/enganabound
Buckling analysis of shear deformable shallow shells by the boundary element method P.M. Baiz, M.H. Aliabadi Imperial College London, South Kensington, SW7 2AZ, London, UK Received 7 February 2006; accepted 5 July 2006 Available online 5 February 2007
Abstract In this work a boundary element (BE) formulation for buckling problem of shear deformable shallow shells is presented. A set of five boundary integral equations are obtained by coupling two-dimensional plane stress elasticity with shear deformable plate bending (Reissner). The domain integrals appearing in the formulation (due to the curvature and due to the domain load) are transferred into equivalent boundary integrals. The BE formulation is presented as an eigenvalue problem, to provide direct evaluation of critical load factors and buckling modes. Several examples are presented. The BE results for a cylindrical shallow shell with different curvatures are compared with other numerical solutions and good agreements are obtained. r 2006 Published by Elsevier Ltd. Keywords: Shallow shell; Curved plate; Eigenvalue; Buckling coefficients; Shear deformable theory
1. Introduction In structural analysis it is common to encounter problems involving stability assessment, particularly in the design of efficient shell structures. Shells are usually designed to carry load primarily as a membrane; but if the membrane state created by the external loading is compressive, there is a possibility that the membrane equilibrium state will become unstable and the structure will buckle. Timoshenko and Gere [1] have presented a comprehensive description of buckling in several classical shell problems. Other important theoretical works dealing with shell buckling are due to Brush and Almroth [2] and Gerard and Becker [3], giving a review of shell buckling. There are numerous papers dealing with the application of the finite element method (FEM) to buckling of shells. Bushnell [4] presents a good example of numerical procedures with FEM for modelling of shell stability problems. Corresponding author. Tel.: +44 20 7594 5077; fax: +44 20 7594 5078.
E-mail address: [email protected] (M.H. Aliabadi). 0955-7997/$ - see front matter r 2006 Published by Elsevier Ltd. doi:10.1016/j.enganabound.2006.07.008
In recent years the boundary element method (BEM) has been successfully applied to the solution of linear stability problems of plate structures. Manolis et al. [5] developed a direct boundary element formulation for linear stability analysis. Syngellakis and Elzein [6] presented an extended boundary element formulation to incorporate different combination of loading and support conditions, while Nerantzaki and Katsikadelis [7] presented a boundary element formulation for buckling of plates with variable thickness. A more general boundary element formulation for linear buckling with different boundary conditions and arbitrary planar shapes was presented by Lin et al. [8]. In all these cases classical or Kirchhoff–Love theory were assumed. The deficiencies of the classical theory for certain problems has been highlighted by Reissner [9]. He showed the classical theory is not in agreement with experimental results for problems with stress concentration, when the hole diameter is of the order of magnitude of the plate thickness; or also in the case of composite shells, where the ratio of Young’s modulus to shear modulus can be very large. In these cases shell theories accounting for transverse shear deformation overcome the problem associated with the application of the classical theory.
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
362
Boundary element formulations with shear deformable theory for the linear buckling problem of thin walled structures are very limited. Purbolaksono and Aliabadi [10] presented a formulation for shear deformable plates with general boundary conditions and arbitrary planar shapes. For the case of shallow shells, Baiz and Aliabadi [11] introduced the first boundary element approach, using a shear deformable formulation in which the domain integrals appearing in the formulation were treated using quadratic isoparametric cells (domain integration), buckling modes and buckling coefficients were obtained for different boundary conditions and geometries. The method presented in [11] although accurate has the disadvantage of requiring domain discretization. This paper presents a boundary only formulation for the buckling analysis of shear deformable shallow shells. In contrast to the buckling formulation [11] the domain integrals due to the shell curvature and due to the pseudo domain load (due to buckling) are transformed into
x3 x2
x1
r1 r2
equivalent boundary integrals. Several examples are presented to demonstrate the accuracy of the proposed BEM formulation. The BEM results are compared with FEM results, and a good agreement is obtained. 2. Boundary integral equations for shear deformable shells Consider an isotropic linear elastic shallow shell, as shown in Fig. 1. The shell has uniform thickness h, Young’s modulus E, Poisson ratio v and shear modulus G ¼ E=2ð1 þ vÞ, with a quadratic middle surface defined by r1 and r2 , which are principal radius of curvature for the shell in the x1 - and x2 -directions, respectively. The index notation used throughout this paper is as follows: the Greek indices ða; b; gÞ will vary from 1 to 2 and Roman indices ði; j; kÞ from 1 to 3. Generalized displacements are represented as wi and ua , where wa denote rotations of the middle surface, w3 denotes the out-of-plane displacement, and ua denote inplane displacements. The generalized tractions are denoted as pi and ta , where pa denote tractions due to the stress couples, p3 denotes the traction due to shear stress resultant and ta denote tractions due to membrane stress resultants as shown in Fig. 2. Equilibrium equations for plate bending and 2-D elasticity can be written in index notation as follows [12]: M ab;b Qa ¼ 0,
(1)
Qa;a kab N ab þ q3 ¼ 0
(2)
and N ab;b þ qa ¼ 0,
(3)
where k11 ¼ 1=r1 , k22 ¼ 1=r2 , k12 ¼ k21 ¼ 0; and ð Þ;b ¼ qð Þ=qxb . N ab represent membrane stresses, M ab represent bending moments, Qa represent shear forces and qi are the body forces, as shown also in Fig. 2.
Fig. 1. Quadratic shallow shell geometry.
p3 u1 t1
u2
x3
w3
w1
t2
w2 p1
p2
x2
x1
Q1 q1
Q2
q2 q3
M22 N11
N12
N12
N22
M12
M11
M12
Fig. 2. Sign convention for generalized displacements, tractions, stresses and body forces.
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
The constitutive equations based on Reissner’s variational theorem of elasticity [13] can be written as follows: 1n 2n wa;b þ wb;a þ wg;g dab , M ab ¼ D (4) 2 1n Qa ¼ Cðwa þ w3;a Þ
(5)
and N ab ¼ B
1n 2n ua;b þ ub;a þ ug;g dab 2 1n
þ B½ð1 nÞkab þ ndab kff w3 ,
ð6Þ
where B ¼ Eh=ð1 n2 Þ is the tension stiffness; D ¼ Eh3 =½12ð1 n2 Þ is the bendingpffiffiffiffiffi stiffness; C ¼ ½Dð1 nÞl2 =2 is the shear stiffness; l ¼ 10=h is called the shear factor; and dab is the Kronecker delta function. The membrane stress N ab is separated into N ðiÞ ab , which is due to ðiiÞ in-plane displacements, and N ab , which is due to curvature and out-of-plane displacement, as follows: 1n 2n ðiÞ ua;b þ ub;a þ ug;g dab , N ab ¼ B 2 1n N ðiiÞ ab ¼ B½ð1 nÞk ab þ ndab kff w3 , this notation for N ab makes the representation more convenient for the derivation of the boundary integral equations. If the constitutive equations (4)–(6) are introduced into Eqs. (1)–(3), the equilibrium equations for shear deformable shallow shell can be rewritten as follows: D q qw1 qw2 qw3 Dr2 w1 þ ð1 þ nÞ þ ¼ 0, Cw1 C 2 qx2 qx2 qx1 qx1 (7) D q qw1 qw2 qw3 ð1 þ nÞ ¼ 0, þ Dr2 w2 Cw2 C 2 qx1 qx2 qx1 qx2 (8) qw1 qw2 þC þ q3 ¼ 0, qx1 qx2
(9)
B q qu1 qu2 Br u1 þ ð1 þ nÞ þ þ q1 ¼ 0, 2 qx2 qx2 qx1
(10)
Cr2 w3 þ C
2
363
and qu1 qu2 Bðnk11 þ k22 Þ qx1 qx2 Bðk211 þ k222 þ 2nk11 k22 Þw3 ,
q3 ¼ q3 Bðk11 þ nk22 Þ
ð14Þ
The integral equations for shear deformable shallow shell problems are derived by considering the integral representations of the governing Eqs. (1)–(3) from the following integral identities: Z ½ðM ab;b Qa ÞW a þ ðQa;a þ q3 ÞW 3 dO ¼ 0 (15) O
and Z ðN ðiÞ ab;b þ qa ÞU a dO ¼ 0,
(16)
O
where U a and W i ði ¼ a; 3Þ are weighting functions and O is the projected domain of a shell on x1 x2 plane, bounded by boundary G. Eq. (15) is an integral representation related to the governing equations for bending and transverse shear stress resultants, while Eq. (16) is an integral representation related to the governing equations for membrane stress resultants. 3. Governing integral equations for buckling of linear elastic shallow shells Eigenvalue solutions are generally used to estimate the critical buckling loads of stiff structures. The main feature of stiff structures is that they carry their design loads primarily by axial or membrane action, rather than by bending action. Therefore, their response usually involves very little deformation prior to buckling. An eigenvalue analysis usually represent the first step in the stability analysis, because the eigenmodes obtained from this procedure can be used in the investigation of sensitivity of the structure to imperfections. In this section, a set of bucking integral equations for solving the linear eigenvalue problem of shallow shell structures is presented. To reduce the linear buckling problem of shallow shells to an eigenvalue formulation, the distribution of membrane stresses N ab should be known when the structure is subjected to the applied load. Later, it is assumed that membrane stresses in the buckling state are l times the value under the initial reference state.
and B q qu1 qu2 ð1 þ nÞ þ Br2 u2 þ q2 ¼ 0, 2 qx1 qx2 qx1 where
q1 ;
q2 ;
q3
3.1. Integral formulations for the linear buckling problem (11)
are equivalent body forces:
q1 ¼ q1 þ Bðk11 þ nk22 Þ
qw3 , qx1
(12)
q2 ¼ q2 þ Bðnk11 þ k22 Þ
qw3 qx2
(13)
The integral representation related to the governing equations for bending and transverse shear stress resultants of a boundary source point x0 can be derived by using the weighted residual method as shown in Dirgantara and Aliabadi [14] from Eq. (15). With a different arrangement (applying the Green’s identity), in order to avoid the use of derivatives of the kernel W ij , Eq. (36) from [14] can be rewritten
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
364
as follows:
Z cij ðx0 Þwj ðx0 Þ þ Pij ðx0 ; xÞwj ðxÞ dGðxÞ G Z Z 1n 0 ¼ W ij ðx ; xÞpj ðxÞ dGðxÞ kab B 2 G O 2n uf;f ðXÞdab ua;b ðXÞ þ ub;a ðXÞ þ 1n Z W i3 ðx0 ; XÞ dOðXÞ kab B½ð1 nÞkab þ ndab kgg w3 ðXÞ ZO 0 W i3 ðx ; XÞ dOðXÞ þ W i3 ðx0 ; XÞq3 ðXÞ dOðXÞ, ð17Þ O
W ij ðx0 ; xÞ
Pij ðx0 ; xÞ
and are the fundamental soluwhere tions for rotations and out-of-plane displacements, and bending and shear tractions, respectively, and can be found in [12] and [14]. cij ðx0 Þ are the jump terms and they are equal to 12dij when x0 is located on a smooth boundary. Z denotes a Cauchy principal value integral, and x0 ; x 2 G, X 2 O are source and field points, respectively. Buckling equations are obtained when a pseudo-transversal body force and a critical load factor l are introduced into the out of plane integral Eq. (17). The result will be a group of integral equations in terms of prebuckling membrane stresses N ab , and buckled shell displacements w3 , as follows: Z cij ðx0 Þwj ðx0 Þ þ Pij ðx0 ; xÞwj ðxÞ dGðxÞ G Z Z 1n 0 ¼ W ij ðx ; xÞpj ðxÞ dGðxÞ kab B 2 G O 2n uf;f ðXÞdab W i3 ðx0 ; XÞ dOðxÞ ua;b ðXÞ þ ub;a ðXÞ þ 1n Z kab B½ð1 nÞkab þ ndab kgg w3 ðXÞW i3 ðx0 ; XÞ dOðXÞ O Z þ l W i3 ðx0 ; XÞq3 ðXÞ dOðXÞ ZO þ l W i3 ðx0 ; XÞðN ab w3;b Þ;a ðXÞ dOðXÞ, ð18Þ
Because the membrane stresses N ab will be already known at this stage; only the derivatives w3;b ðXÞ and w3;ab ðXÞ need to be expressed in terms of w3 ðXÞ (see Appendix A) in order to arrange the eigenvalue problem. After derivatives of the out-of-plane displacement have been expressed in terms of w3 ðXÞ, the last term in Eqs. (18) and (19) is modified so that these equations can be rewritten as follows: Z cij ðx0 Þwj ðx0 Þ þ Pij ðx0 ; xÞwj ðxÞ dGðxÞ G Z W ij ðx0 ; xÞpj ðxÞ dGðxÞ ¼ G Z 1n W i3 ðx0 ; XÞkab B 2 O 2n uf;f ðXÞdab dOðXÞ ua;b ðXÞ þ ub;a ðXÞ þ 1n Z W i3 ðx0 ; XÞkab Bðð1 nÞkab O Z þ ndab kff Þw3 ðXÞ dOðXÞ þ l W i3 ðx0 ; XÞf b ðXÞ dOðXÞ O
ð20Þ and w3 ðX0 Þ ¼
O
W 3j ðX0 ; xÞpj ðxÞ dGðxÞ Z P3j ðX0 ; xÞwj ðxÞ dGðxÞ ZG 1n ua;b ðXÞ þ ub;a ðXÞ kab B 2 O 2n þ uf;f ðXÞdab W 33 ðX0 ; XÞ dOðxÞ 1n Z kab B ð1 nÞkab þ ndab kgg G
O
w3 ðXÞW 33 ðX0 ; XÞ dOðXÞ Z þ l W 33 ðX0 ; XÞf b ðXÞ dOðXÞ,
ð21Þ
f b ¼ q3 þ N ab;a fðrÞ;b F1 w3 þ N ab fðrÞ;a fðrÞ;b F1 w3 .
(22)
O
where ðN ab w3;b Þ;a is a pseudo body force term due to the large deflection of w3 ðXÞ when buckling occurs. The deflection equation w3 at domain points X0 is required as the additional equation to arrange the eigenvalue equation: Z Z w3 ðX0 Þ ¼ W 3j ðX0 ; xÞpj ðxÞ dGðxÞ P3j ðX0 ; xÞwj ðxÞ dGðxÞ G G Z 1n ua;b ðXÞ þ ub;a ðXÞ kab B 2 O 2n uf;f ðXÞdab W 33 ðX0 ; XÞ dOðxÞ þ 1n Z kab B½ð1 nÞkab þ ndab kgg w3 ðXÞ O Z W 33 ðX0 ; XÞ dOðXÞ þ l W 33 ðX0 ; XÞq3 ðXÞ dOðXÞ O Z 0 þ l W 33 ðX ; XÞðN ab w3;b Þ;a ðXÞ dOðXÞ. ð19Þ
Z
O
where
3.2. Membrane integral equations In a similar way as Eq. (17), the integral equations related to the governing equations for membrane stress resultants of a boundary source point x0 can be written as Z 0 0 0 cya ðx Þua ðx Þ þ T ðiÞ ya ðx ; xÞua ðxÞ dGðxÞ G Z þ U ya ðx0 ; xÞB ½kab ð1 nÞ G
þ ndab kff w3 ðxÞnb ðxÞ dGðxÞ Z U ya ðx0 ; XÞB ½kab ð1 nÞ O
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
þ ndab kff w3;b ðXÞ dOðXÞ Z Z ¼ U ya ðx0 ; xÞta ðxÞ dGðxÞ þ U ya ðx0 ; XÞqa dOðXÞ, G
O
ð23Þ 0 T ðiÞ ya ðX ; xÞ
0
where U ya ðX ; xÞ and are the fundamental solutions for in-plane displacements and membrane tractions, respectively, and can be found in [12] and [14]. As in the case of Eqs. (17) and (23) is constructed is such a way that derivatives of the kernels U ya are avoided. This is a useful condition in order to apply the dual reciprocity method. After Eqs. (17) and (23) are used to obtain the boundary values tg ðxÞ, ug ðxÞ, w3 ðxÞ and domain value w3 ðXÞ, membrane stresses at domain points N ab ðXÞ can be calculated. In the absence of membrane body forces (qa ¼ 0), the membrane stress resultants at domain points X0 can be obtained when the displacement integral equation (23) is introduced in the constitutive expression (6): Z U abg ðX0 ; xÞtg ðxÞ dGðxÞ N ab ðX0 Þ ¼ G Z ðiÞ ðX0 ; xÞug ðxÞ dGðxÞ T abg G Z U abg ðX0 ; xÞB ½kgy ð1 nÞ G
þ ndgy kff w3 ðxÞny ðxÞ dGðxÞ Z þ U abg ðX0 ; XÞB ½kgy ð1 nÞ
evaluation of the derivative terms several points uniformly distributed all over the domain O are necessary. Eq. (17) can be written in a discretized form: Z x¼þ1 Ne X 3 X cij ðx0 Þwj ðx0 Þ þ wnl Pij ðx0 ; xðxÞÞFl ðxÞJ n ðxÞ dx j ¼
þ B ½ð1 nÞkab þ ndab kff w3 ðX Þ.
4. Numerical implementation In the present work continuous and semi-discontinuous boundary elements are used; this is because the semidiscontinuous boundary elements help to avoid difficulties with discontinuity of the tractions at corners. In the next subsections a discretized integral equation will be presented, followed by a description for the treatment of the singularities encountered on the developed formulation.
MT X
MT X
x¼1
W ij ðx0 ; xðxÞÞFl ðxÞJ n ðxÞ dx
Bðk11 þ nk22 ÞI Dk 2i
MT X
Bðnk11 þ k22 ÞI Dk 3i
k¼1
kab B½ð1 nÞkab þ ndab kgg I Dk 1i þ
M X
I Dk 4i ,
ð25Þ
k¼1
where N e is the number of boundary elements, Fl are the shape functions, x is the local coordinate and J n is the Jacobian of transformation. M is the number of DRM domain points and M T is the total number of boundary Dk Dk Dk and DRM domain points. I Dk 1i , I 2i , I 3i , I 4i are domain integrals transformed into boundary integrals by the dual reciprocity method. In total, from the integral equations presented in the previous section and in the absence of membrane body forces (qa ¼ 0), the following domain integrals have to be transferred to the boundary: Z Z qu1 D D I 1i ¼ W i3 w3 dO; I 2i ¼ W i3 dO, (26) qx 1 O O Z ¼ O
ID 5a ¼
Z
ID 7ab ¼
O
Z
W i3
qu2 dO; qx2
ID 4i
U a1
qw3 dO; qx1
ID 6a ¼
Z O
U ab1
qw3 dO; qx1
¼ O
W i3 f b dO,
Z O
ID 8ab ¼
U a2 Z O
(27)
qw3 dO, qx2
U ab2
(28)
qw3 dO. qx2
(29)
The transformation of these domain integrals, follows the same procedure as explained in Wen, Aliabadi and Young [16] and Dirgantara and Aliabadi [15]. The discretization of I Dk 1i is given by " Ne X MT 3 X X I Dk cij ðx0 Þw^ 3mj ðx0 Þ p^ 3nl mj 1i ¼ m¼1
Z
n¼1 l¼1 x¼þ1
x¼1 x¼þ1
x¼1
Quadratic isoparametric boundary elements are used to describe the geometry and the function along G, while for the implementation of the dual reciprocity method and the
x¼þ1
k¼1
Z 4.1. Discretization
Z
k¼1
ð24Þ
The kernels U abg and T ðiÞ abg in (24) are linear combination of the first derivatives of U ab and T ðiÞ ab , and can be found in [15]. Due to the presence of the curvature terms, Eq. (20) has to be used simultaneously with the membrane integral Eq. (23).
pnl j
n¼1 l¼1
ID 3i 0
x¼1
n¼1 l¼1 Ne X 3 X
O
þ ndgy kff w3;y ðXÞ dOðXÞ
365
W ij ðx0 ; xðxÞÞFl ðxÞJ n ðxÞ dx þ
Ne X 3 X
w^ 3nl mj
n¼1 l¼1
0 l Pij ðx ; xðxÞÞF ðxÞJ n ðxÞ dx F1 wk3 .
ð30Þ
Dk Dk Integrals I Dk 2i , I 3i , I 4i , and Eqs. (23), (24), (20), (21) have also to be expressed in a similar way as Eqs. (25) and (30), but for the sake of space it will not be shown here.
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
366
The quadratic continuous shape functions for the boundary are defined as F1 ðxÞ ¼ 12xðx 1Þ, F2 ðxÞ ¼ ð1 xÞð1 þ xÞ, F3 ðxÞ ¼ 12xðx þ 1Þ.
ð31Þ
For the case of semi-discontinuous boundary elements, used at the corners, F1S1 ðxÞ F2S1 ðxÞ F3S1 ðxÞ
¼ ¼ ¼
9 10xðx 1Þ; 32ðx 1Þðx 6 2 10xðx þ 3Þ;
obtained (see Appendix A). Finally, the boundary integral equations of the buckling problem (17), (20) and (21) are solved, obtaining buckling modes and buckling load factors. In the following subsections this procedure is explained in detail.
6 F1S3 ðxÞ ¼ 10 xðx 23Þ, þ 23Þ; F2S3 ðxÞ ¼ 32ðx 9 F3S3 ðxÞ ¼ 10 xðx þ 1Þ,
þ 1Þðx 23Þ, ð32Þ
where FlS1 correspond to nodes placed at x ¼ 23; 0; þ1, while FlS3 is for nodes placed at x ¼ 1; 0; þ23. The position of the internal node in semi-discontinuous element is chosen arbitrarily at 23 or þ23; not very close to the element end point to avoid near singularity problems. Finally, the Jacobian of transformation is defined as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qxy ðxÞ qxy ðxÞ J n ðxÞ ¼ , (33) qx qx where qxy ðxÞ=qx is the derivative of the global coordinates xy with respect to the local coordinate x. 4.2. Treatment of the integrals Depending on the integrands, integrals could be regular or singular (when the collocation point belongs to the element over which the integration is performed). In the first case, the integral can be evaluated using a standard gauss quadrature integration; while in the second case special techniques have to be employed. Based on the order of the singularity, different techniques could be used to deal with the singular integral. In this work, near singular integrals (when the collocation node is close to the integration element) are treated with the element subdivision technique. Weakly singular integrals Oðln rÞ are treated using a nonlinear coordinate transformation. Strong singular integrals Oð1=rÞ are computed indirectly by considering the generalized rigid body motion. All these techniques are well established and together with a compilation of some others could be found in Aliabadi [12]. 5. Numerical procedure The numerical procedure to solve the eigenvalue buckling problem of shallow shells is described in this section. First, boundary integral equations (17) and (23) are solved. After displacements and tractions are known for the reference state, the membrane stresses at domain nodes N ab ðXÞ are calculated with Eq. (24). Next, the derivatives of membrane stresses and out of plane displacement are
5.1. Shell in plane stress When Eqs. (17) and (23) are discretized (Section 4.1), and point collocation is performed, the following linear system of equation is obtained for every node on the boundary and domain: 8 9 .> > > > > .. > > > > > > > " # > > p u < H H w=
"
> >u> > > > > > > > . ; :.> . 5ðNþLÞ 8 9 .. > > > > > > >.> > > > > > =
Hw
Hs
> 55ðNþLÞ > >
Gp
0
0
Gs
¼
#
>t> > > > > > >.> > ; :.> . 15NE
> 515NE > >
( ) b þ
0
, ð34Þ
where the vector b is the product of the bending displacement fundamental solutions with the domain load, which in this study are set to zero (q3 ¼ 0 ! b ¼ 0). H and G are boundary element influence for tractions and displacements fundamental solutions, respectively. The indexes p and s on H and G refer to plate bending and plane stress formulations, respectively; while the indexes u and w are coupled terms between plate bending and plane stress formulations. N, L and NE are number of boundary nodes, domain points and boundary elements, respectively. In other words: M ¼ L and M T ¼ L þ N. Once the collocation process has completed, the known and unknown quantities in Eq. (34) are arranged into a set of linear algebraic equation: ½A5ðNþLÞ5ðNþLÞ fXg5ðNþLÞ ¼ fFg5ðNþLÞ ,
(35)
where ½A is the system matrix, fXg contains the displacements and tractions unknowns on the boundary, as well as the displacements in the domain. The vector ½F is obtained by multiplying the related matrices of H or G by the known values of wi , ua or pi , ta . After Eq. (35) is solved, the membrane stresses N ab ðXÞ are calculated from Eq. (24). 5.2. Solving the shell buckling problem Similarly as Section 5.1, after Eqs. (20) and (23) are discretized and point collocation is performed, the following linear system of equations for every collocation node
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
367
on the boundary is obtained:
2
cðx0 Þ þ Hp
Hp
Hpw
Hu
Hu
Hp
cðx0 Þ þ Hp
Hpw
Hu
Hu
Hp
Hp
cðx0 Þ þ Hpw
Hu
Hu
0
0
Hw
cðx0 Þ þ Hs
Hs
0
0
Hw
Hs
cðx0 Þ þ Hs
6 6 6 6 6 6 6 6 4
2
6 6 6 6 ¼6 6 6 6 4 2
Gp
Gp
Gp
0
0
p
p
G
p
G
0
0
Gp
Gp
0
0
G
Gp
6 6 6 6 þ l6 6 6 6 4
s
Gs Gs
0
0
0
G
0
0
0
Gs
Q1 Q2 Q3 0 0
3
5Nþ3L
15NE
8 9 .. > > > > > . > > = < > w3 > > > > > > > ; : .. > . L
7 7 7 7 7 7 7 7 5
8 9 .. > > > . > > > > > > > 3 > > > > > > p > > 1 > > > > 7 > > > > >p > 7 > > 2 > 7 = < > 7 7 p 3 7 > > > > 7 > > > > 7 > t 1 > > > 5 > > > > > > > > > 515NE > t 2 > > > > > > > > > ; : .. > .
8 9 .. > > > . > > > > > > > 3 > > > > > > w > > 1 > > > > > > 7 > > > w2 > 7 > > > 7 = < > 7 7 w 3 7 > > > > 7 > > > > 7 > u 1 > > > 5 > > > > > > > > > > u 2 > > 5ð5Nþ3LÞ > > > > > > > ; : .. > .
ð36Þ
5L
and for every collocation node on the domain X0 , the linear system of equations is given by only the last three equations in (36); which correspond to w3 , u1 and u2 displacements. Qi is formed by the integrals I D 4i , with f b ¼ N ab;a fðrÞ;b F1 þ N ab fðrÞ;a fðrÞ;b F1 . Eq. (36) can be arranged in a similar manner as Eq. (35), to give ½Bð5Nþ3LÞð5Nþ3LÞ fYg5Nþ3L ¼ k½Kð5Nþ3LÞL fw3 gL .
(37)
and (38), homogenous boundary conditions are prescribed (the known values of ua ; wi ; ta and pi on G are zero), resulting in a system where only the pseudo-transversal body force ðN ab w3;b Þ;a is acting. Eq. (37) can be rearranged in term of the unknown vector fYg5Nþ3L , fYg5Nþ3L ¼ k½B1 ð5Nþ3LÞð5Nþ3LÞ ½Kð5Nþ3LÞL fw3 gL ,
(39)
The extra equation (21) has also to be written in a matrix form, similar to Eq. (37), in order to arrange the eigenvalue expression:
where matrix ½B1 is the inverse of matrix ½B. After Eq. (39) is introduced into Eq. (38), the following expression is obtained:
½Ifw3 gL ¼ ½BBLð5Nþ3LÞ fYg5Nþ3L þ k½KKLL fw3 gL ,
½Ifw3 gL ¼ k½BBLð5Nþ3LÞ ½B1 ð5Nþ3LÞð5Nþ3LÞ ½Kð5Nþ3LÞL fw3 gL
(38)
where the matrices ½B and ½BB contain coefficient matrices H and G that are related to the fundamental solutions. Matrix ½I is the identity matrix. Vector fYg represents the unknown boundary conditions and the unknown domain displacements ðua ðXÞ; w3 ðXÞÞ. Vector fw3 g contains the unknown out of plane displacement w3 ðXÞ. Matrices ½K and ½KK are obtained by multiplication of the fundamental solutions with the in-plane stresses N ab ðXÞ and approximation functions f ðrÞ. As it can be seen in Eqs. (37)
þ k½KKLL fw3 gL .
ð40Þ
Eq. (40) can be written as a standard eigenvalue problem equation as follows: 1 ½w ½I fw3 gL ¼ 0, (41) k where ½w ¼ ½BBLð5Nþ3LÞ ½B1 ð5Nþ3LÞð5Nþ3LÞ ½Kð5Nþ3LÞL þ ½KKLL .
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
368
Buckling modes fw3 g and buckling load factors k can be obtained by solving Eq. (41). Each distinct set of values (k ! fw3 g) corresponds to a linear perturbation of the base state. Among these perturbed states we seek the special case for which k has the minimum value that allow the existence of a nontrivial displacement field fw3 g. The freely available source code for Linear Analysis PACKage [17] was used to solve this eigenvalue problem. 6. Numerical examples In this section, some examples will illustrate the capabilities of the present formulation. Cylindrical shallow shells with different curvature, dimensions and boundary conditions will be presented. Analytical solutions for the linear buckling problem of cylindrical shells under the action of uniform axial compression have been presented by several authors. Timoshenko and Gere’s solution was based on a set of three equilibrium equations [1] and Donnell [18] gave a single eighth order partial differential equation in the radial displacement. For the case of curved plates (see Fig. 4) under axial compression, the same method as in the case of a circular cylinder axially compressed has been used for calculating the critical stresses. Basically, these solutions are based on the substitution of trigonometric functions which satisfy some specific boundary conditions and assumed buckling mode. Also, analytical solutions are based on classical shell theory, where shear deformable effects are not considered. Because of these limitations it seems more convenient to compare the results presented here with other numerical solutions, based on the finite element method and the boundary domain element formulation [11]. In order to help the interpretation of the results, the following nondimensional parameters will be used throughout the examples [3]: K¼
N c 12ð1 v2 Þ b2 p2 E h2
(42)
and pffiffiffiffiffiffiffiffiffiffiffiffiffi 1 v2 Z¼ , (43) rh where K represents the buckling coefficient, Z represents the curvature parameter for curved plates, N c represents the critical in plane stress, obtained from the multiplication of the buckling load factor k with the actual applied stress, and b denotes the length of the curved side of the shallow shell. The boundary conditions considered here are simply supported and clamped. In the case of simply supported edges the out-of-plane displacement is set to zero ðw3 ¼ 0Þ while for clamped edges the out-of-plane displacement and the rotations are set to zero (w3 ¼ 0 and wa ¼ 0). The inplane displacements on all cases are set free (ua a0); because these are the directions in which the load will be applied. b2
6.1. Convergence study In this first example, a study of convergence is performed for the case of a simply supported square shallow shell subjected to uniform axial compression, as shown in Fig. 4. The dimensions considered are: a=b ¼ 1, thickness h=b ¼ 0:015, Young’s modulus E and Poisson ratio v ¼ 0:3. The radius of the shell is r=b ¼ 12:5. Tables 1 and 2 present the buckling coefficients for meshes with different number of boundary elements and domain nodes, respectively. The percentages given in the last column, are with respect to the most refined case of 40 boundary elements and 10 10 domain points which gives a buckling coefficient of 4.187. It is clear from Tables 1 and 2 that the increasing number of boundary elements have a small influence on the solution, compared with the number of DRM domain points. For example, in Table 1 an increase from 20 boundary elements to 40 boundary elements produce less than 0.5% improvement, while for Table 2 an increase from 5 5 DRM domain points to 10 10, produce a change of almost 5%. The influence of the number of DRM domain points on the buckling solution could be attributed to the fact (as it was shown during the derivation of the present formulation) that the effect of the pseudo-transversal body force is introduced into the system by the DRM domain points. This means that a higher number of domain points will improve the modelling of the pseudo-transversal body force, and therefore will improve the final solution. Fig. 3 shows 3 different meshes, with boundary elements, as well as uniformly distributed domain points.
Table 1 Buckling coefficients for different boundary meshes
No.
Boundary mesh
DRM points
K
%Diff.
1 2 3 4 5 6
20 24 28 32 36 40
77 77 77 77 77 77
4.226 4.219 4.218 4.217 4.216 4.216
0.923 0.758 0.735 0.711 0.688 0.688
elements elements elements elements elements elements
Table 2 Buckling coefficients for different number of domain points (DRM)
No.
Boundary mesh
DRM points
K
%Diff.
1 2 3 4 5 6
32 32 32 32 32 32
55 66 77 88 99 10 10
4.405 4.342 4.217 4.201 4.190 4.188
4.949 3.570 0.711 0.333 0.072 0.024
elements elements elements elements elements elements
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
20 Elements & 7X7 DRM
32 Elements & 5X5 DRM
369
40 Elements & 10X10 DRM
Fig. 3. Boundary meshes and DRM distribution.
h/2 h/2 r
Nc
a
Nc
b
Fig. 5. Finite element model with simply supported boundary conditions. Fig. 4. Rectangular shallow shell subjected to uniform compression load.
6.2. Cylindrical shallow shell with different curvature In this example rectangular cylindrical shallow shells with different curvature and boundary conditions are subjected to uniform axial compression, as shown in Fig. 4. The shells have a constant aspect ratio a=b ¼ 2, thickness h=b ¼ 0:015, Young’s modulus E and Poisson ratio v ¼ 0:3. The radii considered are r=b ¼ 5; 7:5; 11:5; 17:5; 27:5; 63; plate. Fig. 5 shows the ANSYS [19] model with the 6 12 quadrilateral elements (SHEL93), the applied compressive load and the simply supported boundary conditions for a shallow shell of radius r=b ¼ 11:5. It can be seen from Fig. 5 that several displacement constrains in the membrane directions ðua Þ are also used for stability of the model. The FEM buckling mode for this geometry is presented in Fig. 6.
The present BEM model has 24 boundary elements and 7 14 domain points (98 domain nodes and 48 boundary nodes). The finite element mesh has 72 quadratic quadrilateral elements (253 nodes), and the BEM model from [11] has 18 boundary elements and 18 domain cells (91 internal nodes and 36 boundary nodes). The results are plotted in Fig. 7, together with finite element and boundary domain element solutions [11]. As expected, Fig. 7 shows that an increase in the curvature of the plate also increase the stiffness of the structure, and conversely the decrease of the curvature converge to the plate solution. The results of Fig. 7 are presented in Table 3. From examining the values presented in Table 3 and Fig. 7, it is clear the good agreement between the present formulation and FEM and the boundary domain results [11]. It can also be seen from Table 3 that increasing the curvature of the
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
370
Table 4 Buckling coefficients for the plate solution
BEM FEM BEM [11]
Fig. 6. FEM buckling mode for a simply supported shallow shell of 23 in radius.
3.999 3.937 3.952
Simply supported
Clamped
Exact
%Diff.
Exact
%Diff.
4.0 4.0 4.0
0.025 1.200 1.575
8.0 8.0 8.0
0.550 0.938 1.850
8.044 7.925 7.852
plate results in larger differences with the other solutions (FEM and BEM [11]). This behaviour can be attributed to the way in which the domain integrals in the formulations are treated: FEM and the boundary domain element formulation [11] directly evaluate the domain integrals while the present formulation transforms the domain integrals to equivalent boundary integrals by the use of the dual reciprocity method. Therefore, it is clear that an increment of the curvature (which are mainly domain terms) will increase the difference between the domain (FEM and BEM [11]) and the only boundary element solutions. As a final comment on the results in Table 3, it is important to keep in mind that all these results are obtained by numerical procedures, and that although FEM and BEM [11] agree better between them than with the present formulation, the results from this work cannot be necessarily considered to be less accurate. A good example of this can be observed with the buckling coefficients for the flat plate case, which has the well-known analytical solutions of 4:0 and 8:0 for simply supported and clamped boundary conditions. It is clear from Table 4 that results from the present formulation agree better with the analytical solution than FEM and BEM [11]. 6.3. Cylindrical shallow shell with different aspect ratios
Fig. 7. Buckling coefficients for simple supported and clamped rectangular shallow shells under uniform compression. Table 3 Buckling coefficients for different curvatures and boundary conditions Radius ðr=bÞ Z
Simply supported
Clamped
BEM FEM BEM [11] BEM 5 7.5 11.5 17.5 27.5 63 Plate
12.719 8.479 5.530 3.634 2.313 1.009 0.000
5.982 4.910 4.392 4.170 4.068 4.012 3.999
5.621 4.716 4.274 4.084 3.997 3.948 3.937
5.642 4.734 4.291 4.100 4.012 3.964 3.952
FEM
BEM [11]
10.671 10.453 10.402 9.212 9.048 8.985 8.541 8.403 8.334 8.259 8.131 8.060 8.131 8.009 7.936 8.061 7.943 7.868 8.044 7.925 7.852
In this last example, cylindrical shallow shells with different aspect ratios (a=b) are considered. The material properties and dimensions of the shells are as follows: Young’s modulus E, Poisson ratio v ¼ 0:3, thickness h=b ¼ 0:015, and radius r=b ¼ 17:5. It can be seen from Fig. 8 the good agreement between BEM and FEM results for the whole range of aspects ratios and both set of boundary conditions. Some buckling modes for the simply supported boundary conditions are shown in Fig. 9. As it is expected the buckling modes exhibit the typical number of half waves for each specific case. 7. Conclusions A new boundary element formulation for the solution of the linear eigenvalue buckling problem of shear deformable shallow shells was developed. Buckling integral equations were obtained after a pseudo-transversal body force and a
ARTICLE IN PRESS P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
371
should be focused not only in the linear stability problem, but also in the postbuckling regime (nonlinear analysis). Acknowledgement The first author would like to thank the Department of Aeronautical Engineering at Imperial College London for the scholarship granted. Appendix A. Evaluation of the derivative terms The derivative terms w3;b ðXÞ and w3;ab ðXÞ are approximated by the use of a radial basis function f ðrÞ as follows: w3 ðx1 ; x2 Þ ¼
M X
f ðrÞm Cm .
(44)
m¼1
Fig. 8. Buckling coefficients for different aspect ratios ða=bÞ of rectangular shallow shell under uniform load.
The radial basis function for the present approximation was chosen as pffiffiffiffiffiffiffiffiffiffiffiffiffiffi (45) f ðrÞ ¼ c2 þ r2 , qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 m 2 m 2 where, r ¼ ðx1 x1 Þ þ ðx2 x2 Þ and c ¼ 2. M is the total number of selected points in the domain, which are the same domain points used for the applications of the dual reciprocity method (DRM points) in the domain. Cm are coefficients determined by values at the M domain points: W ¼ F1 w3 .
(46)
The first derivative of deflection w3;b can be expressed as the product of the first derivative of the radial basis function and the coefficients Cm , as follows: w3;b ðx1 ; x2 Þ ¼ fðrÞ;b F1 w3 . Fig. 9. Buckling modes for simply supported shallow shells with different aspect ratios.
loading factor l, were introduced in the governing linear equations for bending and transverse shear stress resultant. The membrane stress (N ab ) is obtained from the prebuckling state, making possible present the equations as a standard eigenvalue problem in terms of the buckling deflection (w3 ) and the buckling load factor (l). All the domain integrals were transferred to the boundary by the use of the dual reciprocity method, therefore, modelling requires only a boundary mesh and several uniformly distributed domain points. From the examples presented it can be concluded that for the range of validity of the present formulation (Zt12) good agreement with other numerical solutions is obtained (FEM and BEM [11]). It can also be seen from the results that the boundary element formulation require fewer number of degrees of freedom compared to FEM solutions in order to achieve a good level of accuracy. It is well-known that buckling load in shell structures is sensitive to imperfections, therefore, more BEM research
(47)
And the second derivative of deflection w3;ab can be written in the same way: w3;ab ðx1 ; x2 Þ ¼ fðrÞ;a fðrÞ;b F1 w3 .
(48)
Similar to Eq. (47), the derivative of in-plane stress resultants N ab;a can be expressed as N ab;a ðx1 ; x2 Þ ¼ fðrÞ;a F1 N ab .
(49) 1
In the above expressions, F represent the inverse matrix of F which is formed with the values of f ðrÞ; and fðrÞ;a ispffiffiffiffiffiffiffiffiffiffiffiffiffiffi a matrix formed with the values of f ðrÞ;a ( ¼ ra =ð c2 þ r2 ). References [1] Timoshenko S, Gere JM. Theory of elastic stability. 2nd ed. New York: McGraw-Hill; 1961. [2] Brush DO, Almroth BO. Buckling of bars, plates and shells. New York: McGraw-Hill; 1975. [3] Gerard G, Becker H. Handbook of structural stability part III— buckling of curved plates and shells. NACA TN 3783, Washington, 1957. [4] Bushnell D. Computerized buckling analysis of shells. Dordrecht: Martinus Nijhoff Publishers; 1985.
ARTICLE IN PRESS 372
P.M. Baiz, M.H. Aliabadi / Engineering Analysis with Boundary Elements 31 (2007) 361–372
[5] Manolis GD, Beskos DE, Pineros MF. Beam and plate stability by boundary elements. Comput Struct 1986;22:917–23. [6] Syngellakis S, Elzein A. Plate buckling loads by the boundary element method. Int J Numer Method Eng 1994;37:1763–78. [7] Nerantzaki MS, Katsikadelis JT. Buckling of plates with variable thickness—an analog equation solution. Eng Anal Boundary Elem 1996;18:149–54. [8] Lin J, Duffield RC, Shih HR. Buckling analysis of elastic plates by boundary element method. Eng Anal Boundary Elem 1999;23: 131–7. [9] Reissner E. On bending of elastic plates. Q Appl Math 1947;5: 55–68. [10] Purbolaksono J, Aliabadi MH. Buckling analysis of shear deformable plates by the boundary element method. Int J Numer Methods Eng 2005;62:537–63. [11] Baiz PM, Aliabadi MH. Linear buckling analysis of shallow shells by the boundary domain element method. Comput Modelling Eng Sci 2006;13:19–34.
[12] Aliabadi MH. The boundary element method, vol II: applications to solid and structures. Chichester: Wiley; 2002. [13] Reissner E. On a variational theorem in elasticity. J Math Phys 1950;29:90–5. [14] Dirgantara T, Aliabadi MH. A new boundary element formulation for shear deformable shells analysis. Int J Numer Methods Eng 1999;45:1257–75. [15] Dirgantara T, Aliabadi MH. Dual boundary element formulation for fracture mechanic analysis of shear deformable shells. Int J Solids and Struct 2001;28:7769–800. [16] Wen PH, Aliabadi MH, Young A. Application of dual reciprocity method to plates and shells. Eng Anal Boundary Elem 2000;24:583–90. [17] Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J, et al. LAPACK users’ guide. 3rd ed. Society for Industrial and Applied Mathematics, 1999. [18] Donnell LH. Stability of thin walled tubes under torsion. NACA Report 479, 1933. [19] ANSYS, Inc. Version 9.0, 1994.