Exact shape functions of imperfect beam element for stability analysis

Advances in Engineering Software 38 (2007) 576–585 www.elsevier.com/locate/advengsoft

Exact shape functions of imperfect beam element for stability analysis R. Adman a

a,*

, H. Afra

b

Built Environment Research Laboratory, Faculty of Civil Engineering, U.S.T.H.B., Algiers, Algeria b National Center of Building Research, Souidania, Algiers, Algeria Received 17 November 2004; accepted 14 August 2006 Available online 27 December 2006

Abstract In this paper, a versatile and efficient solution for the displacement function of a geometrical imperfect beam element under high axial load is derived. The initial imperfection is included by a sine curve approximation. The aim of this work is to propose a natural extension of the element by stability functions. For this purpose, exact displacement shape functions are solved from the differential equilibrium equation, by using arbitrary boundary conditions. These shape functions offer the possibility to express the equivalent nodal forces for any arbitrary transverse loads, and are further extended to obtain the equivalent nodal forces relating to the effect of the initial curvature without affecting the boundary conditions of the element. The accuracy of this element makes the convergence rate better than the conventional cubic element or other used elements. Euler’s critical load is computed readily with no more than two iterations for each step loading with a single element per member. By using these exact shape functions in computer software, the second-order analysis and the nonlinear integrated design become reliable and easy to use for practical purposes.  2006 Elsevier Ltd. and Civil-Comp Ltd. All rights reserved. Keywords: Stability; Second-order; Geometrical imperfection; Beam element; Nonlinear

1. Introduction In the past few decades there has been a large debate about the best element for the modelling of the geometrical nonlinear behaviour of an elastic beam element. In this frame, we can distinguish two main approaches in element formulations: – The element by stability approach. – The element by finite element approach. In the stability function approach the element matrix is developed by solving the governing differential equilibrium equation of a beam-column under the action of an axial load. Different solutions will be obtained for tensile, compressive and zero loads. In the finite element approach, the effect of the geometrical nonlinearity is accounted by *

Corresponding author. Tel.: +213 21 247 914; fax: +213 21 247 992. E-mail address: [email protected] (R. Adman).

including higher order terms in the strain–displacement relationship. The displacement function is approximated by polynomial shape functions which are associated to the nodal displacements and rotations. For instance, we can name the cubic element, the ‘Pointwise Equilibrium Polynomial’ (PEP) proposed by Chan and Zhou [1], etc. Based on the first variation of the energy functional, the element by finite element approach is established. The exactness of the element obtained through the stability approach unfortunately does not take into account the initial geometrical imperfection as well as the effect of any arbitrary load along the member. Its application to design and analysis is therefore limited. On the other hand, from the writers’ experience, the element by finite element approach should only be used when the axial force is small and less than, for instance, half of the buckling load. Although using more elements will improve the accuracy, it will not only require heavy computational time, but divergence problems may also result in some cases due to the inaccurate stiffness matrix used.

0965-9978/$ - see front matter  2006 Elsevier Ltd. and Civil-Comp Ltd. All rights reserved. doi:10.1016/j.advengsoft.2006.08.020

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

Clearly it appears the need to formulate an accurate element with an initial geometrical imperfection, which is adopted in various national design codes as an initial curvature of which the maximal magnitude is commonly taken as 0.1% of the member length. In this frame, Lui and Chen [2] have proposed a more rigorous cubic lateral displacement function expressed with the usual cubic Hermite polynomial, combined with two new shape functions. These new shape functions represent respectively a symmetric first buckling mode and an asymmetric second buckling mode of a fixed ends beam. Investigating the same problem, Chan and Zhou [1] have proposed a more relevant fifth order polynomial function for the displacement field. This element, commonly known as the (PEP), has been adopted in a computer software NIDA (Nonlinear Integrated Design Analysis) [3]. However, to rely on recent works [4,5], it has been noted the limitation of the finite element approach in the geometrical nonlinear analysis because of the inaccuracy of the new proposed elements. In a recent paper, Teh Lip [4] suggested to use the well-documented cubic element in the nonlinear analysis and design of steel frames. According to the author, the concern expressed in recent papers regarding the inaccuracy of the cubic element is ‘‘largely misplaced’’, owing to the pitfalls that are more likely with the newly proposed elements. According to this discussion, the finite element approach, if it is used in the geometrical nonlinear analysis, needs a more precise and closed-form expression of the displacement field for a beam element under high axial load. Recently, Chan and Gu [5], have proposed an exact solution for an imperfect beam-column by using the stability functions approach. In this study the authors have used Thimoshenko’s theory [6], and extended this theory to take into account the effect of the initial curvature along the element length. In the present paper, a versatile and efficient displacement shape functions expressed with respect to nodal generalized coordinates are derived analytically from the differential equilibrium equation. These shape functions are used in the context of the finite element method to formulate the stiffness matrix and the equivalent load vector for any arbitrary transverse load [7]. These shape functions are further used to introduce the effect of the initial curvature in the global equivalent load vector. This combined approach allows to link up the accuracy of the exact stability methods to the simplicity of the finite element approach. These new shape functions offer significant practical advantages to overcome limitations arising from the cubic Hermite element and other formulations [8–10].

