A co-Rotational beam element for geometrically nonlinear analysis of plane frames

Vietnam Journal of Mechanics, VAST, Vol. 35, No. 1 (2013), pp. 51 – 65 A CO-ROTATIONAL BEAM ELEMENT FOR GEOMETRICALLY NONLINEAR ANALYSIS OF PLANE FRAMES Nguyen Dinh Kien1, Trinh Thanh Huong1, Le Thi Ha2 1Institute of Mechanics, VAST, 18 Hoang Quoc Viet, Hanoi, Vietnam 2Hanoi University of Transport and Communication, Vietnam Abstract. A co-rotational beam element for geometrically nonlinear analysis of plane frames is presented. On the base of the shallow arch expression for local strain, the ele- ment is formulated by using exact polynomials to interpolate the transversal displacement and rotation. Using the formulated element, the critical load and equilibrium path are computed with the aid of the bracketing procedure and the arc-length control method, respectively. Numerical examples show that the proposed element is capable to give ac- curate results with a smaller number of elements comparing to the elements previously used in the examples. The effect of the nonlinear term used in the local strain expression on the numerical results is also investigated and highlighted. Keywords: Plane frame, geometric nonlinearity, co-rotational beam element, critical load, equilibrium path. 1. INTRODUCTION Frame structure is widely used in civil engineering, and nonlinear analysis of this structure is of interest from both academic and practical points of views. The essential for a nonlinear analysis is to have a nonlinear element which is capable to model accurately nonlinear behavior of the structure. Development of nonlinear elements, in general, and beam elements, in particular, is an important topic in the field of structural mechanics. Many beam elements which included geometrical nonlinearity for the large displacement analysis of frames can be found in the literature, and some of them have been documented in various well-known textbooks, e.g. [1, 2]. Depending on the choice of reference configu- ration, the nonlinear beam elements can be classified into three main types, namely the total Lagrangian formulation, the updated Lagrangian formulation and the co-rotational formulation. In the co-rotational formulation, which will be adopted in the present work, a local coordinate system attached to the element is introduced, and the finite element formulation is firstly formulated in the local coordinates and then transferred to the global system by the aid of transformation matrices [2, 3]. Literature survey shows that the co-rotational beam elements can be formulated by using different local strain definitions, beam theories 52 Nguyen Dinh Kien, Trinh Thanh Huong, Le Thi Ha as well as the way to define a local coordinate system. In addition to the plane Bernoulli element using a linear strain definition described by Crisfield in [2], the formulations proposed by Pacoste and Eriksson [3], Hsiao and co-workers [4], Meek and Xue [5], are some main co-rotational elements available in the literature for geometrically nonlinear analysis of plane elastic beam and frame structures. In a recent work, Jafari and co- workers [6] also adopted an element attached coordinate system for construction of the flexibility matrix for geometrically nonlinear Timoshenko beam element, which can be viewed as a form of the co-rotational formulation. It is necessary to note that the total Lagrangian formulation and the co-rotational formulation, as demonstrated in [3], though adopted different reference systems for the kinematic description, they are identical. In the previous work, the first author and his co-worker have conducted some re- search on large displacement analysis of elastic and elasto-plastic frames by developing the co-rotational Bernoulli beam elements [7, 8]. The effect of shear deformation on the large displacement response of frames has also investigated by the first author by formulation of the co-rotational Timoshenko element [9]. In the context of the finite element analy- sis, due to the independence of the transversal displacement and rotation, a Timoshenko beam element for geometrically nonlinear analysis of frames, as discussed by Pacoste and Eriksson in [3] or by the first author in [9, 10], can be formulated by linearly interpolat- ing all the kinematic variables, provided that some techniques such as the reduced-order integration are employed to prevent the formulated element from the shear locking. An alternative way for preventing a shear deformable beam element from the shear locking problem, as discussed by Reddy in [11], is to employ appropriate interpolation functions for the kinematic variables. The exact cubic and quadratic polynomials obtained by solv- ing the homogenous equilibrium equations of a beam element proposed by Luo in [12], which will be adopted in the present work, are examples of such interpolation functions. The present paper formulates a co-rotational element for geometrically nonlinear analysis of plane frame structure. The element is adopted a shallow arch expression for local strain and employing exact polynomials obtained from the field consistent approach to interpolate the transversal displacement and rotation. Using the formulated element, nonlinear equilibrium equations in the finite element context are constructed and solved by an incremental and iterative method procedure. The bracketing procedure and the arc- length control method are adopted in computing the critical loads and tracing equilibrium paths, respectively. A number of numerical examples are given to demonstrate the accuracy and efficiency of the formulated element. The effect of the nonlinear term in the local strain expression on the accuracy of numerical results is also investigated and highlighted. 2. CO-ROTATIONAL METHOD As above mentioned, in the co-rotational method an element attached coordinate system which continuously translates and rotates with the element during its deformation process is introduced. By using such a local system, the element deformation is decom- posed into rigid body and pure deformation parts, and the geometric nonlinearity induced by large rigid body motion is incorporated into transformation matrices which can be derived from the relationship between the local and global quantities [2,3]. The kinematic A co-rotational beam element for geometrically nonlinear analysis of plane frames 53 relationship and detail expressions for the transformation matrices are described in this section. Fig. 1. A plane beam element and its kinematics in local and global coordinate systems Fig. 1 shows a two-node plane beam element and its kinematics in the local and global coordinate systems. The element is assumed initially straight and the beam material is linearly elastic. Global coordinates (X,Z) are fixed in space. A local system (x, z) is attached to the element, and it continuously translates and rotates with the element. Denoting (Xi, Zi) and (Xj, Zj) are respectively the coordinates of node i and node j with respective to the global system. The local system (x, z) is chosen as well as its original is always at node i, and x axis directs to node j. By choosing such a local system, the axial displacement at node i and transversal displacements at both the nodes in the local system are always zero. Thus, vectors of the nodal displacements in the local and global systems are given by, respectively, dl = {ul θli θlj} T (1) and d = {ui wi θi uj wj θj} T (2) Here and afterwards a superscript T denotes the transpose of a vector or a matrix, and subscript “l” stands for the “local”. The local displacements in Eq. (1) can be computed as [2, 3] ul = ln − l , θli = θi − θr , θlj = θj − θr (3) where l, ln and θr are the initial length, current length and rigid rotation of the element, respectively. These quantities can be computed by geometrical consideration as l = √ l2X + l 2 Z , ln = √ (lX + d41)2 + (lZ + d52)2 θr = acrtan [ lXd52 − lZd41 lX(lX + d41) + lZ(lZ + d52) ] (4) with lX = Xj −Xi, lZ = Zj −Zi are the projections of undeformed element on the global axes X and Z, respectively; d41 = uj − ui, d52 = wj − wi. 54 Nguyen Dinh Kien, Trinh Thanh Huong, Le Thi Ha Denoting f inl and fin are the element internal force vectors corresponding to the nodal displacement vectors defined by Eq. (1) and Eq. (2), respectively. These force vectors contain the axial, shear forces and moments at the nodes finl = {Nl Mli Mlj} T (5) fin = {Ni Qi Mi Nj Qj Mj} T (6) Assuming the strain energy expression for the element U has been derived in terms of the nodal displacements, the global internal force vector can be computed by differentiation of this strain energy (which is invariant of the coordinate systems) with respect to the global displacement vector as fin = ∂U ∂d = ∂U ∂dl ∂dl ∂d = TT1 ∂U ∂dl = TT1 finl (7) where T1 is the transformation matrix, which can be computed from Eq. (3) as T1 = ∂dl ∂d (8) The element tangent stiffness matrix in the global system kt can be computed as kt = ∂2U ∂d2 = ∂fin ∂d = TT1 ktlT1 +NlT2 +MliT3 +MljT4 (9) where ktl = ∂ 2U/∂d2l is the local tangent stiffness matrix; T2, T3 and T4 are the trans- formation matrices, defined from Eq. (3) as T2 = ∂2ul ∂d2 , T3 = ∂2θli ∂d2 , T4 = ∂2θlj ∂d2 (10) From Eq. (3), one gets T3 = T4 = Tr = − ∂2θr ∂d2 (11) with θr is defined by Eq. (4). Thus, Eqs. (7)-(11) completely define the element formulation provided that the local internal force vector finl and tangent stiffness matrix ktl are known. 3. LOCAL FORCE VECTOR AND STIFFNESS MATRIX This section derives the local internal force vector and tangent stiffness matrix for the element. Adopting the Timoshenko beam theory, the displacements in the local coordinate system are given by u(x, z) = u0(x)− zθ(x) , w(x, z) = w(x) , v(x, z) = 0 (12) where u0(x) and w(x) are the mid-plane axial and transversal displacements; θ(x) is the cross sectional rotation, independent of the transversal displacement w(x). The shallow arch expression, a degenerated form of the Green strain, deduced from Eq. (12), is adopted for the local strain as x = ∂u0 ∂x + 1 2 ( ∂w ∂x )2 − z ∂θ ∂x = m + zχ , γxz = ∂w ∂x − θ (13) A co-rotational beam element for geometrically nonlinear analysis of plane frames 55 where χ = −∂θ/∂x is the beam curvature. The strain energy for the beam element asso- ciated with the strain components defined by Eq. (13) has the form U = 1 2 ∫ l 0 EA2mdx+ 1 2 ∫ l 0 EIχ2dx+ 1 2 ∫ l 0 κGAγ2dx (14) where EA, EI and GA are the axial, bending and shear rigidities, respectively; κ is the shear correction factor, depends on the cross section geometry. In order to derive the local internal force vector and tangent stiffness matrix from the strain energy expression, we need to express the train energy in term of the local nodal displacements. To this end, interpolation functions are necessary to introduce for the axial and transverse displacements and the cross section rotation. In the present work, the field consistence approach discussed by Luo [12] is adopted, in which the interpolation functions for the transversal displacement w(x) and the rotation θ(x) are found by solving the homogeneous equilibrium equations of the beam element  κGA ( − ∂θ ∂x + ∂2w ∂x2 ) = 0 EI ∂2θ ∂x2 + κGA ( −θ + ∂w ∂x ) = 0 (15) with the boundary conditions w|x=0 = wli , θ|x=0 = θli , w|x=l = wlj , θ|x=l = θlj (16) By introducing a dimensionless parameter φ = 12 l2 EI κGA (17) and noticing that wli = wlj = 0, Eqs. (15) and (16) give the solutions w = Nw1θli +Nw2θlj θ = Nθ1θli +Nθ2θlj (18) where Nw1 = l 1 + φ [(x l )3 − ( 2 + φ 2 )(x l )2 + ( 1 + φ 2 ) x l ] Nw2 = l 1 + φ [(x l )3 − ( 1− φ 2 )(x l )2 − φ 2 x l ] (19) and Nθ1 = 1 1 + φ [ 3 (x l )2 − (4 + φ) x l + (1 + φ) ] Nθ2 = 1 1 + φ [ 3 (x l )2 − (2− φ) x l ] (20) The axial displacement is still linearly interpolated as u0 = x l ul, and thus to prevent the formulated element from the membrane locking, which might be occurred from the 56 Nguyen Dinh Kien, Trinh Thanh Huong, Le Thi Ha unbalance interpolation, the membrane strain m in Eq. (13) is necessary to replace by the effective strain as [2] eff. = 1 l ∫ l 0 mdx = 1 l ∫ l 0 [ ∂u0 ∂x + 1 2 ( ∂w ∂x )2] dx (21) Using Eqs. (19) and (20) we can write the effective strain in terms of the local nodal displacements as eff. = 1 l ul + 1 (1 + φ)2 [ φ(2 + φ) 24 (θli − θlj) 2 + 1 30 (2θ2li − θliθlj + θ 2 lj) ] (22) Substituting Eqs. (19), (20) into Eq. (14), and using Eq. (22) one can express the strain energy in term of the local nodal displacements as U = lEA 2 { 1 l ul + 1 (1 + φ)2 [ φ(2 + φ) 24 (θli − θlj) 2 + 1 30 (2θ2li − θliθlj + 2θ 2 lj) ]}2 + EI 2l(1 + φ)2 [ φ(2 + φ)(θli − θlj) 2 + 4(θ2li + θliθlj + θ 2 lj) ] + lφ2κGA 8(1 + φ)2 (θli + θlj) 2 (23) Expressions for the local internal force vector finl and tangent stiffness matrix ktl are then obtained by successive differentiation of the strain energy expression, Eq. (23), with respective to the local degrees of freedom as finl = ∂U ∂dl and ktl = ∂2U ∂d2l (24) From Eqs. (23) and (24), one can write the components of the local internal force vector, finl , in the explicit forms as Nl = EAeff . Mli = lEA (1 + φ)2 [ φ(2 + φ) 12 (θli − θlj) + 1 30 (4θli − θlj) ] eff . + EI l(1 + φ)2 [φ(2 + φ)(θli − θlj) + 2(2θli + θlj)] + lφ2κGA 4(1 + φ)2 (θli + θlj) Mlj = lEA (1 + φ)2 [ φ(2 + φ) 12 (θlj − θli) + 1 30 (4θlj − θli) ] eff . + EI l(1 + φ)2 [φ(2 + φ)(θlj − θli) + 2(θli + 2θlj)] + lφ2κGA 4(1 + φ)2 (θli + θlj) (25) where eff . is given by Eq. (22). The local tangent stiffness matrix, ktl , can be written in the form ktl = ka + kb + ks (26) in which kb = EI l(1 + φ)2   0 0 0 0 (φ2 + 2φ+ 4) −(φ2 + 2φ− 2) 0 −(φ2 + 2φ− 2) (φ2 + 2φ+ 4)   , ks = lφ2κGA 4(1 + φ)2   0 0 0 0 1 1 0 1 1   (27) A co-rotational beam element for geometrically nonlinear analysis of plane frames 57 are the constant stiffness matrices stemming from the beam bending and the shear de- formation, respectively; ka is the nonlinear stiffness matrix resulted from the membrane strain, and having the form ka = EAl   1 l2 1 l ti 1 l tj t2i + (5φ2 + 10φ+ 8) 60(1 + φ)2 eff . titj − (5φ2 + 10φ+ 2) 60(1 + φ)2 eff . sym. t2j + (5φ2 + 10φ+ 8) 60(1 + φ)2 eff .   (28) where ti = ∂eff . ∂θli = 1 60(1 + φ)2 [5φ(2 + φ)(θli − θlj) + 2(4θli − θlj)] tj = ∂eff . ∂θlj = 1 60(1 + φ)2 [5φ(2 + φ)(θlj − θli) + 2(4θlj − θli)] (29) The element formulated from the strain energy defined by Eq. (23) is based on the nonlinear shallow arch expression for the local axial strain, Eq. (13). However, a linear expression for the local strain can be assumed for the local strain [2, 3], and the element based on the linear strain definition can be obtained by deleting the terms in the Eq. (23) corresponding to the nonlinear term 1 2 (∂w/∂x)2. In order to investigate the effect of the nonlinear term on the numerical results, we denote B2CL and B2CS as the elements based on the linear and shallow arch expressions of the local strain, respectively. It is necessary to note that when φ → 0, the strain energy given by Eq. (23) deduces to the strain energy expression of the shallow arch based co-rotational Bernoulli beam element, previously derived by Pacoste and Eriksson in [3] (conf. Eq. (44) of Ref. [3]). Thus, with the definition of the coordinate systems as that in [3], when φ → 0 the Timoshenko beam element of the present work deduces exactly to the 2D co-rotational Bernoulli element in [3]. In this regards, the element of the present work is able to model the slender beams and frames. In other work, the Timoshenko beam element formulated in present paper is free of the shear locking. 4. NUMERICAL PROCEDURES The equilibrium equations for non-linear analysis of structures are obtained by set- ting the out-of-balance forces to zeros [2] g(p, λ) = qin(p)− λfef = 0 (30) where the out-of-balance force vector g is a function of the current structural nodal dis- placements p, and the load-level parameter λ; qin is the structural internal force vector, constructed from the element force vector fin in standard way of the finite element method; fef is the fixed external loading vector. The non-linear equations (30) can be solved by an incremental/iterative technique based on the Newton-Raphson method, in which the norm of vector g(p, λ) is guided towards zero. 58 Nguyen Dinh Kien, Trinh Thanh Huong, Le Thi Ha In order to trace complete equilibrium paths for possible complex behavior such as snap-through and snap-back, the arc-length method developed by Crisfield [2, 13] is employed. In the method, the load parameter λ becomes a variable, and a constraint equation is introduced as a = ( ∆pT∆p+∆λ2ψ2fTef fef ) −∆l2 = 0 (31) where ∆l is a fixed value chosen by the analyst; the vector ∆p and scalar ∆λ are the incremental displacement and load parameter, defined from the last converged state; ψ is the scaling parameter. With ψ = 1 , Eq. (31) defines a spherical surface with a radius ∆l in the load-displacement space, and Eq. (31) define a cylindrical surface when ψ = 0. The methods corresponding to ψ = 1 and ψ = 0 are known in the literature as the “spherical arc-length method” and “cylindrical arc-length method”, respectively. Thus, with the introduction of the constrain equation, instead of solving the equilibrium equations (30) directly, we find the intersection of the equilibrium path with the surface defined by Eq. (31). The iterative procedure is obtained from the truncated Taylor expansion of the vector g and a around the current equilibrium point as g = g|0 + ∂g ∂p |0δp+ ∂g ∂λ |0δλ = g|0 +Kt|0δp+ fefδλ = 0 a = a0 + 2∆p T δp+ 2∆λδλψ2fTef fef = 0 (32) where g|0 and a|0 are defined from the current equilibrium point. Eq. (32) gives the iterative displacements δp and iterative load parameter δλ as{ δp δλ } = − [ Kt|0 −fef 2∆pT 2∆λψ2fT ef fef ] −1 { g|0 a|0 } (33) As we see the augmented Jacobian matrix within the square brackets in Eq. (33) remains non-singular even when theKt is singular. However, the Jacobian matrix in (33) is neither symmetric nor banded, and the computation based on the arc-length method is expensive. The cylindrical arc-length method is adopted in the present work. More detail on the arc-length method and its numerical implementation are described in [2]. In addition, the bracketing procedure [14, 15] is adopted in computing critical and limit loads of the structures. To this end, a load-control strategy is employed in solving the equilibrium equation (30), and the lowest eigenvalue τ of the structural tangent stiffness matrix Kt, assembled from the element matrix kt, is employed as indicator for predicting a critical or a limit point. When passing the lowest critical or limit point, predicted by the negative sign of τ , the bracketing procedure is active to compute an intermediate point λi, by using an interpolation scheme for the control parameter λ as λi = λl − λr − λl τr − τl τl (34) where λl, λr are the values of the load-level parameter at the equilibrium points before and after passing the critical or limit point, respectively; τl and τr are the lowest eigenvalues of Kt at the points. The bracketing procedure is terminated when a convergency criteria A co-rotational beam element for geometrically nonlinear analysis of plane frames 59 is satisfied. The following criterion is adopted herewith(∣∣∣∣λr − λlλi τ¯ ∣∣∣∣ )1/2 < β , with τ¯ = τi (|τl τr|)1/2 (35) where β is the tolerance factor, taken by 10−4 in the present work. The obtained value λi from the bracketing procedure is taken as the critical load-level parameter. 5. NUMERICAL EXAMPLES In this section, various beam and frame structures shown in Fig. 2 are analyzed to verify the formulated element. A tolerance factor of 10−4 is also chosen for judging convergency of the iterative numerical procedure described in section 4. The shear modulus is taken as G = E/2 for all the problem, and if not mentioned otherwise, the shear correction factor κ = 5/6 is chosen for all the analyses reported below. 5.1. Euler buckling load The cantilever beam under a centric axial force P with geometric and material data shown in Fig. 2(a) adopted fromRef. [3] is employed to study the convergency and accuracy of the formulated element in computing critical loads. The buckling coefficient of the beam, fB = PcL 2/EI with Pc is the critical load, computed by various numbers of the B2CL and B2CS elements is listed in Tab. 1. The analytical solution [16], Pc = Pe/(1+Pe/κGA) with Pe = pi 2/4EI/L2, gives the buckling coefficient fB = 2.4674. As seen from Tab. 1, the accuracy of the critical load obtained by the B2CS can be accepted at a very coarse mesh. The accuracy of the critical load computed by the B2CL element is steady improved by using more elements in the analysis, but the load computed by twenty B2CL elements just gives the same accuracy as that obtained by two B2CS elements. Thus, it is hardly improved the accuracy of the critical load by using a large number of elements when the nonlinear term is ignored from the local strain expression. Fig. 2. Test problems for verifying formulated element Following the work by Pacoste and Eriksson in [3], the critical load of the cantilever beam is computed by different slenderness ratios, namely s = 4 and s = 100, where 60 Nguyen Dinh Kien, Trinh Thanh Huong, Le Thi Ha s = L/ √ I/A. Keeping all the data as in Fig. 2(a), the computation is performed with the beam cross-sectional area A = 2.88× 10−5 and A = 0.018, respectively. Table 1. Buckling coefficient for cantilever beam computed by different element numbers Number of element Element B2CL B2CS 1 3.0000 2.4860 2 2.5966 2.4687 4 2.4993 2.4675 6 2.4815 2.4674 8 2.4753 2.4674 10 2.4725 2.4674 20 2.4687 2.4674 Table 2. Critical load for cantilever beam at different slenderness ratios Present Ref. [3] Ref. [17] s B2CL B2CS b2dtt b2dcs 4 Pc 550.7782 548.1241 553.0 548.2 548.64 uc 19.1242 19.0321 14.18 19.03 100 Pc 444.9728 444.0597 446.8 444.6 444.13 uc 0.0247 0.0247 0.0248 0.02477 In consistent with the work in [3], the beam is discretized by 8 elements. The critical load and the corresponding axial displacement at free end of the beam obtained by the B2CL and B2CS elements are listed in Tab. 2. For comparison, the critical loads obtained from a numerical method [3], and from a finite strain theory [17] are also listed in the table. In the table, the ‘b2dtt’ and ‘b2dcs’ are the Timoshenko and Bernoulli beam el- ements, formulated in [3] by using the total Lagrangian and co-rotational formulations, respectively. The numerical result of the present work is very good in agreement with the result computed by the ‘b2dcs’ of Ref. [3], and with that of [17]. As has been shown in [3], the ‘b2dcs’ is very accurate in computing the critical load and large displacement response of planar beam and frame structures. It is noted that for the case s = 4, a shear correction factor ψ = 1200 as adopted in [3] is chosen in the computation to inhibit the shear effect. 5.2. Cantilever beam under transversal point load The cantilever beam subjected to a transversal point load as shown in Fig. 2(b) is employed herewith to test the ability of the elements in tracing the equilibrium path for large displacement. The problem has been analytically analyzed by Mattiasson [18], and it has been employed by several authors to test the accuracy of nonlinear beam elements. A co-rotational beam element for geometrically nonlinear analysis of plane frames 61 The load-displacement curves computed by only two elements formulated in the present work are shown in Fig. 3, where the analytical solution derived by Mattiasson [18] and the finite element result computed by five consistent Bernoulli beam elements by Nanakorn and Vu [19] are also depicted. As seen from the figure, the numerical result computed by 0 0.2 0.4 0.6 0.8 1 0 2 4 6 8 10 Normalized displacements, u/L & w/L No rm al ize d loa d, P L2 /E I Mattiasson Nanakorn & Vu 2 B2CL elements u/L w/L 2 B2CS elements Fig. 3. Load-displacement relation cantilever beam under transversal point load only two B2CS elements of the present work agrees very well with the analytical solution in [18]. The curves computed by the B2CL element slightly differs from that obtained by the B2CS element. The difference in the numerical result computed by the B2CL and B2CS elements in this example can be eliminated by using 4 or more elements. 5.3. Cantilever beam under eccentric axial point load The cantilever beam with the data shown in Fig. 2(c) under an eccentric axial point load is considered herewith. The problem has been previously analyzed by Wood and Zienkiewicz by using five paralinear elements [20]. The load-displacement curves obtained by using a mesh of three elements formulated in the present work are depicted in Fig. 4. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 1 2 3 4 Normalized displacements, u/L & w/L No rm al ize d lo ad , P L2 /pi 2 E I w/L u/L Wood & Zienkiewicz 3 B2CL elements 3 B2CS elements Fig. 4. Load-displacement relation for cantilever beam under eccentric axial load 62 Nguyen Dinh Kien, Trinh Thanh Huong, Le Thi Ha The good agreement between the result computed by the B2CS element of the present work with the numerical result obtained in [20] is observed. Some difference between the curves computed by the B2CL element with that obtained by the B2CS element is seen in the figure. This difference can be diminished by using more elements in the analysis, and with a mesh of 5 elements, no difference is observed. 5.4. Williams toggle frame The toggle frame in Fig. 2(d) was previously employed by several authors to test the behavior of nonlinear finite elements in large displacement analysis of frames [3,4,19]. The frame data are: L = 12.943,H = 0.386, b = 0.753, h = 0.243, E = 10.3×106. The previous work shows snap-through behavior of the frame under the downwards point load P . The load-deflection curves for the frame computed by various numbers of the B2CL and B2CS elements are depicted in Fig. 5. The numerical result obtained by using ten paralinear 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 10 20 30 40 50 60 70 Transverse displacement, w (in.) Ap pli ed lo ad , P (lb .) Wood & Zienkiewicz 2 B2CL elements 2 B2CS elements 8 B2CL elements Fig. 5. Load-displacement curve for Williams frame Table 3. Limit load around snap-through situation of Williams toggle frame Number of element Element B2CL B2CS 2 39.9904 33.0079 4 39.9253 32.9903 6 35.8876 32.8587 8 34.4337 32.8409 10 33.9811 32.8270 12 33.6169 32.8066 elements of Wood and Zienkiewicz [20] is also shown in the figure by small triangles. Very good agrement between the curve computed by two B2CS element of the present work and A co-rotational beam element for geometrically nonlinear analysis of plane frames 63 that in [20] is noted. The two B2CL elements gives the curve, which largely differs from the corresponding curve obtained by the two B2CS and the numerical result of Ref. [20]. To improve the accuracy, a larger number of the B2CL is required, and as seen from Fig. 5 the curve computed by eight B2CL elements has been improved considerately. To examine the effect of the nonlinear term 1 2 ( ∂w ∂x )2, the limit load around the snap-through situation of the frame is computed with different numbers of the the B2CL and B2CS elements by using the bracketing procedure. The obtained result listed Tab. 3 shows that the limit load is overestimated by using the B2CL element, even with a fine mesh of twelve B2CL elements. 5.5. Lee’s frame An asymmetric frame under a concentrated load shown in Fig. 2(e), which also known in the literature as the Lee’s frame, is investigated. The problem has been previously analyzed by Pacoste and Eriksson [3], Hsiao and Huo [4], and Cichón [21]. The frame shows highly nonlinear with snap-through and snap-back behavior. The frame data are L = 120 cm, A = 6 cm2, I = 2 cm4, E = 720 t/cm2, where L, A, I and E are the length, cross section area, moment of inertia and elastic modulus, respectively. A model of five elements with different length as shown in Fig. 2(e) is employed in the present work. The displacements in both horizontal and vertical directions at the loaded point are computed by using the B2CS element. 0 20 40 60 80 100 −1.5 −1 0 1 2 Displacements at loaded point (cm) Ap pl ie d lo ad (to n) u w 5 B2CL elements Hsiao & Huo Fig. 6. Load-displacement relation of Lee’s frame Fig. 6 shows the relation between the applied load and the displacements at loaded point of the frame. The computed curves of the present work agrees well with the numerical result obtained by using 10 elements in Ref. [4], which shown in the figure by small circles. As in case of the Williams toggle frame, the curves computed by the B2CL (not shown in the figure) differ from the curves in Fig. 6. To improve the numerical result, more elements are required when the B2CL element is used. In order to show the ability in modelling the shear deformation, the computation was performed with various values of the slenderness ratio, namely s = 150, 50, and 30 by 64 Nguyen Dinh Kien, Trinh Thanh Huong, Le Thi Ha 0 20 40 60 80 100 −1 0 1 2 Displacements at loaded point (cm) Ap pli ed lo ad (to n) u w s = 150 s = 60 s = 30 Fig. 7. Load-displacement relation of Lee’s frame with different slenderness ratios using five B2CS elements. Keeping the rest data as above, the computation is performed with A = 3.125, 0.3472 and 0.1250 cm2, respectively. The computed curves depicted in Fig. 7 show a reduction in the limit load which the frame can withstand when reducing the slenderness ratio. Thus, the shear deformation effect on the large displacement behavior of the frame can be well modelled by the formulated element. 6. CONCLUSIONS In this paper, a co-rotational beam element for geometrically nonlinear analysis of plane frames has been formulated. The shallow arch expression was adopted for the local strain and the exact polynomials obtained from the field consistence approach were employed as interpolation functions for the transversal displacement and sectional rotation. Using the formulated element, the critical loads and equilibrium paths are computed with the aid of the bracketing procedure and the arc-length method. A number of examples have been given to demonstrate the accuracy and efficiency of the formulated element. The numerical results have shown that the proposed element is capable to give accurate results by using a smaller number of the elements comparing to the elements previously used in the examples. The effect of nonlinear term in the local strain expression on the critical loads and large displacement behavior of the structure has also been numerically investigated. It has been shown that a larger number of elements is required to ensure the accuracy when the nonlinear term in the expression of the local strain is ignored. In addition, ignoring the nonlinear term leads to an overestimation of the critical and limit loads, and in some cases this shortcoming is hardly improved by increasing number of the elements used in analysis. ACKNOWLEDGEMENT The financial support from Vietnam NAFOSTED to the present work is gratefully acknowledged. A co-rotational beam element for geometrically nonlinear analysis of plane frames 65 REFERENCES [1] T. Belytschko, W.K. Liu and B. Moran, Nonlinear finite elements for continua and structures (John Wiley & Sons, Chichester, 2000). [2] M.A. Crisfield, Non-linear finite element analysis of solids and structures, Vol. 1: Essentials (John Wiley & Sons, Chichester, 1991). [3] C. Pacoste and A. Eriksson, Beam elements in instability problems, Comput. Method. Appl. Mech. Eng., 144, (1997), 163-197. [4] K.M. Hsiao and F.Y. Huo, Nonlinear finite element analysis of elastic frames, Comput. Struct., 26, (1987), 693-701. [5] J.L. Meek and Q. Xue, A study on the instability problem for 2D-frames, Comput. Method. Appl. Mech. Eng., 136 (1996), 347-361. [6] V. Jafari, S.H. Vahdani and M. Rahimian, Derivation of the consistent flexibility mtrix for geometrically nonlinear Timoshenko frame finite element, Finite Elem. Anal. Des., 46, (2010), 1077-1085. [7] N.D.Kien and D.Q. Quang, Large deflections analysis of frames by elements contaning higer- order terms, Vietnam Journal of Mechanics, 25, (2003), 243-254. [8] N.D.Kien and D.Q. Quang, Beam element for large displacement analysis of elasto-plastic frames, Vietnam Journal of Mechanics, 26, (2004), 39-54. [9] N.D.Kien, Effect of shear defromation on large deflection behavior of elastic frames, Vietnam Journal of Mechanics, 26, (2004) 167-181. [10] D.K. Nguyen, Post-buckling behavior of beam on two-parameter elastic foundation, Int. J. Struct. Stabil. Dynam., 4, (2004), 21-43. [11] J.N. Reddy, On locking-free shear deformable beam finite elements, Comput. Method. Appl. Mech., Eng., 149, (1997), 113-132. [12] Y.H. Luo, Explanation and elimination of shear locking and membrane locking with field consistence approach, Comput. Method. Appl. Mech. Eng., 162, (1998), 249-269. [13] M.A. Crisfield, A fast incremental/iterative solution procedure that handles ‘snap-through’, Comput. Struct., 13, (1981), 55-62. [14] M.A. Crisfield, Non-linear finite element analysis of solids and structures. Volume 2: Advanced topics, John Wiley & Sons, Chichester, (1997). [15] J. Shi, Computing critical points and secondary paths in nonlinear structural stability analysis by the finite element method, Comput. Struct., 58, (1996), 203-220. [16] S.P. Timoshenko and J.M. Gere, Theory of elastic stability (McGraw-Hill, New York, 1961). [17] Y. Goto, T. Yoshimitsu and M. Obara, Elliptic integrals solution of plane elastica with axial and shear deformations, Int. J. Solids Struct., 26, (1990), 375-390. [18] K. Mattiasson, Numerical results for large deflection beam and frame problems analysed by means of elliptic integrals, Int. J. Numer. Methods Eng., 17, (1981), 145-153. [19] P. Nanakorn and L.N. Vu, A 2D field-consistent beam element for large displacement analysis using the total Lagrangian formulation, Finite Elem. Anal. Des., 42, (2006), 1240-1247. [20] R.D. Wood and O.C. Zienkiewicz, Geometrically nonlinear finite element analysis of beams, frames, arches and axisymmetric shells, Comput. Struct., 7, (1977), 725-735. [21] C. Cichón, Large displacement in-plane analysis of elastic-plastic frames, Comput. Struct., 19, (1984), 737-745. Received November 14, 2011 VIETNAM ACADEMY OF SCIENCE AND TECHNOLOGY VIETNAM JOURNAL OF MECHANICS VOLUME 35, N. 1, 2013 CONTENTS Pages 1. Dao Huy Bich, Nguyen Xuan Nguyen, Hoang Van Tung, Postbuckling of functionally graded cylindrical shells based on improved Donnell equations. 1 2. Bui Thi Hien, Tran Ich Thinh, Nguyen Manh Cuong, Numerical analysis of free vibration of cross-ply thick laminated composite cylindrical shells by continuous element method. 17 3. Tran Ich Thinh, Bui Van Binh, Tran Minh Tu, Static and dynamic analyses of stiffended folded laminate composite plate. 31 4. Nguyen Dinh Kien, Trinh Thanh Huong, Le Thi Ha, A co-rotational beam element for geometrically nonlinear analysis of plane frames. 51 5. Nguyen Chien Thang, Qian Xudong, Ton That Hoang Lan, Fatigue perfor- mance of tubular X-joints: Numberical investigation. 67 6. Hoang H. Truong, Chien H. Thai, H. Nguyen-Xuan, Isogeometric analysis of two–dimensional piezoelectric structures. 79 7. Pham Chi Vinh, Do Xuan Tung, Explicit homogenized equations of the piezo- electricity theory in a two-dimensional domain with a very rough interface of comb-type. 93