- Email: [email protected]
A corotational procedure that handles large rotations of spatial beam structures
Cm$~ers
6 S~~~rures Vol. 27. No. 6. pp. 769-781.
1987 0
Printed in Great Britain.
aa-7949187 13.00 + oxa 1987 Pcrgamon Journals Lid.
A COROTATIONAL PROCEDURE THAT HANDLES LARGE ROTATIONS OF SPATIAL BEAM STRUCTURES Kuo-MO HSIAO,HORNG-JANNHORNG and YEH-REN CHEN ~~rtment
of Mechanic& Engineering, National Chiao Tung University, Hsinchu. Taiwan, Republic of China (Received 5 January 1987)
Abstract-A practical motion process of the three dimensional beam element is presented to remove the restriction of small rotations between two successive increments for large displacement and large rotation analysis of space frames using incremental-iterative methods. In order to improve convergence properties of the equilibrium iteration, an n-cycle iteration scheme is introduced. The nonlinear formulation is based on the corotational formulation. The transformation of the element coordinate system is assumed to be accomplished by a translation and two sueesssive rigid body rotations: a transverse rotation followed by an axial rotation. The element formulation is derived based on the small deflection beam theory with the inclusion of the effect of axial force in the element coordinate system. The membrane strain along the deformed beam axis obtained from the elongation of the arc length of the beam element is assumed to be constant. The element internal nodal forces are cakulated using the total defo~atioM1 nodal rotations. Two methods, referred to as direct method and incremental method, are proposed in this paper to calculate the total defo~~ona~ rotations. An incremental-iterative method based on the Newto~Rap~n method combined with arc length control is adopted. Numerical studies are presented to demonstrate the accuracy and efficiency of the present method.
1. INTRODUCTION
The development of new and efficient formulations for the nonlinear analysis of beam structures has attracted the study of many researchers in recent years, and different alternative formulation strategies and procedures to accomodate large rotation capability during the large defo~ation process have been presented [I-lo]. These formulations can be divided into three categories: Total Lagrangian (TL) formulation, Updated Lagrangian (UL) formulation and Corotational (CR) formulation. It should be noted that within the corotating system either a TL or a UL formulation, or even a formulation based on a small deflection theory may be employed. The large number of publications on the nonlinear analysis of beam structures is, at least partially, due to the fact that various kinematic nonlinear formulations can be employed. It seems that large rotations in plane frames pose no major problem. Hsiao and Hou [IO] introduced a simple and effective corotational formulation of beam element and numerical procedure, which can remove the restriction of small rotations between two successive increments for the large displacement analysis of plane frames using incremental-iterative methods. Unfortunately, the method presented in [lo] cannot be applied to three dimensional frames. The difficulty of obtaining effective solutions is particularly pronounced in the analysis of spatial beam structures; a general three dimensional nonlinear formulation is not a simple extension of a two dimensional fo~ulation, because large rotations in three dimensional analysis are not true vector
quantities; that is, they do not comply with the rules of vector operations and the result will in general depend on the order in whch the rotations are taken. This point has been throughly discussed by Argyris [1 I] and Wempner [ 121. The problem of large rotations on space structures has received wide attention in the ~terature; many different strategies based on the TL, the UL, or the CR formulations have been reported, those of [Z-S, 6, 1l-241 being only a small fraction of the total. In 1191Hughes and Liu developed a specialized shell element which can handle arbitrarily large rotations. Argyris has covered the subject of corotational coordinates extensively including a lengthy discourse on the subject of large rotations [20,21]. Belytschko er al. [3,4, 14 have applied corotational formulation to the dynamic analysis of space frames where arbitrarily large rotations can be expected. Horrigmoe and Bergan (IS] have su~~fully applied a corotational approach to their shell elements Rankin and Brogan [23] have introduced a corotational procedure which may enable existing shell element formulations to be used in problems that contain arbitrarily large rotations. Recently Hsiao 1241 has proposed a motion process for triangular shell elements to remove the restriction of small rotations between two successive increments for nonlinear shell analysis using incremental-iterative methods. The wide range of numerical examples studied in [3 4,7,9, 11 , 13, 15, 16,18,20,‘21,23-261 indicates that the corotation approach, first described by Argyris et al. [25], may be very useful in the analysis of spatial structures containing arbitrarily large
770
Kuo-MO HSIAOet al.
rotations. However, most strategies based on the corotational formulation suffer from one inherent drawback: they are restricted to small rotations between two successive load increments during the deformation process. This limitation arises because the incremental nodal rotations are considered to be vector quantities. Although the method introduced in [24] may remove this restriction for triangular shell elements, unfortunately, this method cannot be applied to the space beam elements, because, unlike the shell elements, the element coordinate of the space beam elements cannot be determined using only nodal coordinates. The objective in this paper is to present a practical motion process of the three dimensional beam element which can remove the restriction of small rotations between two successive increments for large displacement and large rotation analysis of space frames using incremental-iterative methods. Here the motion process presented in [24] is modified to accomodate the characteristics of the motion of beam elements. In this paper the transformation between the element coordinate systems is assumed to be accomplished by a translation and two successive finite rotations: a lateral rotation about an axis perpendicular to the current beam axis followed by an axial rotation about the current beam axis. Two methods, termed direct method and incremental method, are introduced to describe the motion process of the beam element and to determine the total deformational nodal rotations. The dominant factors in the geometrical nonlinearities of space structures are attributable to finite rotations, the strains remaining small. For space frames discretized by finite elements, this implies that the motion of the individual elements with proper size will, to a large extent, consist of rigid body motion.
If the rigid body motion part is eliminated from the total displacements, the deformational part of the motion is always a small quantity relative to the local element axes; thus, incorporated with the corotational formulation, the small deflection beam theory with the inclusion of the effect of axial force is adopted here to deal with the large rotations but small strains problems. The numerical algorithm used here is an incremental-iterative method based on the NewtonRaphson method combined with constant arc length of the incremental displacement vector [27,28]. In order to improve the convergence properties of the equilibrium iterations, an n-cycle iteration scheme, which is an extension of the two-cycle iteration scheme proposed in [lo], is introduced. Numerical examples are presented and compared with the results reported in the literature to demonstrate the accuracy and efficiency of the proposed method. 2. COORDINATESYSTEMS One of the basic considerations in formulating nonlinear structural problems is the selection of a description of motion. In the present study, the corotational (CR) formulation is adopted. In this formulation, each element is associated with a local Cartesian element coordinate system 2, (i = 1,2, 3) as shown in Fig. 1, that rotates and translates with the element but does not deform with the element. The element coordinate system used herein is a right handed one and is defined as follows. The origin of the element coordinate system is located at the local node 1, the 2, axis passes through the centroids of member end sections, and the Z2 and 5, axes are parallel to the principal directions of the undeformed end cross section.
Fig. 1. Coordinate systems, member deformations and associated forces.
Corotational
771
procedure for rotations of beam structures
Also shown in Fig. 1 is a fixed global coordinate system xi (i = 1,2,3) used to define the location of the nodal points. We note that the equilibrium equations of a structure are written in the fixed global coordinate system, and the incremental nodal parameters of the system of equations are also calculated in this coordinate system. The element equations are first formulated in the element coordinate system, and then transformed to the global coordinate system using standard procedure [29] prior to the element assemblage process. If we consider a vector B (or B if measured in the _Zicoordinate system) with global components Bi (i = 1,2,3), and element coordinate components Si (i = 1,2,3), we have the following transformation:
P !!
9
R’
n
Fig. 2. Finite rotation of vector.
respectively; 4 is the angle of counterclockwise rotation, and n is the unit vector along the axis of rotation. 3. CR-FORMULATION OF BEAM ELEMENT
where clti= cos(xi, Xj). The beam element employed here has two nodes with six degrees of freedom per node (Fig. 1): these are the translations 0, in the fi (i = 1,2,3) directions at nodes j (j = 1,2), and the rotations flj about the Zi axes at nodes j (j = 1,2). The global nodal parameters for the system of equations associated with the individual elements are chosen to be the translations Vii in the xi (i = 1,2,3) directions at nodes j (j = 1,2), and rotations 0, about the xi (i = 1,2,3) axes at nodes j (j = 1,2). In this study an incremental-iterative method is used to solve the nonlinear equilibrium equations. Both the incremental nodal translations and the incremental nodal rotations are regarded as vector quantities. Thus, using eqn (l), the nodal vectors, AUj= {AU,j, AU,, AUjj) and AtIj = {A8ij, A&, Ae,}, referred to the global coordinate system can be transformed to AUj = {AO,j, AUq, AOJj} and A4 = {A&, A&, Ae,}, referred to the element coordinate system, respectively. It should be noted that Aej cannot be interpreted as component rotations about Cartesian axes &. In this study, Ae, = {Ae,, 0,0}, the components of Ah along the 2, axis, and Afj”,= (0, A&, AB,), the components of A4 perpendicular to the R, axis, are considered to be rotation vectors to define rotations, details of which will be discussed later. For convenience of the later discussion, the term ‘rotation vector’ is used to represent a finite rotation. Figure 2 shows a vector R which as a result of the application of a rotation vector 4n is transported to a new position R’. The relation between R and R’ may be expressed as [30] R=cos+R+(l
The formulation of the beam element developed here is applicable to arbitrarily large rotations but restricted to small rotations relative to the element axis. The beam element is formulated in the element coordinate system based on the small deflection beam theory with the inclusion of the effect of axial force. The element, as shown in Fig. 1, has two nodes with six degrees of freedom per node, and can transmit an axial force, two shear forces, two bending moments and a torque. Herein the beam element is assumed to be straight and of constant cross section. The cross section is doubly symmetric, thus excluding coupling of the torsional stiffness to that of bending and axial stiffness. Shearing deformations and warping effects are neglected. The material is assumed to be linearly elastic. The element stiffness matrix is obtained by superimposing its bending, geometric, torsional and axial stiffness matrices. The element internal nodal forces are evaluated using the total deformations. 3.1. Kinematics of beam element It is assumed that the lateral deflection curves of the beam member are the cubic Hermitian polynomials in the & and .?s directions of the element coordinates (the principal directions of the undeformed end cross section), and that the axial rotation varies linearly along the member. The membrane strain along the deformed element axis is assumed to be constant. Thus, the membrane strain can be evaluated from the elongation of the arc length of the member. The lateral deflection curves of the beam member may be given by
-cos4)(n.R)n + sin 4(n x R),
(2)
where * and x denote the dot and the cross product
where a bar over a quantity denotes that it is defined in the element coordinate system. P and Ware lateral
IWO-MOHstao et al stant membrane strain along the deformed beam axis, the membrane strain of the beam axis can be written as Cl = (S - SJ)/&,
“2
/
undeformed
end section
Fig. 3. Deformational nodal rotations.
deflections in the & and i, directions, respectively. iii, and ej (i = 2,3; j = f , 2) are the nodal displacements and rotations shown in Fig. 1. N, (i = i,4) are shape functions and are given by N, = l/4(1 - #(2
+ t)
N2 = c/8(1 - r’)(l - r)
(6)
where S, = L is the initial arc length of the beam axis. Figure 3 shows that the normals of the undeformed element end sections at nodes j (j = 1,2), 6, are rotated to g,,,, the current deformed normals of the element end sections by the rotation vectors e)l,, which is ~~ndicular to the f,axis of the element coordinate system. The representations of the lateral deformational nodal rotations are based on the assumptions that these rotations are small. On the basis of this assumption e,,, the 5, (i = 2,3) components of the rotation vectors 6,,j are chosen to be the Iateral defo~ational nodal rotations about Zi axes at nodes j. Note that the direction of the undeformed normal of the element end sections coincides with the positive direction of the 5, axis. Thus, the second subscript j of gq is omitted throughout this paper. 3.2. Determination of element coordinate system and element de~ormat~o~al rotations
Assume that the incremental-iterative method is used for the solution of nonlinear equilibrium equaNl = c/8(1 + c*)(l + r), (4) tions and the equilibrium configuration of the Ith increment is known. Let AUj and A@,(j = 1,2) be where c = Zi, - %ir is the current chord length of the incremental nodal displacement and rotation the beam member, and f, are the 2, coordinates of vectors of an element at nodes j extracted from the incremental nodal parameters of the system of equanodes j (j = I, 2) in the element coordinate system; tions. At this point, an interesting and relevant 5 = - I + 23,/c is a nondimensional coordinate. question arises. Given the incremental nodal disNote that in this study, the relative displacements placements and rotations, how are the current of the elements are referred to the element coodinate system, Due to the definition of the element co- element coordinate system, the axial and lateral ordinate system, the lateral nodal displa~ments 041 defo~ationa1 nodal rotations dete~ined? Let ‘xj (j = 1,2) denote the node coordinate (i = 2,3) at the nodal points j (j = 1,2) are identical vectors of an element in its Ith equilibrium conto zero. The nonzero deformational nodal displacefiguration; the current node coordinate vectors xj ments of an element can be divided into the axial are obtained by adding the incremental nodal disrelative displacement, the axial relative rotation and lateral deformational rotations. The element defor- placement vectors AUj, so that mation can be decomposed into the membrane deformation, the torsional deformation and the flexural xj = ‘x, + AU,. (7) deformation. Herein the flexural deformation is deThe 3, axis of the current element coordinate termined by the lateral deformational nodal rotations system can then be constructed using x, given in using elementary beam theory and the torsional eqn (7) and the definition of the element coordinate deformation is determined from the axial relative nodal rotation, while the membrane deformation is system. But, unlike the cases of triangular shell obtained from the change of arc length of the beam elements fl5, 18,241, the .?r and 5 axes cannot be determined using only the node coordinates. The axis which can be calculated from the lateral determination of the current element coordinate deflection. If the arc length of the beam axis is expressed by system will be discussed in the process of element motion. 1 For determining the element deformational nodal (1 + p’2 -i- @‘2)“2 d&, S=c/2 (5) rotation and element rigid body rotations, we proi -1 pose two methods to describe the process of the element motion in this paper. In the first method, where c is the current chord length of the beam member, ( )’ denotes x,-derivatives and P and rii are referred to as the direct method, the lateral deformational nodal rotations are determined from the given in eqn (3), then from the assumption of conN,= l/4(1 +S)‘(Z--r)
Corotational
773
procedure for rotations of beam structures
(b)
(d)
Fig. 4. Process of element motion.
orientations of the undeformed and deformed normals of the element end sections. In the second method, referred to as incremental method, the total lateral deformational nodal rotations are calculated by incrementation. For both methods, the total twist nodal rotations are determined by incrementation. The processes of element motions and the methods corresponding to these motion processes to determine the deformational nodal rotations and element coordinate system are described as follows. (a) Direct mefhod. The process of element motion is divided into the following six steps.
1. A rigid body translation by AU,. The whole element is translated by AU,, where AU, is the incremental nodal displacement vector of node 1. The origin of the ‘3, axes is translated to the origin of the Zi axes as shown in Fig. 4(a). 2. A lateral rigid body rotation by the rotation vector 8’. The rotation vector 6’ (referred to *ii coordinate system) is given by &‘=cos-‘(‘e,*e,)
42,x e, II‘el x e,II’
(8)
where ‘e, and e, are unit vectors associated with the
774
Kuo-MO
HSIAOef al.
‘.f, and x’, axes rcspcclivcly. The rotation vector oi’, passing through node 1, is applied to the whole element except ‘tQ, the deformed normals of the element end sections at nodes j (j = 1,2) at the Ith equilibrium configuration. Here it is assumed that ‘e,,, are not rotated by the rotation vector ci’, but are translated with the motion of nodes j (Figs 4(b) and (d)). The ‘& axes are rotated by an angle about the axis perpendicular to the ‘2, and f, axes (Fig. 4(a)), and the resultant coordinate system is labeled 2; axes (Fig. 4(c)). As can be seen, the 2; axis coincides with the 2, axis. 3. Finite rotations of ‘C,+by the rotation vectors Afig. The deformed normals ‘c?~are rotated to &, by the application of the rotation vectors A6:, as shown in Figs 4(c) and (d), where AtiLj are the components of Afij (j = 1,2) (the given incremental nodal rotation vectors referred to 2: coordinate system) perpendicular to the Xi axis. 4. Twist rotations by A4j. The rotation vectors A4j are applied to nodes j (j = 1,2), where AtJj are given by A& = A&, - /I
(9)
/? = l/2(6&, + A&,),
(10)
in which tI,, are the components of A& along the 5; axis as shown in Fig. 4(c). 5. A stretch by (c - ‘c)e,. Node 2 is translated along the 2, axis (Fig. 4(d)), where c and ‘c are chord lengths of the element corresponding to the current configuration and the equilibrium configuration of the Ith increment. 6. An axial rigid body rotation by g. The rotation vector (eqn (10)) passing through node 1 is applied to the whole element. The intermediate axes, ai, are rotated about the 2; axis by an angle 11 fi 1)to produce the & axes of the current element coordinate system as shown in Fig. 4(f). Vectors El and 5; shown in Fig. 4(d) are rotated by this rotation vector to reach their final positions Z, and &, (not shown in Fig. 4(f)). ‘I’he orientations of the current element coordinate axes fi may be obtained from the rigid body rotations caused by the rotation vectors given in steps 2 and 6 of the above process. The orientation of the deformed normals gdi and the undeformed normal 5, of the element are also determined by the above process of motion. Thus, the lateral deformational nodal rotation vectors, e,, referred to the current element coordinate system may be expressed as
ii, 11
0 e;I,= IT4 cos-‘(E:~~)
x *dj
II& x ~,$II’
0 e;.=
f$, = ‘4 + (A8, - a).
(13)
1103, For simplicity of computation, only eqn (13) is used to calculate the total lateral deformational nodal rotations for the numerical examples studied in this paper. It is believed that identical results will be obtained when eqn (11) is used.
(11)
O3j
The current axial deformational may be obtained by
whcrc ‘4, is the axial dcformational nodal rotation vector of the Zth equilibrium configuration at nodes j (j = 1,2); A$j is given in eqn (9). (b) Incremental method. In this method, the motion process is also divided into six steps. Only steps 2 and 3 of this process are different from those in the direct method, and are described below. 2. A rigid body rotation by o?‘.The rotation vector 6i’ (eqn (8)) passing through node 1 is applied to the whole element. As can be seen in Figs 4(b) and (e), the deformed normals of the Ith increment ‘gdl are rotated to an intermediate position I&;. 3. Deformational rotations by (68;, - a’). The intermediate deformed normals *EL are rotated to other intermediate positions 5; as shown in Fig. 4(e) by the application of the rotation vectors (A&, - o?‘), in which A#Aj are the components of AJj (j = 1,2) perpendicular to 2; axis, and OS’is the rotation vector given in eqn (8) (referred to 2: coordinate system). Vectors AgAj, oi’ and ‘6; (Fig. 4(e)) are rotated by the rotation vector [ (eqn (10)) to their final positions Afimj,6 and ‘4 (referred to the current element coordinate system &), which are not shown in Fig. 4(f). The element coordinate obtained from this process is the same as that constructed by the direct method. The orientations of the deformed and undeformed normals of the element end sections can be determined from this process of motion as well. Thus, the deformational nodal rotations can be calculated using eqn (11). Due to the assumption of small deformational nodal rotations, the vector operations might be valid for the deformational nodal rotation vectors. Thus, an alternative, referred to as incremental method, is introduced here. The concept of this method is similar to that in [l&24]. The total deformational rotations of the current configuration referred to the current element coordinate are obtained by adding the incremental nodal rotations (A&j-oS) to the deformational nodal rotations of the equilibrium configuration of the Ith increment, and are expressed as
3.3. Element st1~ne.s.smatrix
nodal rotations
The total element stiffness matrix is formulated by superimposing the bending stiffness matrix i& and geometric stiffness matrix $ of the basic beam element, and the axial stiffness matrix g,,, and the (12) torsional stiffness matrix g, of the linear bar element. The derivation of these matrices is well documented
115
Corotational procedure for rotations of beam structures in the textbooks and thus will not be repeated here. However, these matrices are given as follows. (a) Bending stiffness matrix K,:
where GJ is the torsional rigidity and L is the initial arc length of the beam axis. 3.4. Element nodal force vectors
(14) where r
12 -6L -12
-6L 4L= 6L
-12 6L 12
2~~
6L
12 6L -12
6L 4L’ -6L
-12 -6L 12
6L
2L2
L-6L
-6L 1 2L’ 6L
(1%
4~2J
and f
6L 1 2L2 -6L 9 (16)
The nodal force vectors of the elements corresponding to the global coordinate system are evaluated first in the current element coordinate system, and then transformed to the global coordinate system using standard procedure. Since small deformations are assumed, the element nodal forces can-in the element coordinate system-be evaluated in much the same way as in linear analysis. For linearly elastic material properties the element nodal force vectors can be calculated as follows. (a) Bending nodal force vector t,: &i=@,+K$,,
where L is the initial length of the beam axis, and EI, and E13are the flexural rigidities about the 4 axis and 2, axis respectively. The degrees of freedom corresponding to $ are ob=
{~,d2,,
&2.42*
~2,r&,,
(17)
022,42}r
where uU are nodal translations and gU are nodal rotations as shown in Fig. 1. (b) Geometric stiffness matrix &:
rt,= [
$2 0
0 rt,,
1
(18)
36 -3L
-3L 4L2
-36
3L
-36
-3L -L2
3L 36
3L
3L
4L2
4L2 -3L3L
-3L -3636
-3L3L
-L2
-3L
-L2
-3L and
K*&
3L [ -3636 3L
1 (19)
1
-1
[ -1
9 (20)
4L2 I
1
1’
(21)
where AE is the axial rigidity and L is the initial arc length of the beam axis. (d) Torsional stiffness matrix K,: z(
GJ
‘=L
[ 11’ 1
-1
-1
where ~b={F31rn221,F,2,11-3,,F2,,~,,,F22r1U32} is shown in Fig. 2; Kb is the bending stiffness matrix given in eqn (14); ab is the total bending deformation vector given in eqn (17). Note that due to the definition of the element coordinate system, the only nonzero elements in ab are f$, the deformational nodal rotations at nodes j (j = 1,2) about the fi (i = 2,3) axes, which may be obtained by using eqns (11) or (13). (b) Axial nodal force vector Fm;,: The element internal nodal forces are calculated by the total nodal deformation rotations. The axial nodal force vector F,,, = {I?,], F,2} (see Fig. 1) can be evaluated by introducing nodal virtual displacements 80, = 6 {D,, , u12} at nodes 1 and 2 in the f, direction, and equating the work done by the axial nodal force F,,, going through the virtual displacement JO,,, to the work done by the internal stress resultant T going through the virtual strain St,,, (that corresponds to the imposed virtual displacement) along the deformed beam axis as s
-L2
where L is the initial arc length of the beam axis and F,, is the axial nodal force at node 2. The degrees of freedom corresponding to KE are the same as that corresponding to izb. (c) Axial stiffness matrix KM: it_=:
(23)
4L2
-6L
(22)
cm,Fm =
I
Tsr, dS,
(24)
0
where S is the current arc length of the beam axis. The stress resultant T may be given by T = AEF-,,
(25)
where A is the cross section area, E is Young’s modulus and r,,, is the membrane strain of the current deformation. From eqn (6), the virtual strain a<,,, corresponding to ~50, can be expressed as
arm= as/s, (1 + p’2+ p2)‘nd{
So
(26)
Kuo-MO HSIAOelal.
776
in which 6c is the variation of the chord length of the beam axis with respect to SD,,,. Substituting eqns (S), (25), and (26) into eqn (24) gives
Since the virtual nodal displacements aa, are arbitrary, the axial nodal forces are obtained from eqn (27) as
Because the assumption of the small strain and small deformation, S2/S,c in eqn (28) is approximated by unity in this paper, and eqn (28) is thus reduced to
5. SOLUTION
ALGORITHM
An incremental-iterative method based on the Newton-Raphson method is adopted here. In order to deal with the limit points and snap through, the arc length of the incremental displacement vector is kept constant during the equilibrium iteration using Crisfield’s method [27,28]. An n-cycle iteration scheme is introduced here to improve the convergence characteristics of the equilibrium iteration. If the equiiibrium configuration of the Ith increment is assumed to be known, the system tangent stiffness matrix KT then can be calculated at this configuration and an initial displacement increment Aq for the next increment may be obtained by using Euler predictor as Aq = A%,
(33)
where Al is the initial incremental loading parameter and qT= Kf'P is the tangential displacement of unit loading P.For all increments other than the first, Ai. is obtained in much the same way as that mentioned in [27] and is given by
for numerical computation. (c) Torsional nodal force vector F,:
A1 = fAa(q;q,)“2, F,={;;;}=T{_;},
(30)
where 4 = ($,, - $,,) is the total relative rotation of the member about the 2, axis, r&j are the total twist rotations at nodes j (j = 1,2) about the fi axis and are obtained in eqn (12), and &?,j (j = 1,2) are twist moments shown in Fig. 1.
(34)
where the sign is chosen following an approach due to Bergan and Ssreide [32] in which the sign follows that of the previous increment unless the determinant of the tangent stiffness matrix has changed the sign, in which case a sign reversal is applied. Aa is the incremental arc length used for the next increment, and is determined by Au = C, (JD/J,)“2Au,
4. EQUILIBRIUM EQUATIONS AND CONVERGENCE CRITERION
(35)
and The nonlinear pressed by
equilibrium
equation
may be ex(36)
$=F-IP=O,
(31)
where $ is the unbalanced force between the internal nodal force vector F and the external nodal force 1P; 1 is a loading parameter and P is a normalized loading vector. The internal nodal force vector is obtained by summing up the element nodal force vectors in the global coordinate system. In this paper, a weighted Euclidean norm of the unbalanced force [31] is employed as the error measure of the equilibrium state during the equilibrium iterations, and the convergence criterion is given by (32) where N is the number of degrees of freedom for the disceretized structure and p,& is a prescribed value of error tolerance. Unless it is stated otherwise, the error tolerance is set to lo-’ in this paper.
where: Au, is the arc length used for the Ith increment; J, is the number of iterations required to achieve equilibrium for the Ith increment; JD is the desired number of iterations; ‘the safety factor’, C,, lies between 0.7 and 1.0, and ‘the cut parameters’ C, and C, are chosen to be 0.2 and 1.5, respectively, to prevent yielding of an incremental displacement which is too large or too small. Using the displacement increment obtained in eqn (33) and the method described in the previous section, the internal nodal force vector F in eqn (31) associated with the current configuration can be calculated. The loading parameter corresponding to the current configuration is given by I = A’+ AL, where I’ is the convergent loading parameter at the Ith increment and Al the loading parameter increment. Then the unbalanced force $ can be obtained from eqn (31). If the convergence criterion (eqn (32)) is not satisfied, a displacement correction r and loading parameter
777
Corotationai procedure for rotations of beam structures correction 61 [27,28] are added to the previous Aq and AL respectively to obtain a new incremental displacement and incremental loading parameter for the next iteration. The values of r and 612may be determined by
r=K;‘(-$
+61P)
(37)
and Au* = (Aq + r)‘(Aq + r),
(38)
where Kr may be the tangent stiffness matrix at some known configuration. This procedure is repeated until the convergence criterion is satisfied. It should be mentioned that, during the first few equilibrium iterations, the values of the element axial nodal forces obtained from the current deformation using eqns (6) and (29) may be several orders larger than their convergent values for certain problems. This may cause di~culty in convergence or even divergence for a large increment. In [IO] a two-cycle iteration scheme is introduced to overcome this difficulty. This scheme was proven to be very effective by numerous examples studied in [IO]; however, in the present study, it is found that the accuracy of convergent solutions obtained by this scheme is not sufficient for some problems. This difficulty probably arises due to the fact that the difference between the values of convergent element axial forces at the first and second cycles is not small for some problems. In order to overcome this difficulty an n-cycle iteration scheme, an extension of the two-cycle iteration scheme, is proposed in this study and can be described as follows. Let fdenote the element axial force corresponding to the current deformation, h denote the convergent element axial force of thejth cycle, and& denote the convergent element axial force of the fth increment. For the iterations of thejth cycle the element axial force F,, (eqn (29)) is replaced by
where f, and 4_ , are m x 1 column matrices containing the convergent values of element axial forces 4 andA_, , m is the number of the elements used for the discretization of the structure, and p, is a prescribed parameter. If eqn (41) is not satisfied, the iterations for next cycle are performed. Otherwise, the final cycle of iterations is carried out until eqn (32) is satisfied, At the final cycle, the element axial force f calculated from the current defo~ation is used in eqns (29) and (23) to obtain the element axial and bending forces. The convergent solution of the final cycle is used as the solution of the corresponding increment. 6. NUMERICAL STUDIES
Example 1. Cantilever beam with an end moment This example considered is a cantilever beam sub. jetted to a concentrated moment at the free end as shown in Fig. 5. The beam was d&ret&d by 10 elements. The results shown in Fig. 5 are obtained by
using only three increments. The number of iterations is about six per increment. As can be seen, the agreement with analytical solution is quite good. It should be mentioned that the two-cycle iteration scheme is used for this example because the element axial forces are zero for this problem. This example is extensively studied in the literature to demonstrate the efficiency of numerical methods and the large rotation capability of the beam, plate and shell elements. To the authors’ knowledge, only the present authors have achieved bending of the cantilever beam into a full circle by using only three increments. Example 2. Cantilever 4Megree bend with an end Iwd The bend as illustrated in Fig. 6 is curved in the horizontal pfane and subjected to a vertical load. The
+=====&=I”F W
F;,==(l-CJJ;_,+CJ
j21,
(39)
E = l.2X10skN/m2 V=O L = 10.0 m b = 1.0 m h = 0.1 m
where Ca f [0,11, is a prescribed parameter; the FL2 required for the evaluation of $ in eqn (18) and Fb in eqn (23) are replaced by J$=
“fit*
j=l
(I--C&,+C&,,
ja2,
I
(40)
where Cbe [0, I] is a prescribed parameter. The equilibrium iterations of the jth cycle are performed until the Euclidean norm of the unbalanced force vector in eqn (31) is smaller than a prescribed value, which may be chosen to be larger than the error tolerance given in eqn (32). Then the following inequality is checked:
I
I
-Analytical solution Present analysis l
12 6 8 10 Displacement (ml Fig. 5. Cantilever beam with an end moment. 0
2
4
Kuo-MO HSIAOet (11.
I2
R = 100. in u = 0. E = 1O’psi
t
c5jgft$;* 69.26
BEAM
l-l”-1 CROSS SECllON
+
I 61.44
E=4.32XlO’ (GJ), -
5 & .6
Ref *
t
[ 6 ]
Present
69.26
A,=.500
=4.15x105
(GJ),=1.66X10s
analysis
I
L=2OO.Oh A,=.100
(I,),
=.400
(Ix),
=.05
(ly),
=.133
(Iy)2=.05
PL2mJ,
2*o I 7---1.5
-
t/
Ref
“‘1 ;0.
1.
2. 3. 4. LOAD PARAMETER
5. k =
6. Pi’/El
7.
0.0
1.0
Example 3. Space arch frame
Figure 8 shows the structure load system and the load displacement curve. In addition to four vertical
(16.4
, 46.3
A (24.7
Fig.
7. Deformed
]
w,Lx,oJ
p”“““‘”
2.0
3.0
4.0
5.0
Fig. 8. Space arch frame.
Fig. 6. Forty-five-degree circular bend with an end force.
bend has an average radius of 100 in. and cross section area 1 it?. The bend is idealized using eight equal beam elements. At each increment, the stiffness matrix updating is only performed at the first five iterations of each cycle. Only four increments are used in this analysis. Three cycles of iteration are used per increment. The average number of iterations per increment is about 12. The results for the end displacements versus applied load are shown in Fig. 6 together with the solutions given by Bathe and Bolourchi [6] using eight beam elements and 60 equal load increments. Very good agreement between these two solutions is observed. The deformed configurations of the bend at various load levels are shown in Fig. 7.
[ 6
loads P, the structure is subjected to two lateral loads equal to 0.001 P. For all members the major principal axis of inertia x is normal to the plane of either arch rib. The symbols GJ, 1, and 1, denote, respectively, the torsional rigidity, the major and the minor principal moments of inertia, and the subscripts 1 and 2 denote the member groups as noted in the figure. Each member of the structure is idealized by four equal elements. The results of the present study shown in Fig. 8 are obtained by using two increments with the error tolerance p,,,, = lo-‘. Three cycles of iteration are used for both increments and the total number of iterations used for each increment are four and five, respectively. The present results are in excellent agreement with the solutions given in [8] which are obtained using 17 equal load increments (transcribed by the authors). Case A
:
Case
:
B
6 boundary nodes free in trondatlonal movsrnent 6 boundary nodes restrained against tronslotlonol movsment
, 53.6)
, 60.6
, 35.6)
shapes of a 45-degree circular bend.
-439600 lb/in2 =159000 lb/in2 10.494 iI+ PO.02 inz -0.02 it? =0.0331 in’
Fig. 9. Geometry of 12-member hexagonal frame.
Corotational
179
procedure for rotations of beam structures
& 1
V
I
x ,,w
2-
x,,u
E-3.03X1 0’ N/cm* G=l.O96Xl@ N/cd
J
5
DEfTLECllON
Fig. 10. Load-deflection
V (in)
curves for hexagonal frame.
Example 4. Twelve -member hexagonal frame
The hexagonal frame depicted in Fig. 9 is subjected to a concentrated load at the crown; two boundary conditions are considered: (a) six boundary nodes are free in translational movement and (b) six boundary nodes are restrained against translational movement. The frame was idealized using 12 (one element per member) and 36 (three equal elements per member) beam elements. For both boundary conditions, the results (not shown) using 12 elements are in close agreement with the solutions of Meek and Tan [9], who used a similar element but did not mention the number of elements used for discretixation. The results of case (a) using 36 elements are shown in Fig. 10, together with those reported by Meek and Tan [9] and Papadrakakis [7]. It is observed that the present results are in agreement with that of Papadrakakis, which is in exact agreement with experimental results given by Griggs [33]. The present results are obtained using five increments; the number of cycles used is about four per increment, and the
250-
Case B o Preset7 t -
Ref [S]
200-
Fig. 11. Load-deflection
curves for hexagonal frame.
Fig. 12. Geometry of 24-member shallow dome.
total number of iterations used per increment is about 11. The stiffness matrix is updated only at the first two iterations of each cycle. The results of case (b) using 36 elements are shown in Fig. 11. The present results are obtained using six increments; the number of cycles used per increment is about four, and the total number of iterations used per increment is about 13. The dashed curve shown in Fig. 11 is also obtained by the present study using 16 increments. The stiffness matrix is updated only at the first two iterations of each cycle. Also shown in Fig. 11 are the solutions reported by Meek and Tan [9]. The discrepancies between these two solutions may be explained by suggesting that the number of elements used in [9] is insufficient. Example 5. Twenty-four-member shaped shallow dome
hexagonal
star-
Figure 12 shows the geometry of a 24member hexagonal star-shaped shallow dome. The supports of the dome are assumed to be pinned and restrained against translational motion. The dome is idealized using 24 (one element per member) and 72 (three equal elements per member) beam elements. For all loading conditions, the results (not shown) using 12 elements are in close agreement with the solutions of Meek and Tan[9] who used a similar element but did not mention the number of elements used for discretixation. The first loading condition considered is that of a concentrated vertical load at the apex of the dome. The present results using 72 elements, shown in Fig. 13, are obtained by using three increments. The number of iterations for each cycle is also given in parenthesis beside each point on the graph. The average number of iterations required for one increment is about eight. The stiffness matrix is only updated at the Crst two iterations at each cycle. By keeping the same member cross sectional area but decreasing the tlexural stiffness in the vertical plane, the structure is reanalyzed. The graph of the load-
780
Kuo-MO HSMO er al. 2.0
2.5? .
Present
-
2.0
(I 2.1, I ) 0
I 0
Present
-Ref
Rsf [ 9 ]
[9]
l
A
P*V
OK
. Loaded node
Fig.
1.0
2.0 DEFIJXTION
13. Load~efl~tion
3.0 V (cm)
co
curves for con~ntrat~ Ioad.
central
2.0 0
Prsssnt
-Rsf
[S]
1.5 *
a0.
12=2.377cm' 13=0.295cm4
9 s
/ l
Loaded node
/
“‘l/---j Of
Fig.
’
J PO.91 Bcm’
1.0,
1.0 2.0 DEFLECTION
14. ~ad~efl~tion
curves load.
3.0
‘4.0
5.0
V (cm)
for ~n~trat~
central
2.0
0 -
curves
12=2.377cm4 13=0.295cm4 J =0.918cm4
3.0 4.0 V (cm)
3.0
for unsymmetrical loading.
5.0
deflection curves using 72 elements is shown in Fig. 14. As can be seen, only three increments are used. The average number of iterations per increment is nine. The second load condition is al1 nodes loaded symmetrically. The results using 72 elements are shown in Fig. 15. Four increments are used and the
F
1.0 2.0 DEFLECTION
Fig. 16. Load-deflection
0
!wc&ed
Present Ref [ 9 ]
0
average number of iterations used per is eight. For the unsymmetrical loading shown in Fig. 16, five increments are the average number of iterations used increment.
increment condition used and is 10 per
7. CONCLUSIONS
A practical motion process of the three dimensional beam element is presented to remove the restriction of small rotations between two successive increments for large displacement and large rotation analysis of space frames using incremental-iterative methods. The nonlinear fo~uiation is based on the corotational formulation by which the major geometric nonlinearities were shown to be embodied in the coordinate transformation when forming the element assemblage. The transformation of the element coordinate system is assumed to be accomplished by a ~anslation and two successive rigid body rotations: a transverse rotation followed by an axial rotation. The element formulation is derived based on the small deflection beam theory with the inclusion of the effect of axial force in the element coordinate system. The element internal nodal forces are calculated using the total defo~ational nodal rotations. Two methods, referred to as direct method and incremental method, are proposed in this paper to calculate the total deformational rotations. Despite the fact that the formulation of the beam element is very simple, highly accurate solutions are obtained. It is believed that the use of a simple element combined with the corotational formulation and the process of element motion proposed in this paper may represent a valuable engineering tool for the solution of nonlinear spatial beam problems. REFERENCES
0’
H. Chu and R. H. Rampetsreiter, Large deflection buckling of space frames. J. Srrucr. Do., AXE St%, 2701-2722 (1972). 2. C. Oran, Tangent stillness in space frames. J. Srrucf. Diu., ASCE 99, 987-l 101 (1973). t. K.
1.0
2.0 DENCTION
Fig. 15. Load-deflection
3.0 V (cm)
4.0
5.0
6.0
curves for symmetrical loading.
Corotational
procedure for rotations of beam structures
3. T. Belytschko and B. J. Hsieh, Non-linear transient finite element analysis with convected co-ordinates. Inl. .I. Nwner. Meth. Engng 7, 255-271 (1973). 4. T. Belytschko, L. Schwer and M. J. Klein, Large displacement transient analysis of space frames. Inl. J. Numer. Meth. Engng 11, 65-84 (1977). 5. S. N. Remseth, Nonlinear static and dynamic analysis of framed structures. Camp. Struct. 10,879-897 (1979). 6. K. J. Bathe and S. Bolourchi, Large displacement analysis of th~e~imen~onal beam structures. Inr. J. Numer. Meth. Engng 14, 96I-986 (1979). 7. M. Papadrakakis, Post-buckling analysis of spatial structures by vector iteration methods. Compur. Struct. 14, 393-402 (1981). 8. R. K. Wen and J. Rahimzadeh, Nonlinear elastic frame analysis by finite element. J. Srrucf. &gng, ASCE 109, 1952-1971 (1983). 9. J. L. Meek and H. S. Tan, Geometrically nonlinear analysis of space frames by an incremental iterative technique. Compur. Meth. appl. Mech. Engng 47, 261-282 (1984). 10. K. M. Hsiao and F. Y. Hou, Nonlinear finite element analysis of elastic frames. Compur. Srruct. 26, 693-701 (1987). 11. J. H. Argyris and P. C. Dunne, On the application of the natural mode technique to small strain large displacement problems. World Congress on Finite Element Methods in Structural Mechanics, Boumemouth, October 1975. 12. G. Wempner, Finite elements, finite rotations and small strains of flexible shells. Inr. J. Solids Srrucf. 5, 117-153 (1969). 13. J. H. Argyris, P. C. Dunne, G. A. Malejannakis and E. Schelkle, A simple triangular facet shell element with application to linear and nonlinear equilibrium and eiastic stability problems. Comput. Meih. appl. Mech. Engng IO, 371-403; 11, 97-131 (1977). 14. E. Ramm, A plate/shell element for large deflections and rotations. In Forn~ulotionsand ComputationalAfgorithmsin Ftiite Element Analysis,U. S.-Germany Symp. (Edited by K. J. Bathe, J. T.Oden and W. Wun~eriic~), DO. 264-293. MIT Press. Cambridae. MA (1977). IS. i;‘. Horrigmoe and P. G.. Bergan, Nonlinear-analysis of free-form shells by flat finite elements. Comput. Meth. appl. Mech. Engng 16, I l-35 (1978). 16. T. Belytschko and L. Glaum, Application of higher order corotational stretch theories to nonlinear finite element analysis. Comput. Struct. 10, 175-182 (1979). 17. K. J. Bathe and S. Bolourchi, A geometric and material nonlinear plate and shell element. Compur. Srrucr. 11, 23-48 (1980). 18. K. J. Bathe and L. W. Ho, A simple and effective element for analysis of general shell structures. Comput. Srrucr 13, 673-681 (1981).
781
19. T. J. R. Hughes and W. K. Liu, Nonlinear finite element analysis of shells: Part I. Three-dimensional shells; Part II. Two-dimensional shells. Comput. Meth. appl. Mech. Engng 26, 331-362 (1981); 27, 167-181 (1981). 20. J. H. Argyris et of., Finite element method-the natural approach. Comput. Meth. appi. Mech. Engng 17/l& l-106 (1979). 21. J. H. A&is, An excursion into large rotations. Cornmu. Meth. aDDI.Mech. Enpna 32. 85-155 (1982). 22. K. S: Surana, GiometticaIIy n&l&ear formulation for the curved shell elements. Inr. J. Nwner. Meth. Engng 19, 581-615 (1983). 23. C. C. Rankin and F. A. Brogan, An element independent corotational procedure for the treatment of large rotations. 3. Press. Vessel Technol., ASME 108,165-174 (1986). 24. K. M. Hsiao, Nonlinear analysis of general shell structures by flat triangular shell element. Camp. Struct. 25, 665-675 (I 987). 25. J. H. Argyris, S. Kelsey and H. Kamel, Matrix methods of structural analysis: a p&is of recent developments, In Matrix Methods of Structural Analysis. AGARDogruph 72 (Edited by B. F. de Veubeke), pp. l-165. Pergamon Press, London (1964). 26. K. Mattiasson and A. Samuelsson. Total and undated Lagrangian forms of the co-rotational finite eiement fo~ulation in g~rnet~~Ily and materially nonlinear analysis. In Nu&rical Meth&for Non-line& Problems, Vol. 2 (Edited bv C. Tavlor. E. Hinton and D. R. J. Gwen), ‘pp. 13Ll51. PiGeridge Press, Swansea, U.K. (1984). 27, M. A. Crisfield, A fast incremental/iterative solution procedure that handles “snap-through”. Compul. Sfrucl. 13, 55452 (1981). 28. T. Y. Chang and K. Sawamiphakdi, Large deflection and post-buckling analysis of shell structure. Comput. Meth. appl. Mech. Engng 32, 31 l-326 (1982). 29. D. J. Dawe, Matrix and Finite Element Displacement Analysisof Structures. Oxford University Press, London (1984). 30. H. Goldstein, Classical Mechanjcs, 2nd Edn. Addison-W~lev. Readinn. MA 11980). 31. A. K. Noor &d J. My Peters, Tracing post-hmit point with reduced basis technique. Comput.Meth. appi. Mech. Engng Zs, 217-240 (1981). 32. P. G. Etergan and T. Ssreide, Solution of large displacement and instability problems using the current stiffness parameter. In Finite Elements in Non-linear Mechanics (Edited by P. G. Bergan er al.), pp. 647-669. Tapir Press. Trondheim. Norwav (19781. 33. H. fi. Griggs, Eiperimekal stidy of instability in elements of shallow space frames. Research Report, Dept. of Civil Engineering, MIT, Cambridge, MA (1966).
6 S~~~rures Vol. 27. No. 6. pp. 769-781.
1987 0
Printed in Great Britain.
aa-7949187 13.00 + oxa 1987 Pcrgamon Journals Lid.
A COROTATIONAL PROCEDURE THAT HANDLES LARGE ROTATIONS OF SPATIAL BEAM STRUCTURES Kuo-MO HSIAO,HORNG-JANNHORNG and YEH-REN CHEN ~~rtment
of Mechanic& Engineering, National Chiao Tung University, Hsinchu. Taiwan, Republic of China (Received 5 January 1987)
Abstract-A practical motion process of the three dimensional beam element is presented to remove the restriction of small rotations between two successive increments for large displacement and large rotation analysis of space frames using incremental-iterative methods. In order to improve convergence properties of the equilibrium iteration, an n-cycle iteration scheme is introduced. The nonlinear formulation is based on the corotational formulation. The transformation of the element coordinate system is assumed to be accomplished by a translation and two sueesssive rigid body rotations: a transverse rotation followed by an axial rotation. The element formulation is derived based on the small deflection beam theory with the inclusion of the effect of axial force in the element coordinate system. The membrane strain along the deformed beam axis obtained from the elongation of the arc length of the beam element is assumed to be constant. The element internal nodal forces are cakulated using the total defo~atioM1 nodal rotations. Two methods, referred to as direct method and incremental method, are proposed in this paper to calculate the total defo~~ona~ rotations. An incremental-iterative method based on the Newto~Rap~n method combined with arc length control is adopted. Numerical studies are presented to demonstrate the accuracy and efficiency of the present method.
1. INTRODUCTION
The development of new and efficient formulations for the nonlinear analysis of beam structures has attracted the study of many researchers in recent years, and different alternative formulation strategies and procedures to accomodate large rotation capability during the large defo~ation process have been presented [I-lo]. These formulations can be divided into three categories: Total Lagrangian (TL) formulation, Updated Lagrangian (UL) formulation and Corotational (CR) formulation. It should be noted that within the corotating system either a TL or a UL formulation, or even a formulation based on a small deflection theory may be employed. The large number of publications on the nonlinear analysis of beam structures is, at least partially, due to the fact that various kinematic nonlinear formulations can be employed. It seems that large rotations in plane frames pose no major problem. Hsiao and Hou [IO] introduced a simple and effective corotational formulation of beam element and numerical procedure, which can remove the restriction of small rotations between two successive increments for the large displacement analysis of plane frames using incremental-iterative methods. Unfortunately, the method presented in [lo] cannot be applied to three dimensional frames. The difficulty of obtaining effective solutions is particularly pronounced in the analysis of spatial beam structures; a general three dimensional nonlinear formulation is not a simple extension of a two dimensional fo~ulation, because large rotations in three dimensional analysis are not true vector
quantities; that is, they do not comply with the rules of vector operations and the result will in general depend on the order in whch the rotations are taken. This point has been throughly discussed by Argyris [1 I] and Wempner [ 121. The problem of large rotations on space structures has received wide attention in the ~terature; many different strategies based on the TL, the UL, or the CR formulations have been reported, those of [Z-S, 6, 1l-241 being only a small fraction of the total. In 1191Hughes and Liu developed a specialized shell element which can handle arbitrarily large rotations. Argyris has covered the subject of corotational coordinates extensively including a lengthy discourse on the subject of large rotations [20,21]. Belytschko er al. [3,4, 14 have applied corotational formulation to the dynamic analysis of space frames where arbitrarily large rotations can be expected. Horrigmoe and Bergan (IS] have su~~fully applied a corotational approach to their shell elements Rankin and Brogan [23] have introduced a corotational procedure which may enable existing shell element formulations to be used in problems that contain arbitrarily large rotations. Recently Hsiao 1241 has proposed a motion process for triangular shell elements to remove the restriction of small rotations between two successive increments for nonlinear shell analysis using incremental-iterative methods. The wide range of numerical examples studied in [3 4,7,9, 11 , 13, 15, 16,18,20,‘21,23-261 indicates that the corotation approach, first described by Argyris et al. [25], may be very useful in the analysis of spatial structures containing arbitrarily large
770
Kuo-MO HSIAOet al.
rotations. However, most strategies based on the corotational formulation suffer from one inherent drawback: they are restricted to small rotations between two successive load increments during the deformation process. This limitation arises because the incremental nodal rotations are considered to be vector quantities. Although the method introduced in [24] may remove this restriction for triangular shell elements, unfortunately, this method cannot be applied to the space beam elements, because, unlike the shell elements, the element coordinate of the space beam elements cannot be determined using only nodal coordinates. The objective in this paper is to present a practical motion process of the three dimensional beam element which can remove the restriction of small rotations between two successive increments for large displacement and large rotation analysis of space frames using incremental-iterative methods. Here the motion process presented in [24] is modified to accomodate the characteristics of the motion of beam elements. In this paper the transformation between the element coordinate systems is assumed to be accomplished by a translation and two successive finite rotations: a lateral rotation about an axis perpendicular to the current beam axis followed by an axial rotation about the current beam axis. Two methods, termed direct method and incremental method, are introduced to describe the motion process of the beam element and to determine the total deformational nodal rotations. The dominant factors in the geometrical nonlinearities of space structures are attributable to finite rotations, the strains remaining small. For space frames discretized by finite elements, this implies that the motion of the individual elements with proper size will, to a large extent, consist of rigid body motion.
If the rigid body motion part is eliminated from the total displacements, the deformational part of the motion is always a small quantity relative to the local element axes; thus, incorporated with the corotational formulation, the small deflection beam theory with the inclusion of the effect of axial force is adopted here to deal with the large rotations but small strains problems. The numerical algorithm used here is an incremental-iterative method based on the NewtonRaphson method combined with constant arc length of the incremental displacement vector [27,28]. In order to improve the convergence properties of the equilibrium iterations, an n-cycle iteration scheme, which is an extension of the two-cycle iteration scheme proposed in [lo], is introduced. Numerical examples are presented and compared with the results reported in the literature to demonstrate the accuracy and efficiency of the proposed method. 2. COORDINATESYSTEMS One of the basic considerations in formulating nonlinear structural problems is the selection of a description of motion. In the present study, the corotational (CR) formulation is adopted. In this formulation, each element is associated with a local Cartesian element coordinate system 2, (i = 1,2, 3) as shown in Fig. 1, that rotates and translates with the element but does not deform with the element. The element coordinate system used herein is a right handed one and is defined as follows. The origin of the element coordinate system is located at the local node 1, the 2, axis passes through the centroids of member end sections, and the Z2 and 5, axes are parallel to the principal directions of the undeformed end cross section.
Fig. 1. Coordinate systems, member deformations and associated forces.
Corotational
771
procedure for rotations of beam structures
Also shown in Fig. 1 is a fixed global coordinate system xi (i = 1,2,3) used to define the location of the nodal points. We note that the equilibrium equations of a structure are written in the fixed global coordinate system, and the incremental nodal parameters of the system of equations are also calculated in this coordinate system. The element equations are first formulated in the element coordinate system, and then transformed to the global coordinate system using standard procedure [29] prior to the element assemblage process. If we consider a vector B (or B if measured in the _Zicoordinate system) with global components Bi (i = 1,2,3), and element coordinate components Si (i = 1,2,3), we have the following transformation:
P !!
9
R’
n
Fig. 2. Finite rotation of vector.
respectively; 4 is the angle of counterclockwise rotation, and n is the unit vector along the axis of rotation. 3. CR-FORMULATION OF BEAM ELEMENT
where clti= cos(xi, Xj). The beam element employed here has two nodes with six degrees of freedom per node (Fig. 1): these are the translations 0, in the fi (i = 1,2,3) directions at nodes j (j = 1,2), and the rotations flj about the Zi axes at nodes j (j = 1,2). The global nodal parameters for the system of equations associated with the individual elements are chosen to be the translations Vii in the xi (i = 1,2,3) directions at nodes j (j = 1,2), and rotations 0, about the xi (i = 1,2,3) axes at nodes j (j = 1,2). In this study an incremental-iterative method is used to solve the nonlinear equilibrium equations. Both the incremental nodal translations and the incremental nodal rotations are regarded as vector quantities. Thus, using eqn (l), the nodal vectors, AUj= {AU,j, AU,, AUjj) and AtIj = {A8ij, A&, Ae,}, referred to the global coordinate system can be transformed to AUj = {AO,j, AUq, AOJj} and A4 = {A&, A&, Ae,}, referred to the element coordinate system, respectively. It should be noted that Aej cannot be interpreted as component rotations about Cartesian axes &. In this study, Ae, = {Ae,, 0,0}, the components of Ah along the 2, axis, and Afj”,= (0, A&, AB,), the components of A4 perpendicular to the R, axis, are considered to be rotation vectors to define rotations, details of which will be discussed later. For convenience of the later discussion, the term ‘rotation vector’ is used to represent a finite rotation. Figure 2 shows a vector R which as a result of the application of a rotation vector 4n is transported to a new position R’. The relation between R and R’ may be expressed as [30] R=cos+R+(l
The formulation of the beam element developed here is applicable to arbitrarily large rotations but restricted to small rotations relative to the element axis. The beam element is formulated in the element coordinate system based on the small deflection beam theory with the inclusion of the effect of axial force. The element, as shown in Fig. 1, has two nodes with six degrees of freedom per node, and can transmit an axial force, two shear forces, two bending moments and a torque. Herein the beam element is assumed to be straight and of constant cross section. The cross section is doubly symmetric, thus excluding coupling of the torsional stiffness to that of bending and axial stiffness. Shearing deformations and warping effects are neglected. The material is assumed to be linearly elastic. The element stiffness matrix is obtained by superimposing its bending, geometric, torsional and axial stiffness matrices. The element internal nodal forces are evaluated using the total deformations. 3.1. Kinematics of beam element It is assumed that the lateral deflection curves of the beam member are the cubic Hermitian polynomials in the & and .?s directions of the element coordinates (the principal directions of the undeformed end cross section), and that the axial rotation varies linearly along the member. The membrane strain along the deformed element axis is assumed to be constant. Thus, the membrane strain can be evaluated from the elongation of the arc length of the member. The lateral deflection curves of the beam member may be given by
-cos4)(n.R)n + sin 4(n x R),
(2)
where * and x denote the dot and the cross product
where a bar over a quantity denotes that it is defined in the element coordinate system. P and Ware lateral
IWO-MOHstao et al stant membrane strain along the deformed beam axis, the membrane strain of the beam axis can be written as Cl = (S - SJ)/&,
“2
/
undeformed
end section
Fig. 3. Deformational nodal rotations.
deflections in the & and i, directions, respectively. iii, and ej (i = 2,3; j = f , 2) are the nodal displacements and rotations shown in Fig. 1. N, (i = i,4) are shape functions and are given by N, = l/4(1 - #(2
+ t)
N2 = c/8(1 - r’)(l - r)
(6)
where S, = L is the initial arc length of the beam axis. Figure 3 shows that the normals of the undeformed element end sections at nodes j (j = 1,2), 6, are rotated to g,,,, the current deformed normals of the element end sections by the rotation vectors e)l,, which is ~~ndicular to the f,axis of the element coordinate system. The representations of the lateral deformational nodal rotations are based on the assumptions that these rotations are small. On the basis of this assumption e,,, the 5, (i = 2,3) components of the rotation vectors 6,,j are chosen to be the Iateral defo~ational nodal rotations about Zi axes at nodes j. Note that the direction of the undeformed normal of the element end sections coincides with the positive direction of the 5, axis. Thus, the second subscript j of gq is omitted throughout this paper. 3.2. Determination of element coordinate system and element de~ormat~o~al rotations
Assume that the incremental-iterative method is used for the solution of nonlinear equilibrium equaNl = c/8(1 + c*)(l + r), (4) tions and the equilibrium configuration of the Ith increment is known. Let AUj and A@,(j = 1,2) be where c = Zi, - %ir is the current chord length of the incremental nodal displacement and rotation the beam member, and f, are the 2, coordinates of vectors of an element at nodes j extracted from the incremental nodal parameters of the system of equanodes j (j = I, 2) in the element coordinate system; tions. At this point, an interesting and relevant 5 = - I + 23,/c is a nondimensional coordinate. question arises. Given the incremental nodal disNote that in this study, the relative displacements placements and rotations, how are the current of the elements are referred to the element coodinate system, Due to the definition of the element co- element coordinate system, the axial and lateral ordinate system, the lateral nodal displa~ments 041 defo~ationa1 nodal rotations dete~ined? Let ‘xj (j = 1,2) denote the node coordinate (i = 2,3) at the nodal points j (j = 1,2) are identical vectors of an element in its Ith equilibrium conto zero. The nonzero deformational nodal displacefiguration; the current node coordinate vectors xj ments of an element can be divided into the axial are obtained by adding the incremental nodal disrelative displacement, the axial relative rotation and lateral deformational rotations. The element defor- placement vectors AUj, so that mation can be decomposed into the membrane deformation, the torsional deformation and the flexural xj = ‘x, + AU,. (7) deformation. Herein the flexural deformation is deThe 3, axis of the current element coordinate termined by the lateral deformational nodal rotations system can then be constructed using x, given in using elementary beam theory and the torsional eqn (7) and the definition of the element coordinate deformation is determined from the axial relative nodal rotation, while the membrane deformation is system. But, unlike the cases of triangular shell obtained from the change of arc length of the beam elements fl5, 18,241, the .?r and 5 axes cannot be determined using only the node coordinates. The axis which can be calculated from the lateral determination of the current element coordinate deflection. If the arc length of the beam axis is expressed by system will be discussed in the process of element motion. 1 For determining the element deformational nodal (1 + p’2 -i- @‘2)“2 d&, S=c/2 (5) rotation and element rigid body rotations, we proi -1 pose two methods to describe the process of the element motion in this paper. In the first method, where c is the current chord length of the beam member, ( )’ denotes x,-derivatives and P and rii are referred to as the direct method, the lateral deformational nodal rotations are determined from the given in eqn (3), then from the assumption of conN,= l/4(1 +S)‘(Z--r)
Corotational
773
procedure for rotations of beam structures
(b)
(d)
Fig. 4. Process of element motion.
orientations of the undeformed and deformed normals of the element end sections. In the second method, referred to as incremental method, the total lateral deformational nodal rotations are calculated by incrementation. For both methods, the total twist nodal rotations are determined by incrementation. The processes of element motions and the methods corresponding to these motion processes to determine the deformational nodal rotations and element coordinate system are described as follows. (a) Direct mefhod. The process of element motion is divided into the following six steps.
1. A rigid body translation by AU,. The whole element is translated by AU,, where AU, is the incremental nodal displacement vector of node 1. The origin of the ‘3, axes is translated to the origin of the Zi axes as shown in Fig. 4(a). 2. A lateral rigid body rotation by the rotation vector 8’. The rotation vector 6’ (referred to *ii coordinate system) is given by &‘=cos-‘(‘e,*e,)
42,x e, II‘el x e,II’
(8)
where ‘e, and e, are unit vectors associated with the
774
Kuo-MO
HSIAOef al.
‘.f, and x’, axes rcspcclivcly. The rotation vector oi’, passing through node 1, is applied to the whole element except ‘tQ, the deformed normals of the element end sections at nodes j (j = 1,2) at the Ith equilibrium configuration. Here it is assumed that ‘e,,, are not rotated by the rotation vector ci’, but are translated with the motion of nodes j (Figs 4(b) and (d)). The ‘& axes are rotated by an angle about the axis perpendicular to the ‘2, and f, axes (Fig. 4(a)), and the resultant coordinate system is labeled 2; axes (Fig. 4(c)). As can be seen, the 2; axis coincides with the 2, axis. 3. Finite rotations of ‘C,+by the rotation vectors Afig. The deformed normals ‘c?~are rotated to &, by the application of the rotation vectors A6:, as shown in Figs 4(c) and (d), where AtiLj are the components of Afij (j = 1,2) (the given incremental nodal rotation vectors referred to 2: coordinate system) perpendicular to the Xi axis. 4. Twist rotations by A4j. The rotation vectors A4j are applied to nodes j (j = 1,2), where AtJj are given by A& = A&, - /I
(9)
/? = l/2(6&, + A&,),
(10)
in which tI,, are the components of A& along the 5; axis as shown in Fig. 4(c). 5. A stretch by (c - ‘c)e,. Node 2 is translated along the 2, axis (Fig. 4(d)), where c and ‘c are chord lengths of the element corresponding to the current configuration and the equilibrium configuration of the Ith increment. 6. An axial rigid body rotation by g. The rotation vector (eqn (10)) passing through node 1 is applied to the whole element. The intermediate axes, ai, are rotated about the 2; axis by an angle 11 fi 1)to produce the & axes of the current element coordinate system as shown in Fig. 4(f). Vectors El and 5; shown in Fig. 4(d) are rotated by this rotation vector to reach their final positions Z, and &, (not shown in Fig. 4(f)). ‘I’he orientations of the current element coordinate axes fi may be obtained from the rigid body rotations caused by the rotation vectors given in steps 2 and 6 of the above process. The orientation of the deformed normals gdi and the undeformed normal 5, of the element are also determined by the above process of motion. Thus, the lateral deformational nodal rotation vectors, e,, referred to the current element coordinate system may be expressed as
ii, 11
0 e;I,= IT4 cos-‘(E:~~)
x *dj
II& x ~,$II’
0 e;.=
f$, = ‘4 + (A8, - a).
(13)
1103, For simplicity of computation, only eqn (13) is used to calculate the total lateral deformational nodal rotations for the numerical examples studied in this paper. It is believed that identical results will be obtained when eqn (11) is used.
(11)
O3j
The current axial deformational may be obtained by
whcrc ‘4, is the axial dcformational nodal rotation vector of the Zth equilibrium configuration at nodes j (j = 1,2); A$j is given in eqn (9). (b) Incremental method. In this method, the motion process is also divided into six steps. Only steps 2 and 3 of this process are different from those in the direct method, and are described below. 2. A rigid body rotation by o?‘.The rotation vector 6i’ (eqn (8)) passing through node 1 is applied to the whole element. As can be seen in Figs 4(b) and (e), the deformed normals of the Ith increment ‘gdl are rotated to an intermediate position I&;. 3. Deformational rotations by (68;, - a’). The intermediate deformed normals *EL are rotated to other intermediate positions 5; as shown in Fig. 4(e) by the application of the rotation vectors (A&, - o?‘), in which A#Aj are the components of AJj (j = 1,2) perpendicular to 2; axis, and OS’is the rotation vector given in eqn (8) (referred to 2: coordinate system). Vectors AgAj, oi’ and ‘6; (Fig. 4(e)) are rotated by the rotation vector [ (eqn (10)) to their final positions Afimj,6 and ‘4 (referred to the current element coordinate system &), which are not shown in Fig. 4(f). The element coordinate obtained from this process is the same as that constructed by the direct method. The orientations of the deformed and undeformed normals of the element end sections can be determined from this process of motion as well. Thus, the deformational nodal rotations can be calculated using eqn (11). Due to the assumption of small deformational nodal rotations, the vector operations might be valid for the deformational nodal rotation vectors. Thus, an alternative, referred to as incremental method, is introduced here. The concept of this method is similar to that in [l&24]. The total deformational rotations of the current configuration referred to the current element coordinate are obtained by adding the incremental nodal rotations (A&j-oS) to the deformational nodal rotations of the equilibrium configuration of the Ith increment, and are expressed as
3.3. Element st1~ne.s.smatrix
nodal rotations
The total element stiffness matrix is formulated by superimposing the bending stiffness matrix i& and geometric stiffness matrix $ of the basic beam element, and the axial stiffness matrix g,,, and the (12) torsional stiffness matrix g, of the linear bar element. The derivation of these matrices is well documented
115
Corotational procedure for rotations of beam structures in the textbooks and thus will not be repeated here. However, these matrices are given as follows. (a) Bending stiffness matrix K,:
where GJ is the torsional rigidity and L is the initial arc length of the beam axis. 3.4. Element nodal force vectors
(14) where r
12 -6L -12
-6L 4L= 6L
-12 6L 12
2~~
6L
12 6L -12
6L 4L’ -6L
-12 -6L 12
6L
2L2
L-6L
-6L 1 2L’ 6L
(1%
4~2J
and f
6L 1 2L2 -6L 9 (16)
The nodal force vectors of the elements corresponding to the global coordinate system are evaluated first in the current element coordinate system, and then transformed to the global coordinate system using standard procedure. Since small deformations are assumed, the element nodal forces can-in the element coordinate system-be evaluated in much the same way as in linear analysis. For linearly elastic material properties the element nodal force vectors can be calculated as follows. (a) Bending nodal force vector t,: &i=@,+K$,,
where L is the initial length of the beam axis, and EI, and E13are the flexural rigidities about the 4 axis and 2, axis respectively. The degrees of freedom corresponding to $ are ob=
{~,d2,,
&2.42*
~2,r&,,
(17)
022,42}r
where uU are nodal translations and gU are nodal rotations as shown in Fig. 1. (b) Geometric stiffness matrix &:
rt,= [
$2 0
0 rt,,
1
(18)
36 -3L
-3L 4L2
-36
3L
-36
-3L -L2
3L 36
3L
3L
4L2
4L2 -3L3L
-3L -3636
-3L3L
-L2
-3L
-L2
-3L and
K*&
3L [ -3636 3L
1 (19)
1
-1
[ -1
9 (20)
4L2 I
1
1’
(21)
where AE is the axial rigidity and L is the initial arc length of the beam axis. (d) Torsional stiffness matrix K,: z(
GJ
‘=L
[ 11’ 1
-1
-1
where ~b={F31rn221,F,2,11-3,,F2,,~,,,F22r1U32} is shown in Fig. 2; Kb is the bending stiffness matrix given in eqn (14); ab is the total bending deformation vector given in eqn (17). Note that due to the definition of the element coordinate system, the only nonzero elements in ab are f$, the deformational nodal rotations at nodes j (j = 1,2) about the fi (i = 2,3) axes, which may be obtained by using eqns (11) or (13). (b) Axial nodal force vector Fm;,: The element internal nodal forces are calculated by the total nodal deformation rotations. The axial nodal force vector F,,, = {I?,], F,2} (see Fig. 1) can be evaluated by introducing nodal virtual displacements 80, = 6 {D,, , u12} at nodes 1 and 2 in the f, direction, and equating the work done by the axial nodal force F,,, going through the virtual displacement JO,,, to the work done by the internal stress resultant T going through the virtual strain St,,, (that corresponds to the imposed virtual displacement) along the deformed beam axis as s
-L2
where L is the initial arc length of the beam axis and F,, is the axial nodal force at node 2. The degrees of freedom corresponding to KE are the same as that corresponding to izb. (c) Axial stiffness matrix KM: it_=:
(23)
4L2
-6L
(22)
cm,Fm =
I
Tsr, dS,
(24)
0
where S is the current arc length of the beam axis. The stress resultant T may be given by T = AEF-,,
(25)
where A is the cross section area, E is Young’s modulus and r,,, is the membrane strain of the current deformation. From eqn (6), the virtual strain a<,,, corresponding to ~50, can be expressed as
arm= as/s, (1 + p’2+ p2)‘nd{
So
(26)
Kuo-MO HSIAOelal.
776
in which 6c is the variation of the chord length of the beam axis with respect to SD,,,. Substituting eqns (S), (25), and (26) into eqn (24) gives
Since the virtual nodal displacements aa, are arbitrary, the axial nodal forces are obtained from eqn (27) as
Because the assumption of the small strain and small deformation, S2/S,c in eqn (28) is approximated by unity in this paper, and eqn (28) is thus reduced to
5. SOLUTION
ALGORITHM
An incremental-iterative method based on the Newton-Raphson method is adopted here. In order to deal with the limit points and snap through, the arc length of the incremental displacement vector is kept constant during the equilibrium iteration using Crisfield’s method [27,28]. An n-cycle iteration scheme is introduced here to improve the convergence characteristics of the equilibrium iteration. If the equiiibrium configuration of the Ith increment is assumed to be known, the system tangent stiffness matrix KT then can be calculated at this configuration and an initial displacement increment Aq for the next increment may be obtained by using Euler predictor as Aq = A%,
(33)
where Al is the initial incremental loading parameter and qT= Kf'P is the tangential displacement of unit loading P.For all increments other than the first, Ai. is obtained in much the same way as that mentioned in [27] and is given by
for numerical computation. (c) Torsional nodal force vector F,:
A1 = fAa(q;q,)“2, F,={;;;}=T{_;},
(30)
where 4 = ($,, - $,,) is the total relative rotation of the member about the 2, axis, r&j are the total twist rotations at nodes j (j = 1,2) about the fi axis and are obtained in eqn (12), and &?,j (j = 1,2) are twist moments shown in Fig. 1.
(34)
where the sign is chosen following an approach due to Bergan and Ssreide [32] in which the sign follows that of the previous increment unless the determinant of the tangent stiffness matrix has changed the sign, in which case a sign reversal is applied. Aa is the incremental arc length used for the next increment, and is determined by Au = C, (JD/J,)“2Au,
4. EQUILIBRIUM EQUATIONS AND CONVERGENCE CRITERION
(35)
and The nonlinear pressed by
equilibrium
equation
may be ex(36)
$=F-IP=O,
(31)
where $ is the unbalanced force between the internal nodal force vector F and the external nodal force 1P; 1 is a loading parameter and P is a normalized loading vector. The internal nodal force vector is obtained by summing up the element nodal force vectors in the global coordinate system. In this paper, a weighted Euclidean norm of the unbalanced force [31] is employed as the error measure of the equilibrium state during the equilibrium iterations, and the convergence criterion is given by (32) where N is the number of degrees of freedom for the disceretized structure and p,& is a prescribed value of error tolerance. Unless it is stated otherwise, the error tolerance is set to lo-’ in this paper.
where: Au, is the arc length used for the Ith increment; J, is the number of iterations required to achieve equilibrium for the Ith increment; JD is the desired number of iterations; ‘the safety factor’, C,, lies between 0.7 and 1.0, and ‘the cut parameters’ C, and C, are chosen to be 0.2 and 1.5, respectively, to prevent yielding of an incremental displacement which is too large or too small. Using the displacement increment obtained in eqn (33) and the method described in the previous section, the internal nodal force vector F in eqn (31) associated with the current configuration can be calculated. The loading parameter corresponding to the current configuration is given by I = A’+ AL, where I’ is the convergent loading parameter at the Ith increment and Al the loading parameter increment. Then the unbalanced force $ can be obtained from eqn (31). If the convergence criterion (eqn (32)) is not satisfied, a displacement correction r and loading parameter
777
Corotationai procedure for rotations of beam structures correction 61 [27,28] are added to the previous Aq and AL respectively to obtain a new incremental displacement and incremental loading parameter for the next iteration. The values of r and 612may be determined by
r=K;‘(-$
+61P)
(37)
and Au* = (Aq + r)‘(Aq + r),
(38)
where Kr may be the tangent stiffness matrix at some known configuration. This procedure is repeated until the convergence criterion is satisfied. It should be mentioned that, during the first few equilibrium iterations, the values of the element axial nodal forces obtained from the current deformation using eqns (6) and (29) may be several orders larger than their convergent values for certain problems. This may cause di~culty in convergence or even divergence for a large increment. In [IO] a two-cycle iteration scheme is introduced to overcome this difficulty. This scheme was proven to be very effective by numerous examples studied in [IO]; however, in the present study, it is found that the accuracy of convergent solutions obtained by this scheme is not sufficient for some problems. This difficulty probably arises due to the fact that the difference between the values of convergent element axial forces at the first and second cycles is not small for some problems. In order to overcome this difficulty an n-cycle iteration scheme, an extension of the two-cycle iteration scheme, is proposed in this study and can be described as follows. Let fdenote the element axial force corresponding to the current deformation, h denote the convergent element axial force of thejth cycle, and& denote the convergent element axial force of the fth increment. For the iterations of thejth cycle the element axial force F,, (eqn (29)) is replaced by
where f, and 4_ , are m x 1 column matrices containing the convergent values of element axial forces 4 andA_, , m is the number of the elements used for the discretization of the structure, and p, is a prescribed parameter. If eqn (41) is not satisfied, the iterations for next cycle are performed. Otherwise, the final cycle of iterations is carried out until eqn (32) is satisfied, At the final cycle, the element axial force f calculated from the current defo~ation is used in eqns (29) and (23) to obtain the element axial and bending forces. The convergent solution of the final cycle is used as the solution of the corresponding increment. 6. NUMERICAL STUDIES
Example 1. Cantilever beam with an end moment This example considered is a cantilever beam sub. jetted to a concentrated moment at the free end as shown in Fig. 5. The beam was d&ret&d by 10 elements. The results shown in Fig. 5 are obtained by
using only three increments. The number of iterations is about six per increment. As can be seen, the agreement with analytical solution is quite good. It should be mentioned that the two-cycle iteration scheme is used for this example because the element axial forces are zero for this problem. This example is extensively studied in the literature to demonstrate the efficiency of numerical methods and the large rotation capability of the beam, plate and shell elements. To the authors’ knowledge, only the present authors have achieved bending of the cantilever beam into a full circle by using only three increments. Example 2. Cantilever 4Megree bend with an end Iwd The bend as illustrated in Fig. 6 is curved in the horizontal pfane and subjected to a vertical load. The
+=====&=I”F W
F;,==(l-CJJ;_,+CJ
j21,
(39)
E = l.2X10skN/m2 V=O L = 10.0 m b = 1.0 m h = 0.1 m
where Ca f [0,11, is a prescribed parameter; the FL2 required for the evaluation of $ in eqn (18) and Fb in eqn (23) are replaced by J$=
“fit*
j=l
(I--C&,+C&,,
ja2,
I
(40)
where Cbe [0, I] is a prescribed parameter. The equilibrium iterations of the jth cycle are performed until the Euclidean norm of the unbalanced force vector in eqn (31) is smaller than a prescribed value, which may be chosen to be larger than the error tolerance given in eqn (32). Then the following inequality is checked:
I
I
-Analytical solution Present analysis l
12 6 8 10 Displacement (ml Fig. 5. Cantilever beam with an end moment. 0
2
4
Kuo-MO HSIAOet (11.
I2
R = 100. in u = 0. E = 1O’psi
t
c5jgft$;* 69.26
BEAM
l-l”-1 CROSS SECllON
+
I 61.44
E=4.32XlO’ (GJ), -
5 & .6
Ref *
t
[ 6 ]
Present
69.26
A,=.500
=4.15x105
(GJ),=1.66X10s
analysis
I
L=2OO.Oh A,=.100
(I,),
=.400
(Ix),
=.05
(ly),
=.133
(Iy)2=.05
PL2mJ,
2*o I 7---1.5
-
t/
Ref
“‘1 ;0.
1.
2. 3. 4. LOAD PARAMETER
5. k =
6. Pi’/El
7.
0.0
1.0
Example 3. Space arch frame
Figure 8 shows the structure load system and the load displacement curve. In addition to four vertical
(16.4
, 46.3
A (24.7
Fig.
7. Deformed
]
w,Lx,oJ
p”“““‘”
2.0
3.0
4.0
5.0
Fig. 8. Space arch frame.
Fig. 6. Forty-five-degree circular bend with an end force.
bend has an average radius of 100 in. and cross section area 1 it?. The bend is idealized using eight equal beam elements. At each increment, the stiffness matrix updating is only performed at the first five iterations of each cycle. Only four increments are used in this analysis. Three cycles of iteration are used per increment. The average number of iterations per increment is about 12. The results for the end displacements versus applied load are shown in Fig. 6 together with the solutions given by Bathe and Bolourchi [6] using eight beam elements and 60 equal load increments. Very good agreement between these two solutions is observed. The deformed configurations of the bend at various load levels are shown in Fig. 7.
[ 6
loads P, the structure is subjected to two lateral loads equal to 0.001 P. For all members the major principal axis of inertia x is normal to the plane of either arch rib. The symbols GJ, 1, and 1, denote, respectively, the torsional rigidity, the major and the minor principal moments of inertia, and the subscripts 1 and 2 denote the member groups as noted in the figure. Each member of the structure is idealized by four equal elements. The results of the present study shown in Fig. 8 are obtained by using two increments with the error tolerance p,,,, = lo-‘. Three cycles of iteration are used for both increments and the total number of iterations used for each increment are four and five, respectively. The present results are in excellent agreement with the solutions given in [8] which are obtained using 17 equal load increments (transcribed by the authors). Case A
:
Case
:
B
6 boundary nodes free in trondatlonal movsrnent 6 boundary nodes restrained against tronslotlonol movsment
, 53.6)
, 60.6
, 35.6)
shapes of a 45-degree circular bend.
-439600 lb/in2 =159000 lb/in2 10.494 iI+ PO.02 inz -0.02 it? =0.0331 in’
Fig. 9. Geometry of 12-member hexagonal frame.
Corotational
179
procedure for rotations of beam structures
& 1
V
I
x ,,w
2-
x,,u
E-3.03X1 0’ N/cm* G=l.O96Xl@ N/cd
J
5
DEfTLECllON
Fig. 10. Load-deflection
V (in)
curves for hexagonal frame.
Example 4. Twelve -member hexagonal frame
The hexagonal frame depicted in Fig. 9 is subjected to a concentrated load at the crown; two boundary conditions are considered: (a) six boundary nodes are free in translational movement and (b) six boundary nodes are restrained against translational movement. The frame was idealized using 12 (one element per member) and 36 (three equal elements per member) beam elements. For both boundary conditions, the results (not shown) using 12 elements are in close agreement with the solutions of Meek and Tan [9], who used a similar element but did not mention the number of elements used for discretixation. The results of case (a) using 36 elements are shown in Fig. 10, together with those reported by Meek and Tan [9] and Papadrakakis [7]. It is observed that the present results are in agreement with that of Papadrakakis, which is in exact agreement with experimental results given by Griggs [33]. The present results are obtained using five increments; the number of cycles used is about four per increment, and the
250-
Case B o Preset7 t -
Ref [S]
200-
Fig. 11. Load-deflection
curves for hexagonal frame.
Fig. 12. Geometry of 24-member shallow dome.
total number of iterations used per increment is about 11. The stiffness matrix is updated only at the first two iterations of each cycle. The results of case (b) using 36 elements are shown in Fig. 11. The present results are obtained using six increments; the number of cycles used per increment is about four, and the total number of iterations used per increment is about 13. The dashed curve shown in Fig. 11 is also obtained by the present study using 16 increments. The stiffness matrix is updated only at the first two iterations of each cycle. Also shown in Fig. 11 are the solutions reported by Meek and Tan [9]. The discrepancies between these two solutions may be explained by suggesting that the number of elements used in [9] is insufficient. Example 5. Twenty-four-member shaped shallow dome
hexagonal
star-
Figure 12 shows the geometry of a 24member hexagonal star-shaped shallow dome. The supports of the dome are assumed to be pinned and restrained against translational motion. The dome is idealized using 24 (one element per member) and 72 (three equal elements per member) beam elements. For all loading conditions, the results (not shown) using 12 elements are in close agreement with the solutions of Meek and Tan[9] who used a similar element but did not mention the number of elements used for discretixation. The first loading condition considered is that of a concentrated vertical load at the apex of the dome. The present results using 72 elements, shown in Fig. 13, are obtained by using three increments. The number of iterations for each cycle is also given in parenthesis beside each point on the graph. The average number of iterations required for one increment is about eight. The stiffness matrix is only updated at the Crst two iterations at each cycle. By keeping the same member cross sectional area but decreasing the tlexural stiffness in the vertical plane, the structure is reanalyzed. The graph of the load-
780
Kuo-MO HSMO er al. 2.0
2.5? .
Present
-
2.0
(I 2.1, I ) 0
I 0
Present
-Ref
Rsf [ 9 ]
[9]
l
A
P*V
OK
. Loaded node
Fig.
1.0
2.0 DEFIJXTION
13. Load~efl~tion
3.0 V (cm)
co
curves for con~ntrat~ Ioad.
central
2.0 0
Prsssnt
-Rsf
[S]
1.5 *
a0.
12=2.377cm' 13=0.295cm4
9 s
/ l
Loaded node
/
“‘l/---j Of
Fig.
’
J PO.91 Bcm’
1.0,
1.0 2.0 DEFLECTION
14. ~ad~efl~tion
curves load.
3.0
‘4.0
5.0
V (cm)
for ~n~trat~
central
2.0
0 -
curves
12=2.377cm4 13=0.295cm4 J =0.918cm4
3.0 4.0 V (cm)
3.0
for unsymmetrical loading.
5.0
deflection curves using 72 elements is shown in Fig. 14. As can be seen, only three increments are used. The average number of iterations per increment is nine. The second load condition is al1 nodes loaded symmetrically. The results using 72 elements are shown in Fig. 15. Four increments are used and the
F
1.0 2.0 DEFLECTION
Fig. 16. Load-deflection
0
!wc&ed
Present Ref [ 9 ]
0
average number of iterations used per is eight. For the unsymmetrical loading shown in Fig. 16, five increments are the average number of iterations used increment.
increment condition used and is 10 per
7. CONCLUSIONS
A practical motion process of the three dimensional beam element is presented to remove the restriction of small rotations between two successive increments for large displacement and large rotation analysis of space frames using incremental-iterative methods. The nonlinear fo~uiation is based on the corotational formulation by which the major geometric nonlinearities were shown to be embodied in the coordinate transformation when forming the element assemblage. The transformation of the element coordinate system is assumed to be accomplished by a ~anslation and two successive rigid body rotations: a transverse rotation followed by an axial rotation. The element formulation is derived based on the small deflection beam theory with the inclusion of the effect of axial force in the element coordinate system. The element internal nodal forces are calculated using the total defo~ational nodal rotations. Two methods, referred to as direct method and incremental method, are proposed in this paper to calculate the total deformational rotations. Despite the fact that the formulation of the beam element is very simple, highly accurate solutions are obtained. It is believed that the use of a simple element combined with the corotational formulation and the process of element motion proposed in this paper may represent a valuable engineering tool for the solution of nonlinear spatial beam problems. REFERENCES
0’
H. Chu and R. H. Rampetsreiter, Large deflection buckling of space frames. J. Srrucr. Do., AXE St%, 2701-2722 (1972). 2. C. Oran, Tangent stillness in space frames. J. Srrucf. Diu., ASCE 99, 987-l 101 (1973). t. K.
1.0
2.0 DENCTION
Fig. 15. Load-deflection
3.0 V (cm)
4.0
5.0
6.0
curves for symmetrical loading.
Corotational
procedure for rotations of beam structures
3. T. Belytschko and B. J. Hsieh, Non-linear transient finite element analysis with convected co-ordinates. Inl. .I. Nwner. Meth. Engng 7, 255-271 (1973). 4. T. Belytschko, L. Schwer and M. J. Klein, Large displacement transient analysis of space frames. Inl. J. Numer. Meth. Engng 11, 65-84 (1977). 5. S. N. Remseth, Nonlinear static and dynamic analysis of framed structures. Camp. Struct. 10,879-897 (1979). 6. K. J. Bathe and S. Bolourchi, Large displacement analysis of th~e~imen~onal beam structures. Inr. J. Numer. Meth. Engng 14, 96I-986 (1979). 7. M. Papadrakakis, Post-buckling analysis of spatial structures by vector iteration methods. Compur. Struct. 14, 393-402 (1981). 8. R. K. Wen and J. Rahimzadeh, Nonlinear elastic frame analysis by finite element. J. Srrucf. &gng, ASCE 109, 1952-1971 (1983). 9. J. L. Meek and H. S. Tan, Geometrically nonlinear analysis of space frames by an incremental iterative technique. Compur. Meth. appl. Mech. Engng 47, 261-282 (1984). 10. K. M. Hsiao and F. Y. Hou, Nonlinear finite element analysis of elastic frames. Compur. Srruct. 26, 693-701 (1987). 11. J. H. Argyris and P. C. Dunne, On the application of the natural mode technique to small strain large displacement problems. World Congress on Finite Element Methods in Structural Mechanics, Boumemouth, October 1975. 12. G. Wempner, Finite elements, finite rotations and small strains of flexible shells. Inr. J. Solids Srrucf. 5, 117-153 (1969). 13. J. H. Argyris, P. C. Dunne, G. A. Malejannakis and E. Schelkle, A simple triangular facet shell element with application to linear and nonlinear equilibrium and eiastic stability problems. Comput. Meih. appl. Mech. Engng IO, 371-403; 11, 97-131 (1977). 14. E. Ramm, A plate/shell element for large deflections and rotations. In Forn~ulotionsand ComputationalAfgorithmsin Ftiite Element Analysis,U. S.-Germany Symp. (Edited by K. J. Bathe, J. T.Oden and W. Wun~eriic~), DO. 264-293. MIT Press. Cambridae. MA (1977). IS. i;‘. Horrigmoe and P. G.. Bergan, Nonlinear-analysis of free-form shells by flat finite elements. Comput. Meth. appl. Mech. Engng 16, I l-35 (1978). 16. T. Belytschko and L. Glaum, Application of higher order corotational stretch theories to nonlinear finite element analysis. Comput. Struct. 10, 175-182 (1979). 17. K. J. Bathe and S. Bolourchi, A geometric and material nonlinear plate and shell element. Compur. Srrucr. 11, 23-48 (1980). 18. K. J. Bathe and L. W. Ho, A simple and effective element for analysis of general shell structures. Comput. Srrucr 13, 673-681 (1981).
781
19. T. J. R. Hughes and W. K. Liu, Nonlinear finite element analysis of shells: Part I. Three-dimensional shells; Part II. Two-dimensional shells. Comput. Meth. appl. Mech. Engng 26, 331-362 (1981); 27, 167-181 (1981). 20. J. H. Argyris et of., Finite element method-the natural approach. Comput. Meth. appi. Mech. Engng 17/l& l-106 (1979). 21. J. H. A&is, An excursion into large rotations. Cornmu. Meth. aDDI.Mech. Enpna 32. 85-155 (1982). 22. K. S: Surana, GiometticaIIy n&l&ear formulation for the curved shell elements. Inr. J. Nwner. Meth. Engng 19, 581-615 (1983). 23. C. C. Rankin and F. A. Brogan, An element independent corotational procedure for the treatment of large rotations. 3. Press. Vessel Technol., ASME 108,165-174 (1986). 24. K. M. Hsiao, Nonlinear analysis of general shell structures by flat triangular shell element. Camp. Struct. 25, 665-675 (I 987). 25. J. H. Argyris, S. Kelsey and H. Kamel, Matrix methods of structural analysis: a p&is of recent developments, In Matrix Methods of Structural Analysis. AGARDogruph 72 (Edited by B. F. de Veubeke), pp. l-165. Pergamon Press, London (1964). 26. K. Mattiasson and A. Samuelsson. Total and undated Lagrangian forms of the co-rotational finite eiement fo~ulation in g~rnet~~Ily and materially nonlinear analysis. In Nu&rical Meth&for Non-line& Problems, Vol. 2 (Edited bv C. Tavlor. E. Hinton and D. R. J. Gwen), ‘pp. 13Ll51. PiGeridge Press, Swansea, U.K. (1984). 27, M. A. Crisfield, A fast incremental/iterative solution procedure that handles “snap-through”. Compul. Sfrucl. 13, 55452 (1981). 28. T. Y. Chang and K. Sawamiphakdi, Large deflection and post-buckling analysis of shell structure. Comput. Meth. appl. Mech. Engng 32, 31 l-326 (1982). 29. D. J. Dawe, Matrix and Finite Element Displacement Analysisof Structures. Oxford University Press, London (1984). 30. H. Goldstein, Classical Mechanjcs, 2nd Edn. Addison-W~lev. Readinn. MA 11980). 31. A. K. Noor &d J. My Peters, Tracing post-hmit point with reduced basis technique. Comput.Meth. appi. Mech. Engng Zs, 217-240 (1981). 32. P. G. Etergan and T. Ssreide, Solution of large displacement and instability problems using the current stiffness parameter. In Finite Elements in Non-linear Mechanics (Edited by P. G. Bergan er al.), pp. 647-669. Tapir Press. Trondheim. Norwav (19781. 33. H. fi. Griggs, Eiperimekal stidy of instability in elements of shallow space frames. Research Report, Dept. of Civil Engineering, MIT, Cambridge, MA (1966).