2. Element formulation Consider a beam segment (ij) with initial curvature given in Fig. 1. The equilibrium equation along the element length can be expressed in compressive case as

577

y

x

L i Vi

j Ui

Uj φi

Ni Mi

Vj

W0 (x)

Ti

W1 (x) WF (x)

Nj

φj Tj

Mj

Fig. 1. Beam-segment model.

W 1 ðxÞ;xx þ K 2 ðW 1 ðxÞ þ W 0 ðxÞÞ ¼ MðxÞ=EI

ð1Þ

where (,xx) indicate the second derivative operator; W1(x) is the deflection superimposed to the initially curved beam. W0(x) is the initial curvature given by the following relation: npx ð2Þ W 0 ðxÞ ¼ V 0 sin ¼ N 0V 0 L The parameter n is associated with the shape of the initial curvature (if n = 1 then initial curvature has one maximal magnitude). V0 is the maximal magnitude of the initial curvature (if n = 1 then V0 is at mid-span). rffiffiffiffiffiffiffi jN j K¼ EI  x MðxÞ ¼ M i þ M j  M i L

ð3Þ ð4Þ

Mi and Mj are the nodal moments, N is the axial force. E = Young modulus of elasticity, I = moment of inertia. The general solution of Eq. (1) is W 1 ðxÞ ¼ C 1 cosðKxÞ þ C 2 sinðKxÞ þ C 3 x þ C 4 2 np  ðKLÞ x V0 þ sin L ðnpÞ2  ðKLÞ2

ð5Þ

C1, C2, C3, C4, are constants depending on the boundary conditions of the element (ij). We have to work with four degrees of freedom (d.o.f.). Two displacements Vi and Vj according to OY axis; and two rotations /i and /j around OZ axis (?x0y). By using the following arbitrary boundary conditions: (

x ¼ 0:

W 1 ð0Þ ¼ V i ;

W 1;x ð0Þ ¼ /i

x ¼ L: W 1 ðLÞ ¼ V j ; W 1;x ðLÞ ¼ /j (Subscript (,x) indicate the first derivative operator.)

ð6Þ

578

2

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

We obtain:

1 0 0 6 0 K 1 6 6 4 cosðKLÞ sinðKLÞ L K sinðKLÞ K cosðKLÞ 1 9 8 vi > > > > > > = < /i  /0 ¼ vj > > > > > > ; : /j  /0 cosðn  pÞ

98 9 1> > C1 > > > > > >> = 0 < C2 = 1> > > C3 > >> > ; : > ;> C4 0

The shape function associated with the effect of the initial curvature is given by N 1 ¼ B1 cosðbX Þ þ B2 sinðbX Þ þ B3 ðbX Þ  B1 þ ð7Þ

B1 ¼

ð8Þ

/0 is expressed by b2

V0 /0 ¼ np 2 2 L ðnpÞ  b

ð9Þ

The determinant of the system of Eq. (7) is det ¼ K  D

ð10Þ

D is expressed by D ¼ 2 þ 2 cosðbÞ þ b sinðbÞ

ð11Þ

The first nonzero value that vanishes this determinant is b = 2p, which corresponds to the critical load of a fixed– fixed beam-column given by the following expression: 2

N cr ¼

ð2pÞ EI L2

ð12Þ

The solution of the system of Eq. (7) gives the constants C1, C2, C3, C4. Then, by setting: X ¼ x=L

ð13Þ

Substitution of C1, C2, C3 and C4 into Eq. (5) yields the desired shape functions. The ones associated with the nodal d.o.f. (Vi/iVj/j) are given by 8 N vi ¼ A1 cosðbX Þ þ A4 sinðbX Þ  A4 ðbX Þ þ A5 > > > < N ¼ ðA cosðbX Þ þ A sinðbX Þ þ A ðbX Þ  A ÞL=b /i 2 5 1 2 > N ¼ A cosðbX Þ  A sinðbX Þ þ A ðbX Þ þ A vj 1 4 4 1 > > : N /j ¼ ðA3 cosðbX Þ  A1 sinðbX Þ þ A1 ðbX Þ  A3 ÞL=b ð14Þ where A1, A2, A3, A4 and A5 are expressed by 8 ðcosðbÞ  1Þ > > A ¼ > > 1 D > > > > ð b cosðbÞ  sinðbÞÞ > > A2 ¼ > > > D > < ðsinðbÞ  bÞ A3 ¼ > D > > > > sinðbÞ > > A4 ¼ > > > D > > > > : A ¼ ð1 þ cosðbÞ þ b sinðbÞÞ 5 D

2

ðnpÞ  b2

sinðnpX Þ

ð16Þ

where B1, B2, B3 are respectively expressed by

One poses: b ¼ KL

b2

1 ðb=npÞ ððsinðbÞ  b cosðbÞÞ D 1  ðb=npÞ2

þ cosðnpÞðb  sinðbÞÞÞ 1 ðb=npÞ ðð1  cosðbÞ  b sinðbÞÞ B2 ¼ D 1  ðb=npÞ2

ð17Þ

þ cosðnpÞðcosðbÞ  1ÞÞ

ð18Þ

1 ðb=npÞ ðð1  cosðbÞÞ B3 ¼ D 1  ðb=npÞ2 þ cosðnpÞð1  sinðbÞÞÞ

ð19Þ

Finally we can express the lateral displacement associated with the nodal displacements and rotations (Vi/iVj/j) which corresponds to the situation of a straight beam element, by the following relation: T

W ðxÞ ¼ fN vi N /i N vj N /j gfV i /i V j /j g ¼ fN i gfd i g

T

ð20Þ

Finally, the solution of the differential equilibrium equation Eq. (1) writes as T

W 1 ðxÞ ¼ W ðxÞ þ N 1 V 0 ¼ fN i gfd i g þ N 1 V 0

ð21Þ

N1V0 yields for the lateral displacement associated with the maximal magnitude V0 of the initial curvature. We note that the effect of the initial curvature is decoupled from the solution of the straight beam element. Now, superimposing the deflection W1 to the initial curvature W0 we obtain the final deflection of the element, written as W F ðxÞ ¼ W 1 ðxÞ þ W 0 ðxÞ ¼ W ðxÞ þ ðN 1 þ N 0 ÞV 0

ð22Þ

Let us introduce W v ðxÞ ¼ ðN 1 þ N 0 ÞV 0 ¼ N v V 0

ð23Þ

with: N V ¼ B1 cosðbX Þ þ B2 sinðbX Þ þ B3 ðbX Þ  B1 þ

1 1  ðb=npÞ

2

sinðnpX Þ

ð24Þ

Finally, the expression of the total deflection becomes: W F ðxÞ ¼ W ðxÞ þ W v ðxÞ ð15Þ

ð25Þ

The shape functions Ni (Eq. (20)), bring out exclusively the coupling between the axial force N and the deflection of the beam element. We note that each function Ni depends on the dimensionless parameter b. By allowing b to approach zero in Eq. (14), we obtain the well-known cubic Hermite polynomials given by Eq. (26).

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

(

  N vi ¼ 1  3X 2 þ 2X 3   N /i ¼ X  2X 2 þ X 3 L

  N vj ¼ 3X 2  2X 3   N /j ¼ X 2 þ X 3 L

ð26Þ

This provides a good check of the accuracy of the computations since Eq. (14) can readily be reduced to their corresponding cubic Hermite polynomials when the compressive force and then the parameter b vanish. For any arbitrary transverse load we have to use the Ni shape functions to derive the analytic expression of the equivalent nodal forces. The shape function Nv associated with the maximal magnitude V0 is used to derive the analytic expression for the equivalent nodal forces relating to the initial curvature. By a similar approach presented above, we can establish, the appropriate shape functions Ni and Nv for tensile case which are detailed in Appendix A. 3. Finite element procedure The aim of the following step, consists in establishing the coefficients of the stiffness matrix and also to express the equivalent load vector of a beam element. A fiber at distance z above the neutral axis of bending has the strain [11]: 1 ex ¼ U ðxÞ;x  z  ðW ðxÞ;xx Þ þ ðW F ðxÞ;x Þ2 2

ð27Þ

(Subscript (,x) and (,xx) indicate respectively the first and second derivative operator.) (U(x),x) represents axial stretching. z Æ (W(x),xx) is the strain produced by the curvature W(x),xx. (WF(x),x)2/2 is the strain produced by the lengthening associated with rotation of differential length dx through a small angle (WF(x),x) without x-direction motion of any point. The element strain energy is Z Z 1 P¼ Ee2 dA dx ð28Þ 2 L A x Neglecting terms of fourth order, and considering the following results: Z Z Z 2 dA ¼ A z dA ¼ I z dA ¼ 0 EAU ðxÞ;x  N A

A

A

ð29Þ Eq. (28) writes: Z Z EA EI 2 2 P¼ ðU ðxÞ;x Þ dx þ ðW ðxÞ;xx Þ dx L 2 L 2 Z N 2 ðW F ðxÞ;x Þ dx þ 2 L

ð30Þ

The first integral yields for truss element associated with d.o.f. ui and uj in x-direction. The second integral yields for a classical beam element. The third integral represent

579

the strain energy stored, when differential elements dx are 2 stretched an amount ðW F ðxÞ;x Þ dx=2. Let the strain energy P, written as P ¼ PN þ PF þ PV þ P0 Replacing WF(x) in the Eq. (30), we obtain: Z EA 2 ðU ðxÞ;x Þ dx PN ¼ L 2 Z Z EI N ðW ðxÞ;xx Þ2 dx þ ðW ðxÞ;x Þ2 dx PF ¼ 2 2 L ZL PV ¼ N ðW ðxÞ;x WvðxÞ;x Þdx ZL N 2 ðWvðxÞ;x Þ dx P0 ¼ 2 L

ð31Þ

ð32Þ ð33Þ ð34Þ ð35Þ

The stationary principle states that equilibrium prevails when the generalized coordinates {di} defines a configuration such that dP = 0 for any admissible {ddi}. This is possible only if ðdPN þ dPF þ dPV þ dP0 Þ ¼ 0

ð36Þ

The axial displacement U(x) and its first derivative U(x),x are given by  

ui

ui U ðxÞ ¼ N ui N uj ; U ðxÞ;x ¼ N ui N uj ;x uj uj ð37Þ with: Lx x N uj ¼ ð38Þ L L For the lateral deflection, let us suggest the following notation: N ui ¼

The generalized strain vector:  W ;xx feg ¼ ¼ ½Bfd i g W ;x

ð39Þ

The generalized stress vector: frg ¼ ½D½Bfd i g The displacement matrix: " # fN vi N /i N vj N /j g;xx ½B ¼ fN vi N /i N vj N /j g;x The constitutive matrix: " # 1 0 2 ½D ¼ EI 0  bL2

ð40Þ

ð41Þ

ð42Þ

The minus sign of the term (b/L)2 in the Eq. (42) rises from the compressive axial force N considered in this study. By using the strain energy variation dPF, we obtain the coefficients of the secant stiffness matrix of a straight beam element under compressive axial forces Z BT DB dx ð43Þ k¼ L

580

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

The explicit stiffness matrix k, can be expressed as 3 2 12EI 6EI 12EI 6EI X D X 7 6 L3 D L2 L3 L2 7 6 7 6 6EI 4EI 6EI 2EI 6 H U 7 X 7 6 L2 X 2 L L L 7 k¼6 7 6 12EI 6EI 12EI 6EI 6 X D X7 7 6 L3 D 2 3 2 L L L 7 6 5 4 6EI 2EI 6EI 4EI U H X X 2 2 L L L L

ð44Þ

where 8 > > >
b2 cosðbÞ  bsinðbÞ b sinðbÞ  b2 ; U¼ 2ð2 þ 2cosðbÞ þ b sinðbÞÞ 4ð2 þ 2 cosðbÞ þ b sinðbÞÞ 3 > b sinðbÞ b2 cosðbÞ  b2 > > ; X¼ :D ¼ 6ð2 þ 2cosðbÞ þ b sinðbÞÞ 12ð2 þ 2cosðbÞ þ b sinðbÞÞ

The vector Fv of the virtual equivalent nodal loads relating to the initial imperfection is derived from dPV. This is done by 9 0 8 1 N vi;x > > > > > > Z > @ A L L > N vj;x > > > ; : N /j;x Thus, the initial curvature is taken into account as an equivalent nodal load, which simplify its implementation in any existing static and dynamic nonlinear geometrical analysis software. P0 is independent of the nodal d.o.f. (di) so, dP0 vanishes. 4. Numerical analysis

ð45Þ

4.1. Shape functions H, U, D, and X are the well-known stability functions, here obtained by using the finite element approach. This result provides a good check for the accuracy of the present computations. Furthermore, the equivalent nodal loads for any arbitrarily distributed load Q(x) can readily be evaluated from the shape functions previously expressed. This is obtained from the following integration, which must be performed along the distance length where the distributed load is defined Z q¼ N T QðxÞdx ð46Þ L

In the case of a distributed moment m(x), Eq. (46) is modified and expressed in the form: (47) Z dN T m¼ mðxÞdx ð47Þ L dx

The shape functions previously formulated (Eq. (14)) should be examined numerically in terms of amplification of the deflection versus parameter b. By setting b = 5 which corresponds to the situation of a high axial load, we can see in Figs. 2 and 3 that the shape functions Nvi and Nvj stay close to their corresponding cubic polynomials N vi and N vj (Eq. (26)) which relate only the first order effect corresponding to b = 0. However, in Figs. 4 and 5 we can observe, always for b = 5, that the shift between N/i, N/j and their corresponding cubic polynomials are more pronounced, indicating that the nodal rotations affect the deflection more than the nodal displacement. Figs. 6 and 7 show respectively a comparison between the shape functions N/i and N/j and their corresponding (PEP) proposed by Chan and Zhou [1]. These both shape functions correspond to a unit rotation at the left end for N/i (respectively right end for N/j) of the beam-segment

Fig. 2. Shape function Nvi for b = 0 and b = 5 (compr. and tensile cases).

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

581

Fig. 3. Shape function Nvj for b = 0 and b = 5 (compr. and tensile cases).

Fig. 4. Shape function N/i for b = 0 and b = 5 (compr. and tensile cases).

(ij). We can observe that the (PEP) solution underestimates the deflection in comparison with the present solution. We can observe that the shift between these two solutions becomes more important when b increases mainly for p < b < 2p. Indeed, we can see from the numerical results that for p 6 b 6 1.5p the error is less than 10% between the proposed shape function N/i and their corresponding (PEP) whilst for 1.5p < b < 2p the error increases rapidly and becomes very high for b = 2p. 4.2. Equivalent nodal forces 4.2.1. Case of a concentrated load q In the situation of a concentrated load q at the midspan, we can note that the nodal moment varies with

respect to the parameter b, whereas the equivalent nodal force is constant and equal to q/2. Fig. 8 shows the variation of the equivalent nodal moment versus b indicating that the moment tends to 0.125qL when b is allowed to vanish (b = 0). On the other hand, when b tends to 2p, the value of the moment increases considerably.

4.2.2. Case of a distributed load Q In the situation of a uniformly distributed load Q, we can observe that the nodal moment varies with respect to the parameter b, whereas the equivalent nodal force is constant and equal to QL/2. Fig. 9 shows the variation of the nodal moment versus b indicating that the moment tends to (QL/12) when b is allowed to vanish (b = 0). On the

582

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

Fig. 5. Shape function N/j for b = 0 and b = 5 (compr. and tensile cases).

Fig. 6. Comparison with Chan solution: shape function N/i for b = 0, b = 4, b = 5.

Fig. 7. Comparison with Chan solution: shape function N/j for b = 0, b = 4, b = 5.

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

583

Fig. 8. Equivalent nodal moment versus b for fixed–fixed beam supporting a concentrated load q at the mid-span.

Fig. 9. Equivalent nodal moment versus b for fixed–fixed beam supporting a uniform distributed load Q.

other hand, when b tends to 2p, the value of the moment increases considerably. 5. Numerical applications The element previously formulated is implemented into computer software based on a nonlinear algorithm resolution. A step-by-step loading is adopted as a way to proceed with an iterative substitution method. Example 1. Study of an imperfect column To evaluate the accuracy of the proposed element in terms of static stability, the case of the imperfect column with geometrical and mechanical characteristics shown in Fig. 10 should be examined. The geometrical imperfection of the column is included as an initial lateral curvature, which is approached by a

half-sine curve with magnitude equals to 0.001 and 0.002 of the member length. The column is considered under a concentrated load P at the end B. The relationship between lateral deflection WF at halfheight and axial load P is shown in Fig. 11. The plots shown in Fig. 11 are respectively obtained in the cases of pinned ends and fixed ends by using only one element per member with no more than two iterations for each load step. In Fig. 11 we can readily appreciate that the predicted load for the pinned ends case is about 68 N which is close to Euler’s buckling load (72 N) given by the relation p2EI/ L2, and is about 273 N (4 · 72 N) for the fixed ends case. The CPU time does not exceed 0.22 s for both examples using Pentium(R) 4 CPU 3.00 GHz. On the other hand, Fig. 11 reveals that an increase in the magnitude of the initial imperfection does not lead to any changes in the value of the predicted load.

584

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

P

P B

L

B

P

A

C

B

V0

L

V0

L

L= 1 m S= 6.4516*10-4 m2 I = 3.4684*10-8 m4 E= 200*106 N/m2 V0= .001L and .002L

L= 1 m S= 6.4516*10-4 m2 I = 3.4684*10-8 m4 E= 200*106 N/m2 V0= .001L

V0

A A

Pinned ends

Fixed ends

Fig. 10. Geometrical and mechanical characteristics of column AB.

Pinned end

Fixed end

Fig. 12. Geometrical and mechanical characteristics of framed structure ABC.

4 V0 =.001L

3 End A: pinned

3

End A: fixed

2

Case of pinned ends

1

Case of fixed ends

Lateral deflection WF (cm) at half-height of the column AB

Lateral deflection WF (cm) at half-height of the column AB

V0 =.002L

2

1

0 0

50

100

150

200

250

300

Load P(N) Fig. 11. Buckling analysis of column AB.

Example 2. Study of a framed structure This example is chosen to emphasize on the ability of the present element when it is incorporated into a framed structure (see Fig. 12), which is analysed numerically by using only one element per member. On the other hand, this example is shown to appreciate the sensitivity of the predicted load of this structure. The chosen structure with geometrical and mechanical characteristics illustrated in Fig. 12 is examined in two different boundary condition situations. The first case corresponds to pinned end at A and fixed end at C and the second case corresponds to fixed ends at both A and C. Geometrical imperfection with magnitude equals to 0.001 of the member length is included in the column AB, the element BC is assumed to be initially straight. A vertical concentrated load P is applied at the node B. The relationship between lateral deflection WF at halfheight of the column AB and vertical load P is shown in Fig. 13. The plots are obtained by using only one element per member with no more than two iterations for each load

0 0

50

100

150 Load P(N)

200

250

300

Fig. 13. Buckling analysis of framed structure ABC.

step. When the node A is assumed to be pinned, we can readily appreciate that the predicted load of the element AB is about 100 N which is greater than the Euler’s buckling load (72 N). When the node A is fixed, the predicted load of the element AB is about 200 N, which is smaller than 288 N of a fixed ends column. These results can be explained by the flexural rigidity of the node B. We can see with these examples that the predicted load of the portal frame is mainly sensitive to the boundary condition at the node A. The CPU time does not exceed 0.25 s for both examples using Pentium(R) 4 CPU 3.00 GHz. 6. Conclusion This paper proposes new shape functions of an imperfect beam element under high axial load, for the geometrical nonlinear analysis. These shape functions are exact, since they are obtained analytically from the resolution of

R. Adman, H. Afra / Advances in Engineering Software 38 (2007) 576–585

the differential equilibrium equation. The exactness of these shape functions is in the same sense as the cubic shape functions are exact for a classical beam element. This result offers to the finite element approach the needful accuracy required besides the simplicity of the method. It should also be noted that the equivalent nodal loads for any arbitrarily distributed load become reliable, which is highly desirable for the numerical convergence and to predict the buckling load accurately. Thus, it must be pointed out that Euler’s buckling load is checked by using only one element per member, with few iterations for each load step. Furthermore, these shape functions are advantageous for the analysis of braced and unbraced structures without special numerical treatments [12–15]. The expressions of these shape functions seem to be complex, but in fact, once they have been incorporated into finite element routines of frames analysis software, the second-order analysis and the nonlinear integrated design become reliable and easy to use for practical purpose. Appendix A. Tensile case The displacement’s shape functions for the tensile case are expressed as below: 8 N vi ¼ A4 sinhðbX Þ  A1 coshðbX Þ  A4 ðbX Þ þ A5 > > > > < N /i ¼ ðA5 sinhðbX Þ þ A2 coshðbX Þ  A1 ðbX Þ  A2 ÞL=b > N vj ¼ A4 sinhðbX Þ þ A1 coshðbX Þ þ A4 ðbX Þ  A1 > > > : N /j ¼ ðA1 sinhðbX Þ þ A3 coshðbX Þ  A1 ðbX Þ  A3 ÞL=b ð49Þ with: 8 > > A1 > > > > > > > > A2 > > > > < A3 > > > > > > > A4 > > > > > > > :A 5

ðcoshðbÞ  1Þ D ðb coshðbÞ þ sinhðbÞÞ ¼ D ð sinhðbÞ þ bÞ ¼ D sinhðbÞ ¼ D ð1  coshðbÞ þ b sinhðbÞÞ ¼ D ¼

where rffiffiffiffiffiffiffi jN j b¼ L EI

ðN is the tensile loadÞ

ð50Þ

ð51Þ

585

The stability functions are then expressed by 8 > > >
b2 coshðbÞ  b sinhðbÞ ; 4ð2  2 coshðbÞ þ b sinhðbÞÞ

> > > :D ¼

b3 sinhðbÞ b2 coshðbÞ  b2 ; X¼ 12ð2  2 coshðbÞ þ b sinhðbÞÞ 6ð2  2 coshðbÞ þ b sinhðbÞÞ

U¼

b sinhðbÞ  b2 2ð2  2 coshðbÞ þ b sinhðbÞÞ

ð52Þ The elastic matrix D is modified and then become expressed in the form: " # 1 0 2 ð53Þ ½D ¼ EI 0 bL2 References [1] Chan SL, Zhou ZH. On the development of a robust element for second-order ‘non-linear integrated design and analysis (NIDA)’. J Construct Steel Res 1998;47():169–90. [2] Lui EM, Chen W-F. Analysis and behaviour of flexibly-jointed frames. Eng Struct 1986;8 (April). [3] NIDA—nonlinear integrated design and analysis for framed structures. NAF series software, User’s manual, March 2000. [4] Teh Lip H. Cubic beam elements in practical analysis and design of steel frames. Eng Struct 2001;23:1243–55. [5] Chan S-L, Gu Jian-Xin. Exact tangent stiffness for imperfect beamcolumn members. J Struct Eng 2000(9):1094–107. [6] Timoshenko SP. ‘The´orie de la stabilite´ e´lastique’. 2nd ed. 1966: Dunod; 1966. [7] Admanan R, Afra H. Analytical solution of a beam element for elastic buckling analysis. In: Topping BHV, Mota Soares CA, editors. Proceedings of the seventh international conference on computational structures technology. Stirling (United Kingdom): Civil-Comp Press; 2004, paper 249. [8] Kuo-Mo Hsiao, Horng Horng-Jann, Chen Yeh-Reh. A corotational procedure that handles large rotations of spatial beam structures. Comput Struct 1987;27(6):769–81. [9] Kondoh K, Tanaka K, Atluri SN. An explicit expression for the tangent-stiffness of a finitely deformed 3-D beam and its use in the analysis of space frames. Comput Struct 1986;24(2):253–71. [10] Fujii Fumio. A simple mixed formulation for elastica problems. Comput Struct 1983;17(1):79–88. [11] Petersen C. Statik und Stabilitat der Baukonstruktionen. Braunschweig/Wiesbaden; 1982, Seite 855ff. [12] Ph Jetteur, Cescotto S, de Ville de Goyet V, Frey F. Improved nonlinear finite elements oriented bodies using an extension of marguerre’s theory. Comput Struct 1987;17(1):129–37. [13] Chen WF. Structural stability: from theory to practice. Eng Struct 2000;22:116–22. [14] Chan SL, Shu Gan-Pin, Zhi-Tao Lu¨. Stability analysis and parametric study of pre-stressed stayed columns. Eng Struct 2002;24:115–24. [15] Kalaga S, Adluri SM. Beam columns with finite deflections. J Struct Eng 2000(2):266–9.