- Email: [email protected]
Nonlinear analysis of plane frames using rigid body-spring discrete element method
Computers and Structures 71 (1999) 105±119
Nonlinear analysis of plane frames using rigid body-spring discrete element method Wei-Xin Ren a,*, Xiangguang Tan b, Zhaochang Zheng c a
Department of Civil Engineering, Changsha Railway University, Changsha, 410075, People's Republic of China Department of Aerospace Engineering, Mechanics and Engineering Science, University of Florida, Gainsville, FL 32611, USA c Department of Engineering Mechanics, Tsinghua University, Beijing, 100084, People's Republic of China
b
Received 26 January 1997; accepted 3 October 1998
Abstract Based on the rigid body-spring element discrete model, this paper presents the large displacement and elastic± plastic incremental formulation to analyze the ultimate load-carrying capacity of plane framed structures. A given structure is divided into a number of rigid body ®nite elements mutually connected by spring systems between elements. In such a discrete model, displacements of an element can be completely described by the rigid body motions of its centroid, while the deformation energy of the structure is stored only in the spring systems. The detailed tangential stiness matrix for plane frames has been derived under a global coordinate system. The elastic± plastic spring coecient matrix is also developed in terms of the elastic±plastic incremental theory and the internal force yielding interaction surface equation. An ecient numerical procedure is established by implementing the plastic hinge concept. The formulation has been applied to a variety of nonlinear problems of plane frames involving large displacements, large rotations but small strains and elastic±plasticity. Results obtained from this approach agree with independent analytical and other published ®nite element solutions. Numerical results show that the formulation is considerably eective and time saving. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Frame structures; Geometric nonlinearity; Elastic±plasticity; Ultimate load-carrying capacity; Rigid body-spring element; Discrete element method
1. Introduction It is well known that a structure should be designed for sucient ultimate load-carrying capacity. Due to the increasing use of limit state design methodologies instead of the allowable stress method, there is growing awareness of the need for a reliable, ecient method to analyze the ultimate load-carrying capacity of structures in the design of structures. The ultimate load-car-
* Corresponding author
rying capacity analysis of structures is generally involved in both ®nite deformations and elastic±plasticity. There is no closed-form solution because of the complexity of nonlinear problem > ms. In such a case the numerical method becomes essential. One of the most popular methods for ultimate load-carrying capacity analysis of structures, needless to say, is the nonlinear ®nite element method (NFEM). The most conventional ®nite element displacement method is based on the variational principle of the energy. The inter-element continuity of the displacement ®eld and the strain components is a basic
0045-7949/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 2 3 0 - 2
106
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
Nomenclature {dF } {dDu } {dDud} {dDup} Fi F {G } [k e] [K ] [Kep] K ai K si K yi Li Mi {P } {Pe } {q } {qi } {qn } Qi [se ] [s eT] sij [ST] ui u re u le+1 ni n re n le+1 X l {d } yi y re+1 y le+1 {DP } {Dq } {Du } Dx Dy Dy F
spring force increment vector relative displacement increment vector of the spring relative displacement elastic increment vector of the spring relative displacement plastic increment vector of the spring axial force of ith spring element force vector of rigid body-spring elements gradient vector of the yielding surface element stiness matrix of ®nite elements spring stiness coecient matrix of rigid body-spring elements elastic±plastic spring stiness coecient matrix of rigid body-spring elements axial spring stiness coecient of ith spring shear spring stiness coecient of ith spring bending spring stiness coecient of ith spring element length of ith rigid body-spring element bending moment of ith spring external load vector external load vector of rigid body-spring element e displacement vector of the rigid body-spring elements displacement vector of ith rigid body-spring elements total displacement vector after nth loading step is ®nished shear force of ith spring small deformation stiness matrix of rigid body-spring element e tangent stiness matrix of rigid body-spring element e elements of tangent stiness matrix of rigid body-spring element e global tangent stiness matrix of the structure x-axis translation component of ith rigid body-spring element x-axis translation component at the right end of rigid body-spring element e x-axis translation component at the left end of rigid body-spring element e+1 y-axis translation component of ith rigid body-spring element y-axis translation component at the right end of rigid body-spring element e y-axis translation component at the left end of rigid body-spring element e+1 displacement vector of rigid body-spring elements plastic deformation magnitude displacement vector of the ®nite element plane rotation angle of ith rigid body-spring element plane rotation angle at the right end of rigid body-spring element e plane rotation angle at the left end of rigid body-spring element e+1 external force increment vector of the structure displacement increment vector of the rigid body-spring element relative displacement vector of the spring x-axis relative translation component between rigid body-spring element e and e+1 y-axis relative translation component between rigid body-spring element e and e+1 relative rotation angle between rigid body-spring element e and e+1 internal force yielding surface function
requirement for element shape functions. Practically, the incompatible elements are commonly used in structural analysis. In order to improve FEM solutions, higher order shape functions and small elements are needed. However, the implementation of such treatments in nonlinear structural analysis requires a large amount of computer time.
NFEM formulations are generally obtained by ®nite deformation and elastic±plasticity incremental theory. The ®nite deformation theory rigorously distinguishes the reference con®guration of a deformed body. The geometric nonlinearity of three-dimensional beam elements with large displacements, large rotation and small strains was concisely studied by Bathe and
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
Bolourchi [1]. Based on choosing the reference con®guration, some nonlinear ®nite element formations were proposed, such as: . total Lagrangian formulation (TL); . updated Lagrangian formulation (UL); . co-rotational formulation. Among these formulations, the tangent stiness matrix of the elements or structure is composed of two or three separate matrices. In TL formulation, for example, the tangent stiness matrix consists of the sum of an elastic stiness matrix, an initial stress stiness matrix and an initial displacement stiness matrix. The recalculation of the element tangent stiness matrices for the nonlinear analysis requires considerable numerical integration and computation time. In this paper, a relatively simple numerical treatment method called the rigid body-spring discrete element method (RBSDEM) is developed. The method presented here is based on the rigid body-spring model (RBSM) proposed by Kawai and Toi [2] and rigid body dynamics. A given structure is divided into a proper number of rigid body elements mutually connected by spring systems between elements. In such a discrete model, displacements of each element can be completely described by the rigid body motions of its centroid so that the lowest order (zeroth order) shape function of the element is used. The deformation energy of the structure is stored only in the spring systems between elements. These incompatible rigid body elements are de®ned with concentrated stiness in the spring system between elements. The discrete element method based on the rigid body-spring model (RBSM) has been successfully applied to dierent branches of structural analysis,
107
such as the initiation and propagation cracks [3], the seismic response analysis of structures [4] and the coupled dynamic analysis of deep ocean pipes and risers [5,6]. Further studies and applications of the rigid body-spring discrete element method are still an ongoing project, attracting many researchers. This paper is intended to develop the large displacement and elastic±plasticity formulation of the rigid body-spring discrete element method to analyze the ultimate load-carrying capacity of plane frame structures involved in both geometric and material nonlinearity. The general tangent stiness matrix and incremental equilibrium equations are established under a global coordinate system. For the elastic±plastic analysis, the elastic±plastic spring coecient matrix of the rigid body-spring element is derived in terms of the elastic± plastic incremental theory and the internal force yielding interaction surface equation combining the concept of plastic hinges. Based on this formulation, the computer programs are coded where the solutions are obtained through the use of the Newton±Raphson iteration technique. The formulation can easily be applied to a variety of nonlinear static analyses of frame structures involved in large displacements, large rotations but small strains and elastic±plasticity. 2. Rigid body-spring discrete element formulation of plane frames 2.1. Spring stiness coecients The plane deformation problem of frame structures is discussed here. Within the rigid body-spring discrete element formulation the frame structure is divided into several rigid body-spring elements. The eth element is
Fig. 1. Rigid body-spring element model.
108
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
connected to the adjacent elements by an axial spring stiness K ae , a shear spring stiness K se and a bending spring stiness K ye as shown in Fig. 1. The displacements of rigid body-spring element e are represented by the rigid body motions of its centroid ue, ne, ye as shown in Fig. 2a. Furthermore, if a shear spring in the z-direction, a bending spring in y±z plane and a torsional spring with respect to the axis of the element are introduced, such rigid body-spring discrete element formulation can easily be extended to the space analysis of frame structures. In order to determine the spring stiness coecients of each spring connected to the adjacent elements, the rigid body-spring element (Fig. 2a) is compared with a corresponding ®nite element (Fig. 2b) in the local coordinate. Ke denotes the three springs between rigid body element e and rigid element e+1. Both end points of the ®nite element correspond to the centroid oe and oe+1 of rigid body element e and e+1. In the case of elasticity, the elastic potential energy stored in the spring Ke as shown in Fig. 2a can be written as follows:
1 a K ue1 ÿ ue 2 2 e 2 1 L ye1 ye K se ne1 ÿ ne ÿ 2 2
V1
1
1 K ye ye1 ÿ ye 2 2 Introducing the nodal displacement column matrix of the corresponding ®nite element as shown in Fig. 2(b): fdgT fue ,ne ,ye ,ue1 ,ne1 ,ye1 g
2
the strain energy of the ®nite element is V2
1 T e fdg k fdg 2
3
where [ke ] is a well-known element stiness matrix in FEM ke 2
EA 6 L 6 6 6 0 6 6 6 6 6 0 6 6 6 EA 6ÿ 6 6 L 6 6 6 0 6 6 4 0
3 12EI L3 6EI ÿ 2 L
4EI L
0
0
12EI L3 6EI ÿ 2 L
6EI L2 2EI L
ÿ
EA L 0 0
12EI L3 6EI L2
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 4EI 5 L
4
According to the energy equivalent principle between the elastic potential energy of the rigid body element stored in the spring and the strain energy of the ®nite element stored in the element, namely: V1 V2
5
three spring stiness coecients of spring Ke can be determined as K ae
EA L
K se
12EI L3
K ye
EI L
6
These spring stiness coecient formulas are identical to the results of a direct application of the virtual work principle [6]. However, the direct application of the virtual work principle cannot determine the shear spring stiness coecient K se. Therefore, they de®ne the shear spring stiness coecient Fig. 2. Comparison between rigid body elements and a ®nite element.
K se gK ae where g is the shear stiness factor.
7
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
2.2. The motions of rigid body-spring elements The geometric nonlinearity of frame structures described here is commonly known large displacements, large rotation but small strains. This geometric nonlinearity can be fully described by the rigid body motions of the element centroid in the following rigid body-spring discrete element formulation. The motions of rigid body-spring elements in the global coordinate±plane coordinate system o±xy is shown in Fig. 3. The initial position of the structure is located in the axis of o±x. When the loads are applied the rigid body elements move in the plane of o±xy. Assuming that the displacements of the centroid oe within rigid body element e are translation ue, ne and rotation angle ye, the displacement vector of rigid body element e can be de®ned as fqe gT fue ,ne ,ye g
8
According to the kinetic theory of rigid body, in the plane case, the translations of an arbitrary point on the rigid body-spring element e are given as follows: e
e
Å ÿR Å uÅ e uÅ eoe T ÿ I R oe
oe; I is a unit matrix; T is the coordinate transformation matrix of rigid body-spring element e, cos ye ÿsinye T sinye cosye Therefore, the displacement components at the right end point of rigid body-spring element e can be obtained: ure ue
Le cos ye ÿ 1 2
nre ne
Le sin ye 2
yre ye
10
where Le is the length of rigid body-spring element e. In the same way, the displacement components at the left end point of rigid body-spring element e+1 are ule1 ue1
Le1 cos ye1 ÿ 1 2
nle1 ne1
Le1 sin ye1 2
9
where ueoe is the translation vector of rigid body-spring element e at the centroid oe: ue, ne; Re is the position vector of an arbitrary point on rigid body elementspring e with respect to the global coordinate; ReÿReoe is the relative position vector of an arbitrary point on rigid body-spring element e with respect to the centroid
109
yle1 ye1
11
where Le+1 is the length of rigid body-spring element
Fig. 3. The motion of rigid body-spring beam elements.
110
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
e+1. The relative displacements between rigid bodyspring element e and rigid body-spring e+1 can be obtained as follows: Dx ule1 ÿ ure ue1 ÿ ue 1 1 ÿ Le1 cos ye1 Le cos ye Le1 Le 2 2
13
In the same way the force components at the left end point (spring Keÿ1) of rigid body-spring element e can be obtained: Feÿ1 K aeÿ1 ue ÿ ueÿ1 1 1 ÿ Le cos ye Leÿ1 cos yeÿ1 Le Leÿ1 2 2
Dy nle1 ÿ nre ne1 ÿ ne 1 ÿ Le1 sin ye1 Le sin ye 2
Qeÿ1 K
Dy yle1 ÿ yre ye1 ÿ ye
12
Similarly, the relative displacements between rigid body element spring e and rigid body-spring eÿ1 can also be obtained. 2.3. The forces of rigid body-spring element e In the rigid body-spring element e of Fig. 3, except the external load {Pe }, the free body diagram of rigid body-spring element e is shown as Fig. 4. From Eq. (12), the force components at the right end point (spring Ke ) of rigid body-spring element e are Fe K ae Dx K ae ue1 ÿ ue 1 1 ÿ Le1 cos ye1 Le cos ye Le1 Le 2 2
s eÿ1
Fig. 4. The free body diagram of rigid body-spring element e.
1 ne ÿ neÿ1 ÿ Le sin ye Leÿ1 sin yeÿ1 2
Meÿ1 K yeÿ1 ye ÿ yeÿ1
14
2.4. Equilibrium equations The force equilibrium equations of rigid body-spring element e are 8 9 < Feÿ1 ÿ Fe = Q ÿ Qe 15 fPe g : eÿ1 ; Meÿ1 ÿ Me where {Pe } is the external load vector acting on rigid body-spring element e, Fe, Feÿ1, Qe, Qeÿ1 are the axial forces and shear forces at the both ends of rigid bodyspring element e, the same as in Eqs. (13) and (14), but Me, Meÿ1 should be noted as some extra items are needed to meet the equilibrium requirement of bending moments: Me K ye ye1 ÿ ye ÿ Fe
Qe K se Dy K se ne1 ÿ ne 1 ÿ Le1 sin ye1 Le sin ye 2
Me K ye Dy K ye ye1 ÿ ye
Le Le sin ye Qe cos ye 2 2
Meÿ1 K yeÿ1 ye ÿ yeÿ1 Feÿ1
Le Le sin ye ÿ Qeÿ1 cos ye 2 2
The equilibrium equation (15) can be written in the matrix form of: 9 8 < qeÿ1 = se qe fPe g 16 ; : qe1 where {qi } is the displacement vector of the ith rigid body-spring element as shown in Eq. (8) and [se ] is the stiness matrix of rigid body-spring element e. Because the force components in Eq. (15) involve the sine and cosine functions of rotation angle yi, it is dicult to form the element stiness matrix [se ] explicitly unless some approximations are introduced, such as small deformations [5] or other kind of approximation [7].
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
Similar assembling procedure can be employed to form the global stiness matrix of the structure, and therefore the structural equilibrium equations become S fqg fP g
17
where [S ] is the global stiness matrix of the structure. 2.5. The tangent stiness matrix of rigid body-spring discrete elements In general, the nonlinear problems of structures should be solved by the step by step increment technique where the tangent stiness matrix is needed. In particular, to solve the material nonlinear problems, such a tangent stiness matrix is necessary. The tangent stiness matrix of rigid body-spring discrete elements is derived in the following. If the element force vector F and the displacement vector X at centroids are introduced FÅ Feÿ1 ÿ Fe ,Qeÿ1 ÿ Qe ,Meÿ1 ÿ Me T
18
Å ueÿ1 ,neÿ1 ,yeÿ1 ,ue ,ne ,ye ,ue1 ,ne1 ,ye1 T X
19
the tangent stiness matrix of rigid body-spring elements can be de®ned by ! @ F i e sTij 20 @ X j The generally explicit form of the tangent stiness matrix of rigid body-spring elements is given in Appendix A. It can be found that the element tangent stiness matrix of rigid body-spring discrete element formulation is associated with three rigid body-spring elements. This nonlinear tangent stiness matrix is concentrated only in one matrix, which diers from NFEM formulations. The tangent stiness is so general that all geometrically nonlinear sources of structural deformations, such as large displacements, large rotation but small strains, are included. Some kinds of approximated tangent stiness matrix may be obtained when the limited items of Taylor's expansions of sin yi and cos yi are taken. The most simple case is small deformation where sin yi10 and cos yi11. If it is further assumed that each rigid bodyspring element has the same section geometric properties, the same material and the same element length (i.e. each rigid body-spring element has the same spring stiness coecients), the small deformation element stiness matrix of the rigid body-spring discrete element formation can easily be obtained from the general tangent stiness matrix. The explicit form of the small deformation stiness matrix is also given in Appendix A.
111
According to Eq. (6), it is noted that the spring stiness coecients of each spring are the same in the geometrically nonlinear analysis when each rigid bodyspring element is assumed to have the same section geometric properties, the same material and the same element length. However, in the material nonlinear analysis, the spring stiness coecients of each spring may not be the same because of the material plasti®cation. In this case the spring stiness coecients of each spring must be treated dierently. 2.6. Incremental equilibrium equations Based on above tangent stiness matrix concept, the incremental equilibrium equations of rigid body-spring element e are 9 8 < Dqeÿ1 = 21 fDPe g seT Dqe ; : Dqe1 where [s eT] is the element tangent stiness matrix (see Appendix A). When the tangent stiness matrix of each element is assembled together, the global incremental equilibrium equations of the structure can be established ST fDqg fDP g
22
where {DP } is the load increment vector of the structure and [ST] is the global tangent stiness matrix of the structure. It can be seen that the assemblage procedure is simpler and such a global tangent stiness matrix has more sparseness than FEM. Incremental displacements {Dqn } at the nth loading step are solved using the Newton±Raphson iteration method. When the nth loading step is ®nished the ®nal displacements of the structure are fqn g fqnÿ1 g fDqn g
23
The above result {qn } becomes the initial value of the next loading step.
3. Elastic±plastic spring stiness coecient matrix The ultimate load-carrying capacity analysis of structures should involve both geometric and material nonlinearity. The practical analysis and design techniques ever-increasingly need structural response through the elastic±plastic range. The development of limit state design methodologies, for example, has paid particular attention on more sophisticated analysis methods than the older `allowable stress design method'.
112
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
From the literature reviews on the elastic±plastic analysis or ultimate load-carrying capacity analysis of frame structures, there are two main kinds of analysis formulations. One is the stress formulation where the stresses are the basic variables, while another is the stress-resultant formulation where the structural internal forces (stress-resultants) are the basic variables. The stress formulation is actually the plasticity ®nite element method, which is considered to be the most rigorous method. In this method structures or members are divided into elements and the element elastic± plastic stiness matrix is numerically formed according to the stresses of speci®c Gaussian integration points. In order to represent more closely the actual stress and strain distributions under consideration, those elements are sometimes further divided into segments and the sections of the segments are further divided into small blocks or layers. Although the stress formulation provides the best tool to represent various factors, such as the actual stress±strain behavior during loading and unloading as well as the strain hardening of the material, it requires huge amounts of data and large computer capacity and time. The stress-resultant formulation, on the other hand, is a simpli®ed elastic±plastic analysis method. In this method the sectional behavior of a member is assumed to be elastic±perfectly plastic, and material plasti®cation is assumed to concentrate on a zero's length zone at the ends of the members (i.e. the sections subjected to the referred internal forces). The stress-resultant formulation on the elastic±plastic analysis is sometimes called `natural approach' [8,9]. Because the internal forces are the basic variables instead of stresses, the stress-resultant formulation is always associated with the yielding equation in the form of internal forces. There are several treatments on the stress-resultant formulation, such as . classic plastic hinge method [10]; . segment stiness method based on M±P±F relationships [11]; . plastic hinge method based on internal force yielding surface [12]. The ultimate carrying load capacity analysis of frame structures involved in both geometric and material nonlinearity usually employs the stress-resultant formulation associated with the internal force yielding interactive surface equations [9,13±16]. The dierences among these analyses lie in employing dierent internal force yielding interactive surface equations. In the rigid body-spring discrete element formulation, because the zeroth order shape function is used, the internal forces naturally become the basic variables. The forces and deformations of rigid bodyspring elements are concentrated in the springs on
both ends of each element. Therefore, it is natural to apply internal force yield surface equations associated with the plastic hinge concept to the rigid body-spring discrete element formulation in elastic±plastic analysis. When the elements become elastic±plastic the element will degrade. This degradation of the element stiness is re¯ected in the change of the spring stiness coecients of the rigid body-spring elements. In such a case the elastic±plastic spring stiness coecients become essential. For the rigid body-spring element e, the internal force {F } vector of spring Ke is fF g K fDug K diag K ae ,K se ,K ye
24
where [K ] is the spring stiness coecient matrix of spring Ke, and K ae , K se, K ye are the three spring stiness coecients. {Du } is the relative spring displacement vector as de®ned in Eq. (12). According to Prandtl±Reuss ¯ow theory, increments of the relative spring displacements from a given equilibrium position may be expressed as the sum of the elastic and plastic components fdDug fdDue g fdDup g
25
Let the stress-resultant (internal force {F }) yielding condition be de®ned by the function F F F 1:0
26
By Drucker's normality rule, plastic deformations are normal to the yielding surface and orthogonal to the element force increments. Therefore, the spring force {F } in the plastic range and the plastic displacements of the spring should satisfy fdDup g fdF g 0
27
and the plastic displacement components of the spring must satisfy fdDup g lfG g
28
where l is the plastic deformation magnitude and {G } is a gradient vector of the yielding surface de®ned by T
fG g
@F @F @F , , @F @Q @M
29
Noting that l is arbitrary, one obtains fG gT fdF g 0 From Eq. (24) one may get
30
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
fdDue g K ÿ1 fdF g
31
Substituting Eq. (31) into Eq. (25), the following formulas can be obtained: fdDug K ÿ1 fdF g lfG g
32
T
113
Substituting Eq. (35) into Eq. (36), the relationship between the spring force increment {dF } and the relative spring displacement increment {dDu } is ®nally obtained fdF g Kep fdDug
37
Using {G } [K ] multiply both sides of Eq. (32), one ®nds
where [Kep] are the elastic±plastic stiness coecients of the springs
fG gT K fdDug fG gT fdF g fG gT K fG gl
Kep K ÿ K fG g fG gT K fG gÿ1 fG gT K
33
Applying Eq. (30) to Eq. (33), Eq. (33) becomes fG gT K fdDug fG gT K fG gl
34
Thus l fG gT K fG gÿ1 fG gT K fdDug
35
Again from Eq. (33) one can get fdF g K fdDug ÿ K fG gl
36
38
Using such elastic±plastic stiness coecients of the springs to modify the corresponding tangent stiness matrix of rigid body-spring elements, the elastic±plastic tangent stiness matrix of structures can easily be formed. Unlike the yielding criterion represented by stress components such as the Trasca and Mises yielding criterion, the yielding criterion (26) represented by stressresultant components depends largely on the cross section shapes of members or structures. The design inter-
Fig. 5. Large displacement and large rotation analysis of a cantilever beam.
114
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
action equations for a section have been studied and recommended [11,17]. A single-equation yielding surface for I-shape section suggested by Orbison et al. [12] is used here. The uniaxial strong axis bending interaction equation has the form F 1:15p2 m2x 3:67p2 m2x 1
39
in which p = P/Py and mx=Mx/Mpx where P = the axial force; Py=the limit axial force at whole section plasticity; Mx=the bending moment about the strong x-axis; and Mpx=the limit bending moment about the strong x-axis at whole section plasticity. Based on the above large displacement and elastic± plastic formulation of rigid body-spring discrete elements, the computer programs have been coded. Incremental equilibrium equations are solved by Newton±Raphson iteration. The formulation has been applied to a wide range of the geometric nonlinear and
ultimate load-carrying capacity analysis of plane frame structures.
4. Numerical examples 4.1. Example. The large displacement, large rotation and small strain analysis of a cantilever beam The cantilever beam with rectangle section under uniform distributed load is shown in Fig. 5a. The beam is divided into eight rigid body-spring elements. The parameters are listed as follows: L=10 cm, E=1.2 104 N/cm2, K = PL 3/EI, h=1 cm, n=0.2, b=1 cm. In the computation the load is divided into 12 steps in the form of load factor k. The iteration of Newton± Raphson is used in each load step. The linear and geo-
Fig. 6. The geometrically nonlinear analysis of a shallow arch.
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
metrically nonlinear results are shown in Fig. 5b. It can be found that the results are well agreeable with those of ADINA [18]. The numerical results have shown that the nonlinear solutions converge quickly where only two iteration calculations are generally needed. 4.2. Example. The geometric nonlinear analysis of a shallow arch The shallow arch under loading is shown in Fig. 6a. It is well known that the load±displacement response analysis of a shallow arch is a typically geometrically nonlinear problem. The example arch is modeled by 11 straight rigid body-spring elements. The computational parameters are listed as follows: a=7.3397 8 , A=0.188 cm2, L=34.0 cm, H=1.09 cm, I=0.00055 cm4, E=1. 106 N/cm2, R=133.11 cm, n=0.2.
115
The load±displacement response curves are shown in Fig. 6b. It can be seen that the results of rigid bodyspring discrete element formulation are identical to the NFEM results of TL or UL formulation [19] and show a good convergence. 4.3. Example. The ultimate load analysis of a plane steel frame The test example is a one-bay steel frame with a geometry and loading as shown in Fig. 7a. The following computation parameters are used: A=1.49 10ÿ2m2, I=2.517 10ÿ5m4, Py=3.57 106 N, Mpy=4.483 105 Nm, E=2.1 1011 N/m2, n=0.3. The ultimate load of the steel frame is solved using the present large displacement elastic±plastic formulation of rigid body-spring discrete elements. At points 1 and 4 plastic hinges occur ®rst. Consequently, the other plastic hinges at points 2 and 3 form and the
Fig. 7. The ultimate load analysis of a plane steel frame.
116
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
structure fails. The kinematic failure mechanism (mode of failure) of the frame is shown in Fig. 7b. The load± displacement curves are shown in Fig. 7c. Before the ultimate load the load±displacement relationship behaves linearly. The ultimate load is about 260 kN, which is close to the result 267 kN presented by Argyris [8]. 4.4. Example. The ultimate load-carrying capacity analysis of a plane steel frame (PÿD eect) This example is the same frame as example 4.3 but under dierent loading conditions as shown in Fig. 8a. The load±displacement curve of the frame obtained by the rigid body-spring discrete element formulation is shown in Fig. 8c. The signi®cance of the axial force of the column in the limit state condition is clearly seen
in this case where the geometric nonlinearity (PÿD eect) becomes more important. 4.5. Example. The ultimate load-carrying capacity analysis of a plane steel frame for an I-shaped section. The example plane steel frame with geometry and loading is shown in Fig. 9a. The frame is divided into 30 rigid body-spring elements. Then each member of the frame has an I-shape section. The computation parameters are: A=0.6445 10ÿ2 m2, I=1.08856 10ÿ4 m4, E=2.09 1011 N/m2, n=0.3, Py=1.77785 106 N, Mpy=2.175 105 Nm. The same frame was studied and tested by Hodge [10] where the ultimate load of the test was 133.0 kN. The load±vertical displacement and load±lateral displacement curves obtained by the rigid body-spring discrete element formulation are shown in Fig. 9b. When
Fig. 8. The ultimate load-carrying capacity analysis of a plane steel frame.
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
117
Fig. 9. The ultimate load-carrying capacity analysis of an I-section steel frame.
the load is increased to 130 kN, the plastic hinges occur at points 2 and 3. When the load is increased to 132 kN, the stiness becomes singular and the structure crashes. It is found that the current numerical results reach good agreement with Hodge's experimental results. 5. Conclusions Based on the above theoretical and numerical studies of the rigid body-spring discrete element formulation, it can be concluded that: 1. The rigid body-spring discrete element formulation employs the rigid body displacements at the centroid of each rigid body-spring element to simulate the motions of the structure. The principle of rigid body motions can be used so that the description of
structural displacements is clear. In the method the nonlinear tangent stiness of the structure is concentrated only on one tangent stiness matrix. This tangent stiness matrix has the most general form in which no approximate assumption is introduced. All geometric nonlinear sources, such as large displacements and large rotations, are clearly included. 2. One of the advantages of the rigid body-spring discrete element formulation over other discrete element methods is that it is easy to treat the elastic± plasticity. In the rigid body-spring discrete element formulation, the computation of elastic±plastic stress resultants at the both ends of the elements becomes more convenient because the forces and deformations are concentrated in the springs on both ends of each rigid body-spring element. When some elements yield, it means that the element stiness will degrade. As the tangent stiness matrix of the elements depends on the spring stiness coecients
118
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
of rigid body-spring elements, therefore, adjusting those spring stiness coecients to amend the stiness of structures coincides with the physical concept. The elastic±plastic spring stiness coecients of the rigid body-spring elements are given in this paper. Because the current elastic±plastic spring stiness coecients are based on the incremental plasticity theory and internal force yielding surface equations, the amendment of spring stiness coecients due to the material plasti®cation is clearer than that of Al-Mashary and Chen [20]. 3. Due to the simplicity of the discretion process and the sparseness of the stiness matrix, the present rigid body-spring discrete element formulation shows the superb numerical performance, especially in reducing the computational eort for the nonlinear structure analysis. The numerical studies have shown that the method is eective and agreeable with those obtained from other numerical methods and experiments. The formulation can easily be applied to a variety of nonlinear problems of framed structures involved in large displacements, large rotation but small strains and elastic± plasticity. 4. Because the elastic±plastic analysis of current rigid body-spring discrete element formulation is based on the plastic hinge concept, the formulation also belongs to an approximate method. The formulation cannot describe the graduate evolution of plastic deformation along the section and length direction of the element. However, the formulation can provide the reasonable ultimate load-carrying capacity of the structure, which is sucient to meet the requirements of engineering analysis and design.
s13
1 a K Leÿ1 sin yeÿ1 2 eÿ1
A.1. The general tangent stiness matrix of rigid bodyspring elements seT 39 2
s11 4 0 s31
0 s22 s32
s13 s23 s33
s14 0 s34
0 s25 s35
s16 s26 s36
s17 0 s37
0 s28 s38
where s11 ÿK aeÿ1
s14 K aeÿ1 K ae
s17 ÿK ae
3 s19 s29 5 s39
1 Le sin ye K aeÿ1 ÿ K ae 2
1 s19 ÿ K ae Le1 sin ye1 2 s22 ÿK seÿ1
s25 K seÿ1 K se
s28 K se
1 s23 ÿ K seÿ1 Leÿ1 cos yeÿ1 2 s26
1 Le cos ye K se ÿ K seÿ1 2
s29
1 s K Le1 cos ye1 2 e
1 s31 ÿ K aeÿ1 Le sin ye 2
s32
1 s K Le cos ye 2 eÿ1
1 s33 ÿK yeÿ1 K seÿ1 Leÿ1 Le cos yeÿ1 cos ye 4 1 K aeÿ1 Leÿ1 Le sin yeÿ1 sin ye 4
s34
1 Le sin ye K aeÿ1 ÿ K ae 2
s35
1 Le cos ye K se ÿ K seÿ1 2
s36 K yeÿ1 K ye
Appendix A. The stiness matrix of rigid body-spring elements
s16
Le a K eÿ1 cos ye ue ÿ ueÿ1 2
1 1 ÿ Le cos ye Leÿ1 cos yeÿ1 Le Leÿ1 2 2 Le 1 K ae cos ye ue1 ÿ ue ÿ Le1 cos ye1 2 2 1 Le cos ye Le1 Le 2 Le 1 K seÿ1 sin ye ne ÿ neÿ1 ÿ Le sin ye 2 2 Le Leÿ1 sin yeÿ1 K se sin ye ne1 ÿ ne 2 1 L2 ÿ Le1 sin ye1 Le sin ye e K aeÿ1 sin2 ye 2 4
L2e a 2 L2 L2 K e sin ye e K seÿ1 cos2 ye e K se cos2 ye 4 4 4
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
s37
1 a K Le sin ye 2 e
1 s38 ÿ K se Le cos ye 2 [5]
s39
1 ÿK K ae Le1 Le sin ye1 sin ye 4 y e
[6]
1 K ae Le1 Le cos ye1 cos ye 4 A.2. The small deformation stiness matrix of rigid body-spring elements s e 39 seÿ1
se
se1
[8]
where 2
ÿK a
6 6 6 0 seÿ1 6 6 4 0 2
2K a 4 0 se 0 2
ÿK
0 2K s 0
3
0
s
1 s K Le 2
ÿK a
6 6 6 0 se1 6 6 4 0
0
1 ÿ K s Le 2 ÿK y K s L 2e =4
[9]
7 7 7 7 7 5
[10] [11]
3 0 5 0 2K y K s L 2e =2 0 ÿK s
1 ÿ K s Le 2
[7]
0
1 s K Le 2 ÿK y K s L 2e =4
[12]
[13]
3 7 7 7 7 7 5
References [1] Bathe KJ, Bolourchi S. Large displacement analysis of three-dimensional beam structures. Int J Numer Meth Engng 1979;14:961±86. [2] Kawai T, Toi Y. A new element in discrete analysis of plane strain problems. Seisan Kenkyu 1977;29(4):204±10. [3] Nakazawa S, Kawai T. A rigid element spring method with application to nonlinear problems. In: . Proceedings of 1st Conference on Numerical Methods in Fracture Mechanics, Swansea, UK, 1978. p. 38±51. [4] Kawai T. New discrete models and their application to
[14]
[15] [16] [17]
[18] [19]
[20]
119
seismic response analysis of structures. Nucl Engng Des 1978;48:207±29. Mustoe GGW. A generalized formulation of the discrete element method. Engng Comput 1992;9:181±90. Mustoe GGW, Huttelmaier HP, Chung JS. Dynamic coupled bending-axial analysis of two-dimensional deepocean pipes by the discrete element method. In: . Proceedings of the Second International Oshore and Polar Engineering Conference, San Francisco, 1992. p. 504±511. Tan XG. Discrete element method to nonlinearly static and dynamic analyses of frame structures. Masters thesis, Department of Engineering Mechanics, Tsinghua University, 1995. Argyris JH. Finite element analysis of two- and threedimensional elasto±plastic framesÐthe natural approach. Comput Meth Appl Mech Engng 1982;35:221±48. Conci A. Natural approach for thin-walled beam-column with elastic-plasticity. Int J Numer Meth Engng 1990;29:1653±79. Hodge Jr PG. Plastic analysis of structure. New York: McGraw-Hill, 1959. Chen WF, Atsuta T. Theory of beam-columns. New York: McGraw-Hill, 1977. Orbison JG, Mcguire W, Abel F. Yield surface applications in nonlinear steel frame analysis. Comput Meth Appl Mech Engng 1982;33:557±73. Kondoh K, Atluri SN. Large deformation elastic-plastic analysis of frames under non-conservative loading using explicitly derived tangent stiness based on assumed stresses. Comput Mech 1987;2(1):1±25. Shi G, Atluri SN. Elasto-plastic large deformation analysis of space frame: a plastic-hinge and stress-based explicit derivation of tangent stiness. Int J Numer Meth Engng 1988;26:589±695. Chandra R. Elastic±plastic analysis of steel space structures. J Struct Engng, ASCE 1990;116(4):939±55. Toi Y, Yang HJ. Finite element crash analysis of framed structures. Comput Struct 1991;41(1):137±49. Duan L, Chen WF. Design interaction equation for steel beam±columns. J Struct Engng, ASCE 1989;115(5):1125± 243. Bathe KJ. Static and dynamic geometric and material nonlinear analysis using ADINA, User manual 1978. Mondkard DP, Power GH. Finite element analysis of non-linear static and dynamic response. Int J Numer Meth Engng 1976;11:499±520. Al-Mashary F, Chen WF. Simpli®ed second-order inelastic analysis for steel frames. Struct Engng 1991;69(23):395±9.
Nonlinear analysis of plane frames using rigid body-spring discrete element method Wei-Xin Ren a,*, Xiangguang Tan b, Zhaochang Zheng c a
Department of Civil Engineering, Changsha Railway University, Changsha, 410075, People's Republic of China Department of Aerospace Engineering, Mechanics and Engineering Science, University of Florida, Gainsville, FL 32611, USA c Department of Engineering Mechanics, Tsinghua University, Beijing, 100084, People's Republic of China
b
Received 26 January 1997; accepted 3 October 1998
Abstract Based on the rigid body-spring element discrete model, this paper presents the large displacement and elastic± plastic incremental formulation to analyze the ultimate load-carrying capacity of plane framed structures. A given structure is divided into a number of rigid body ®nite elements mutually connected by spring systems between elements. In such a discrete model, displacements of an element can be completely described by the rigid body motions of its centroid, while the deformation energy of the structure is stored only in the spring systems. The detailed tangential stiness matrix for plane frames has been derived under a global coordinate system. The elastic± plastic spring coecient matrix is also developed in terms of the elastic±plastic incremental theory and the internal force yielding interaction surface equation. An ecient numerical procedure is established by implementing the plastic hinge concept. The formulation has been applied to a variety of nonlinear problems of plane frames involving large displacements, large rotations but small strains and elastic±plasticity. Results obtained from this approach agree with independent analytical and other published ®nite element solutions. Numerical results show that the formulation is considerably eective and time saving. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Frame structures; Geometric nonlinearity; Elastic±plasticity; Ultimate load-carrying capacity; Rigid body-spring element; Discrete element method
1. Introduction It is well known that a structure should be designed for sucient ultimate load-carrying capacity. Due to the increasing use of limit state design methodologies instead of the allowable stress method, there is growing awareness of the need for a reliable, ecient method to analyze the ultimate load-carrying capacity of structures in the design of structures. The ultimate load-car-
* Corresponding author
rying capacity analysis of structures is generally involved in both ®nite deformations and elastic±plasticity. There is no closed-form solution because of the complexity of nonlinear problem > ms. In such a case the numerical method becomes essential. One of the most popular methods for ultimate load-carrying capacity analysis of structures, needless to say, is the nonlinear ®nite element method (NFEM). The most conventional ®nite element displacement method is based on the variational principle of the energy. The inter-element continuity of the displacement ®eld and the strain components is a basic
0045-7949/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 2 3 0 - 2
106
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
Nomenclature {dF } {dDu } {dDud} {dDup} Fi F {G } [k e] [K ] [Kep] K ai K si K yi Li Mi {P } {Pe } {q } {qi } {qn } Qi [se ] [s eT] sij [ST] ui u re u le+1 ni n re n le+1 X l {d } yi y re+1 y le+1 {DP } {Dq } {Du } Dx Dy Dy F
spring force increment vector relative displacement increment vector of the spring relative displacement elastic increment vector of the spring relative displacement plastic increment vector of the spring axial force of ith spring element force vector of rigid body-spring elements gradient vector of the yielding surface element stiness matrix of ®nite elements spring stiness coecient matrix of rigid body-spring elements elastic±plastic spring stiness coecient matrix of rigid body-spring elements axial spring stiness coecient of ith spring shear spring stiness coecient of ith spring bending spring stiness coecient of ith spring element length of ith rigid body-spring element bending moment of ith spring external load vector external load vector of rigid body-spring element e displacement vector of the rigid body-spring elements displacement vector of ith rigid body-spring elements total displacement vector after nth loading step is ®nished shear force of ith spring small deformation stiness matrix of rigid body-spring element e tangent stiness matrix of rigid body-spring element e elements of tangent stiness matrix of rigid body-spring element e global tangent stiness matrix of the structure x-axis translation component of ith rigid body-spring element x-axis translation component at the right end of rigid body-spring element e x-axis translation component at the left end of rigid body-spring element e+1 y-axis translation component of ith rigid body-spring element y-axis translation component at the right end of rigid body-spring element e y-axis translation component at the left end of rigid body-spring element e+1 displacement vector of rigid body-spring elements plastic deformation magnitude displacement vector of the ®nite element plane rotation angle of ith rigid body-spring element plane rotation angle at the right end of rigid body-spring element e plane rotation angle at the left end of rigid body-spring element e+1 external force increment vector of the structure displacement increment vector of the rigid body-spring element relative displacement vector of the spring x-axis relative translation component between rigid body-spring element e and e+1 y-axis relative translation component between rigid body-spring element e and e+1 relative rotation angle between rigid body-spring element e and e+1 internal force yielding surface function
requirement for element shape functions. Practically, the incompatible elements are commonly used in structural analysis. In order to improve FEM solutions, higher order shape functions and small elements are needed. However, the implementation of such treatments in nonlinear structural analysis requires a large amount of computer time.
NFEM formulations are generally obtained by ®nite deformation and elastic±plasticity incremental theory. The ®nite deformation theory rigorously distinguishes the reference con®guration of a deformed body. The geometric nonlinearity of three-dimensional beam elements with large displacements, large rotation and small strains was concisely studied by Bathe and
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
Bolourchi [1]. Based on choosing the reference con®guration, some nonlinear ®nite element formations were proposed, such as: . total Lagrangian formulation (TL); . updated Lagrangian formulation (UL); . co-rotational formulation. Among these formulations, the tangent stiness matrix of the elements or structure is composed of two or three separate matrices. In TL formulation, for example, the tangent stiness matrix consists of the sum of an elastic stiness matrix, an initial stress stiness matrix and an initial displacement stiness matrix. The recalculation of the element tangent stiness matrices for the nonlinear analysis requires considerable numerical integration and computation time. In this paper, a relatively simple numerical treatment method called the rigid body-spring discrete element method (RBSDEM) is developed. The method presented here is based on the rigid body-spring model (RBSM) proposed by Kawai and Toi [2] and rigid body dynamics. A given structure is divided into a proper number of rigid body elements mutually connected by spring systems between elements. In such a discrete model, displacements of each element can be completely described by the rigid body motions of its centroid so that the lowest order (zeroth order) shape function of the element is used. The deformation energy of the structure is stored only in the spring systems between elements. These incompatible rigid body elements are de®ned with concentrated stiness in the spring system between elements. The discrete element method based on the rigid body-spring model (RBSM) has been successfully applied to dierent branches of structural analysis,
107
such as the initiation and propagation cracks [3], the seismic response analysis of structures [4] and the coupled dynamic analysis of deep ocean pipes and risers [5,6]. Further studies and applications of the rigid body-spring discrete element method are still an ongoing project, attracting many researchers. This paper is intended to develop the large displacement and elastic±plasticity formulation of the rigid body-spring discrete element method to analyze the ultimate load-carrying capacity of plane frame structures involved in both geometric and material nonlinearity. The general tangent stiness matrix and incremental equilibrium equations are established under a global coordinate system. For the elastic±plastic analysis, the elastic±plastic spring coecient matrix of the rigid body-spring element is derived in terms of the elastic± plastic incremental theory and the internal force yielding interaction surface equation combining the concept of plastic hinges. Based on this formulation, the computer programs are coded where the solutions are obtained through the use of the Newton±Raphson iteration technique. The formulation can easily be applied to a variety of nonlinear static analyses of frame structures involved in large displacements, large rotations but small strains and elastic±plasticity. 2. Rigid body-spring discrete element formulation of plane frames 2.1. Spring stiness coecients The plane deformation problem of frame structures is discussed here. Within the rigid body-spring discrete element formulation the frame structure is divided into several rigid body-spring elements. The eth element is
Fig. 1. Rigid body-spring element model.
108
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
connected to the adjacent elements by an axial spring stiness K ae , a shear spring stiness K se and a bending spring stiness K ye as shown in Fig. 1. The displacements of rigid body-spring element e are represented by the rigid body motions of its centroid ue, ne, ye as shown in Fig. 2a. Furthermore, if a shear spring in the z-direction, a bending spring in y±z plane and a torsional spring with respect to the axis of the element are introduced, such rigid body-spring discrete element formulation can easily be extended to the space analysis of frame structures. In order to determine the spring stiness coecients of each spring connected to the adjacent elements, the rigid body-spring element (Fig. 2a) is compared with a corresponding ®nite element (Fig. 2b) in the local coordinate. Ke denotes the three springs between rigid body element e and rigid element e+1. Both end points of the ®nite element correspond to the centroid oe and oe+1 of rigid body element e and e+1. In the case of elasticity, the elastic potential energy stored in the spring Ke as shown in Fig. 2a can be written as follows:
1 a K ue1 ÿ ue 2 2 e 2 1 L ye1 ye K se ne1 ÿ ne ÿ 2 2
V1
1
1 K ye ye1 ÿ ye 2 2 Introducing the nodal displacement column matrix of the corresponding ®nite element as shown in Fig. 2(b): fdgT fue ,ne ,ye ,ue1 ,ne1 ,ye1 g
2
the strain energy of the ®nite element is V2
1 T e fdg k fdg 2
3
where [ke ] is a well-known element stiness matrix in FEM ke 2
EA 6 L 6 6 6 0 6 6 6 6 6 0 6 6 6 EA 6ÿ 6 6 L 6 6 6 0 6 6 4 0
3 12EI L3 6EI ÿ 2 L
4EI L
0
0
12EI L3 6EI ÿ 2 L
6EI L2 2EI L
ÿ
EA L 0 0
12EI L3 6EI L2
7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 4EI 5 L
4
According to the energy equivalent principle between the elastic potential energy of the rigid body element stored in the spring and the strain energy of the ®nite element stored in the element, namely: V1 V2
5
three spring stiness coecients of spring Ke can be determined as K ae
EA L
K se
12EI L3
K ye
EI L
6
These spring stiness coecient formulas are identical to the results of a direct application of the virtual work principle [6]. However, the direct application of the virtual work principle cannot determine the shear spring stiness coecient K se. Therefore, they de®ne the shear spring stiness coecient Fig. 2. Comparison between rigid body elements and a ®nite element.
K se gK ae where g is the shear stiness factor.
7
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
2.2. The motions of rigid body-spring elements The geometric nonlinearity of frame structures described here is commonly known large displacements, large rotation but small strains. This geometric nonlinearity can be fully described by the rigid body motions of the element centroid in the following rigid body-spring discrete element formulation. The motions of rigid body-spring elements in the global coordinate±plane coordinate system o±xy is shown in Fig. 3. The initial position of the structure is located in the axis of o±x. When the loads are applied the rigid body elements move in the plane of o±xy. Assuming that the displacements of the centroid oe within rigid body element e are translation ue, ne and rotation angle ye, the displacement vector of rigid body element e can be de®ned as fqe gT fue ,ne ,ye g
8
According to the kinetic theory of rigid body, in the plane case, the translations of an arbitrary point on the rigid body-spring element e are given as follows: e
e
Å ÿR Å uÅ e uÅ eoe T ÿ I R oe
oe; I is a unit matrix; T is the coordinate transformation matrix of rigid body-spring element e, cos ye ÿsinye T sinye cosye Therefore, the displacement components at the right end point of rigid body-spring element e can be obtained: ure ue
Le cos ye ÿ 1 2
nre ne
Le sin ye 2
yre ye
10
where Le is the length of rigid body-spring element e. In the same way, the displacement components at the left end point of rigid body-spring element e+1 are ule1 ue1
Le1 cos ye1 ÿ 1 2
nle1 ne1
Le1 sin ye1 2
9
where ueoe is the translation vector of rigid body-spring element e at the centroid oe: ue, ne; Re is the position vector of an arbitrary point on rigid body elementspring e with respect to the global coordinate; ReÿReoe is the relative position vector of an arbitrary point on rigid body-spring element e with respect to the centroid
109
yle1 ye1
11
where Le+1 is the length of rigid body-spring element
Fig. 3. The motion of rigid body-spring beam elements.
110
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
e+1. The relative displacements between rigid bodyspring element e and rigid body-spring e+1 can be obtained as follows: Dx ule1 ÿ ure ue1 ÿ ue 1 1 ÿ Le1 cos ye1 Le cos ye Le1 Le 2 2
13
In the same way the force components at the left end point (spring Keÿ1) of rigid body-spring element e can be obtained: Feÿ1 K aeÿ1 ue ÿ ueÿ1 1 1 ÿ Le cos ye Leÿ1 cos yeÿ1 Le Leÿ1 2 2
Dy nle1 ÿ nre ne1 ÿ ne 1 ÿ Le1 sin ye1 Le sin ye 2
Qeÿ1 K
Dy yle1 ÿ yre ye1 ÿ ye
12
Similarly, the relative displacements between rigid body element spring e and rigid body-spring eÿ1 can also be obtained. 2.3. The forces of rigid body-spring element e In the rigid body-spring element e of Fig. 3, except the external load {Pe }, the free body diagram of rigid body-spring element e is shown as Fig. 4. From Eq. (12), the force components at the right end point (spring Ke ) of rigid body-spring element e are Fe K ae Dx K ae ue1 ÿ ue 1 1 ÿ Le1 cos ye1 Le cos ye Le1 Le 2 2
s eÿ1
Fig. 4. The free body diagram of rigid body-spring element e.
1 ne ÿ neÿ1 ÿ Le sin ye Leÿ1 sin yeÿ1 2
Meÿ1 K yeÿ1 ye ÿ yeÿ1
14
2.4. Equilibrium equations The force equilibrium equations of rigid body-spring element e are 8 9 < Feÿ1 ÿ Fe = Q ÿ Qe 15 fPe g : eÿ1 ; Meÿ1 ÿ Me where {Pe } is the external load vector acting on rigid body-spring element e, Fe, Feÿ1, Qe, Qeÿ1 are the axial forces and shear forces at the both ends of rigid bodyspring element e, the same as in Eqs. (13) and (14), but Me, Meÿ1 should be noted as some extra items are needed to meet the equilibrium requirement of bending moments: Me K ye ye1 ÿ ye ÿ Fe
Qe K se Dy K se ne1 ÿ ne 1 ÿ Le1 sin ye1 Le sin ye 2
Me K ye Dy K ye ye1 ÿ ye
Le Le sin ye Qe cos ye 2 2
Meÿ1 K yeÿ1 ye ÿ yeÿ1 Feÿ1
Le Le sin ye ÿ Qeÿ1 cos ye 2 2
The equilibrium equation (15) can be written in the matrix form of: 9 8 < qeÿ1 = se qe fPe g 16 ; : qe1 where {qi } is the displacement vector of the ith rigid body-spring element as shown in Eq. (8) and [se ] is the stiness matrix of rigid body-spring element e. Because the force components in Eq. (15) involve the sine and cosine functions of rotation angle yi, it is dicult to form the element stiness matrix [se ] explicitly unless some approximations are introduced, such as small deformations [5] or other kind of approximation [7].
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
Similar assembling procedure can be employed to form the global stiness matrix of the structure, and therefore the structural equilibrium equations become S fqg fP g
17
where [S ] is the global stiness matrix of the structure. 2.5. The tangent stiness matrix of rigid body-spring discrete elements In general, the nonlinear problems of structures should be solved by the step by step increment technique where the tangent stiness matrix is needed. In particular, to solve the material nonlinear problems, such a tangent stiness matrix is necessary. The tangent stiness matrix of rigid body-spring discrete elements is derived in the following. If the element force vector F and the displacement vector X at centroids are introduced FÅ Feÿ1 ÿ Fe ,Qeÿ1 ÿ Qe ,Meÿ1 ÿ Me T
18
Å ueÿ1 ,neÿ1 ,yeÿ1 ,ue ,ne ,ye ,ue1 ,ne1 ,ye1 T X
19
the tangent stiness matrix of rigid body-spring elements can be de®ned by ! @ F i e sTij 20 @ X j The generally explicit form of the tangent stiness matrix of rigid body-spring elements is given in Appendix A. It can be found that the element tangent stiness matrix of rigid body-spring discrete element formulation is associated with three rigid body-spring elements. This nonlinear tangent stiness matrix is concentrated only in one matrix, which diers from NFEM formulations. The tangent stiness is so general that all geometrically nonlinear sources of structural deformations, such as large displacements, large rotation but small strains, are included. Some kinds of approximated tangent stiness matrix may be obtained when the limited items of Taylor's expansions of sin yi and cos yi are taken. The most simple case is small deformation where sin yi10 and cos yi11. If it is further assumed that each rigid bodyspring element has the same section geometric properties, the same material and the same element length (i.e. each rigid body-spring element has the same spring stiness coecients), the small deformation element stiness matrix of the rigid body-spring discrete element formation can easily be obtained from the general tangent stiness matrix. The explicit form of the small deformation stiness matrix is also given in Appendix A.
111
According to Eq. (6), it is noted that the spring stiness coecients of each spring are the same in the geometrically nonlinear analysis when each rigid bodyspring element is assumed to have the same section geometric properties, the same material and the same element length. However, in the material nonlinear analysis, the spring stiness coecients of each spring may not be the same because of the material plasti®cation. In this case the spring stiness coecients of each spring must be treated dierently. 2.6. Incremental equilibrium equations Based on above tangent stiness matrix concept, the incremental equilibrium equations of rigid body-spring element e are 9 8 < Dqeÿ1 = 21 fDPe g seT Dqe ; : Dqe1 where [s eT] is the element tangent stiness matrix (see Appendix A). When the tangent stiness matrix of each element is assembled together, the global incremental equilibrium equations of the structure can be established ST fDqg fDP g
22
where {DP } is the load increment vector of the structure and [ST] is the global tangent stiness matrix of the structure. It can be seen that the assemblage procedure is simpler and such a global tangent stiness matrix has more sparseness than FEM. Incremental displacements {Dqn } at the nth loading step are solved using the Newton±Raphson iteration method. When the nth loading step is ®nished the ®nal displacements of the structure are fqn g fqnÿ1 g fDqn g
23
The above result {qn } becomes the initial value of the next loading step.
3. Elastic±plastic spring stiness coecient matrix The ultimate load-carrying capacity analysis of structures should involve both geometric and material nonlinearity. The practical analysis and design techniques ever-increasingly need structural response through the elastic±plastic range. The development of limit state design methodologies, for example, has paid particular attention on more sophisticated analysis methods than the older `allowable stress design method'.
112
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
From the literature reviews on the elastic±plastic analysis or ultimate load-carrying capacity analysis of frame structures, there are two main kinds of analysis formulations. One is the stress formulation where the stresses are the basic variables, while another is the stress-resultant formulation where the structural internal forces (stress-resultants) are the basic variables. The stress formulation is actually the plasticity ®nite element method, which is considered to be the most rigorous method. In this method structures or members are divided into elements and the element elastic± plastic stiness matrix is numerically formed according to the stresses of speci®c Gaussian integration points. In order to represent more closely the actual stress and strain distributions under consideration, those elements are sometimes further divided into segments and the sections of the segments are further divided into small blocks or layers. Although the stress formulation provides the best tool to represent various factors, such as the actual stress±strain behavior during loading and unloading as well as the strain hardening of the material, it requires huge amounts of data and large computer capacity and time. The stress-resultant formulation, on the other hand, is a simpli®ed elastic±plastic analysis method. In this method the sectional behavior of a member is assumed to be elastic±perfectly plastic, and material plasti®cation is assumed to concentrate on a zero's length zone at the ends of the members (i.e. the sections subjected to the referred internal forces). The stress-resultant formulation on the elastic±plastic analysis is sometimes called `natural approach' [8,9]. Because the internal forces are the basic variables instead of stresses, the stress-resultant formulation is always associated with the yielding equation in the form of internal forces. There are several treatments on the stress-resultant formulation, such as . classic plastic hinge method [10]; . segment stiness method based on M±P±F relationships [11]; . plastic hinge method based on internal force yielding surface [12]. The ultimate carrying load capacity analysis of frame structures involved in both geometric and material nonlinearity usually employs the stress-resultant formulation associated with the internal force yielding interactive surface equations [9,13±16]. The dierences among these analyses lie in employing dierent internal force yielding interactive surface equations. In the rigid body-spring discrete element formulation, because the zeroth order shape function is used, the internal forces naturally become the basic variables. The forces and deformations of rigid bodyspring elements are concentrated in the springs on
both ends of each element. Therefore, it is natural to apply internal force yield surface equations associated with the plastic hinge concept to the rigid body-spring discrete element formulation in elastic±plastic analysis. When the elements become elastic±plastic the element will degrade. This degradation of the element stiness is re¯ected in the change of the spring stiness coecients of the rigid body-spring elements. In such a case the elastic±plastic spring stiness coecients become essential. For the rigid body-spring element e, the internal force {F } vector of spring Ke is fF g K fDug K diag K ae ,K se ,K ye
24
where [K ] is the spring stiness coecient matrix of spring Ke, and K ae , K se, K ye are the three spring stiness coecients. {Du } is the relative spring displacement vector as de®ned in Eq. (12). According to Prandtl±Reuss ¯ow theory, increments of the relative spring displacements from a given equilibrium position may be expressed as the sum of the elastic and plastic components fdDug fdDue g fdDup g
25
Let the stress-resultant (internal force {F }) yielding condition be de®ned by the function F F F 1:0
26
By Drucker's normality rule, plastic deformations are normal to the yielding surface and orthogonal to the element force increments. Therefore, the spring force {F } in the plastic range and the plastic displacements of the spring should satisfy fdDup g fdF g 0
27
and the plastic displacement components of the spring must satisfy fdDup g lfG g
28
where l is the plastic deformation magnitude and {G } is a gradient vector of the yielding surface de®ned by T
fG g
@F @F @F , , @F @Q @M
29
Noting that l is arbitrary, one obtains fG gT fdF g 0 From Eq. (24) one may get
30
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
fdDue g K ÿ1 fdF g
31
Substituting Eq. (31) into Eq. (25), the following formulas can be obtained: fdDug K ÿ1 fdF g lfG g
32
T
113
Substituting Eq. (35) into Eq. (36), the relationship between the spring force increment {dF } and the relative spring displacement increment {dDu } is ®nally obtained fdF g Kep fdDug
37
Using {G } [K ] multiply both sides of Eq. (32), one ®nds
where [Kep] are the elastic±plastic stiness coecients of the springs
fG gT K fdDug fG gT fdF g fG gT K fG gl
Kep K ÿ K fG g fG gT K fG gÿ1 fG gT K
33
Applying Eq. (30) to Eq. (33), Eq. (33) becomes fG gT K fdDug fG gT K fG gl
34
Thus l fG gT K fG gÿ1 fG gT K fdDug
35
Again from Eq. (33) one can get fdF g K fdDug ÿ K fG gl
36
38
Using such elastic±plastic stiness coecients of the springs to modify the corresponding tangent stiness matrix of rigid body-spring elements, the elastic±plastic tangent stiness matrix of structures can easily be formed. Unlike the yielding criterion represented by stress components such as the Trasca and Mises yielding criterion, the yielding criterion (26) represented by stressresultant components depends largely on the cross section shapes of members or structures. The design inter-
Fig. 5. Large displacement and large rotation analysis of a cantilever beam.
114
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
action equations for a section have been studied and recommended [11,17]. A single-equation yielding surface for I-shape section suggested by Orbison et al. [12] is used here. The uniaxial strong axis bending interaction equation has the form F 1:15p2 m2x 3:67p2 m2x 1
39
in which p = P/Py and mx=Mx/Mpx where P = the axial force; Py=the limit axial force at whole section plasticity; Mx=the bending moment about the strong x-axis; and Mpx=the limit bending moment about the strong x-axis at whole section plasticity. Based on the above large displacement and elastic± plastic formulation of rigid body-spring discrete elements, the computer programs have been coded. Incremental equilibrium equations are solved by Newton±Raphson iteration. The formulation has been applied to a wide range of the geometric nonlinear and
ultimate load-carrying capacity analysis of plane frame structures.
4. Numerical examples 4.1. Example. The large displacement, large rotation and small strain analysis of a cantilever beam The cantilever beam with rectangle section under uniform distributed load is shown in Fig. 5a. The beam is divided into eight rigid body-spring elements. The parameters are listed as follows: L=10 cm, E=1.2 104 N/cm2, K = PL 3/EI, h=1 cm, n=0.2, b=1 cm. In the computation the load is divided into 12 steps in the form of load factor k. The iteration of Newton± Raphson is used in each load step. The linear and geo-
Fig. 6. The geometrically nonlinear analysis of a shallow arch.
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
metrically nonlinear results are shown in Fig. 5b. It can be found that the results are well agreeable with those of ADINA [18]. The numerical results have shown that the nonlinear solutions converge quickly where only two iteration calculations are generally needed. 4.2. Example. The geometric nonlinear analysis of a shallow arch The shallow arch under loading is shown in Fig. 6a. It is well known that the load±displacement response analysis of a shallow arch is a typically geometrically nonlinear problem. The example arch is modeled by 11 straight rigid body-spring elements. The computational parameters are listed as follows: a=7.3397 8 , A=0.188 cm2, L=34.0 cm, H=1.09 cm, I=0.00055 cm4, E=1. 106 N/cm2, R=133.11 cm, n=0.2.
115
The load±displacement response curves are shown in Fig. 6b. It can be seen that the results of rigid bodyspring discrete element formulation are identical to the NFEM results of TL or UL formulation [19] and show a good convergence. 4.3. Example. The ultimate load analysis of a plane steel frame The test example is a one-bay steel frame with a geometry and loading as shown in Fig. 7a. The following computation parameters are used: A=1.49 10ÿ2m2, I=2.517 10ÿ5m4, Py=3.57 106 N, Mpy=4.483 105 Nm, E=2.1 1011 N/m2, n=0.3. The ultimate load of the steel frame is solved using the present large displacement elastic±plastic formulation of rigid body-spring discrete elements. At points 1 and 4 plastic hinges occur ®rst. Consequently, the other plastic hinges at points 2 and 3 form and the
Fig. 7. The ultimate load analysis of a plane steel frame.
116
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
structure fails. The kinematic failure mechanism (mode of failure) of the frame is shown in Fig. 7b. The load± displacement curves are shown in Fig. 7c. Before the ultimate load the load±displacement relationship behaves linearly. The ultimate load is about 260 kN, which is close to the result 267 kN presented by Argyris [8]. 4.4. Example. The ultimate load-carrying capacity analysis of a plane steel frame (PÿD eect) This example is the same frame as example 4.3 but under dierent loading conditions as shown in Fig. 8a. The load±displacement curve of the frame obtained by the rigid body-spring discrete element formulation is shown in Fig. 8c. The signi®cance of the axial force of the column in the limit state condition is clearly seen
in this case where the geometric nonlinearity (PÿD eect) becomes more important. 4.5. Example. The ultimate load-carrying capacity analysis of a plane steel frame for an I-shaped section. The example plane steel frame with geometry and loading is shown in Fig. 9a. The frame is divided into 30 rigid body-spring elements. Then each member of the frame has an I-shape section. The computation parameters are: A=0.6445 10ÿ2 m2, I=1.08856 10ÿ4 m4, E=2.09 1011 N/m2, n=0.3, Py=1.77785 106 N, Mpy=2.175 105 Nm. The same frame was studied and tested by Hodge [10] where the ultimate load of the test was 133.0 kN. The load±vertical displacement and load±lateral displacement curves obtained by the rigid body-spring discrete element formulation are shown in Fig. 9b. When
Fig. 8. The ultimate load-carrying capacity analysis of a plane steel frame.
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
117
Fig. 9. The ultimate load-carrying capacity analysis of an I-section steel frame.
the load is increased to 130 kN, the plastic hinges occur at points 2 and 3. When the load is increased to 132 kN, the stiness becomes singular and the structure crashes. It is found that the current numerical results reach good agreement with Hodge's experimental results. 5. Conclusions Based on the above theoretical and numerical studies of the rigid body-spring discrete element formulation, it can be concluded that: 1. The rigid body-spring discrete element formulation employs the rigid body displacements at the centroid of each rigid body-spring element to simulate the motions of the structure. The principle of rigid body motions can be used so that the description of
structural displacements is clear. In the method the nonlinear tangent stiness of the structure is concentrated only on one tangent stiness matrix. This tangent stiness matrix has the most general form in which no approximate assumption is introduced. All geometric nonlinear sources, such as large displacements and large rotations, are clearly included. 2. One of the advantages of the rigid body-spring discrete element formulation over other discrete element methods is that it is easy to treat the elastic± plasticity. In the rigid body-spring discrete element formulation, the computation of elastic±plastic stress resultants at the both ends of the elements becomes more convenient because the forces and deformations are concentrated in the springs on both ends of each rigid body-spring element. When some elements yield, it means that the element stiness will degrade. As the tangent stiness matrix of the elements depends on the spring stiness coecients
118
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
of rigid body-spring elements, therefore, adjusting those spring stiness coecients to amend the stiness of structures coincides with the physical concept. The elastic±plastic spring stiness coecients of the rigid body-spring elements are given in this paper. Because the current elastic±plastic spring stiness coecients are based on the incremental plasticity theory and internal force yielding surface equations, the amendment of spring stiness coecients due to the material plasti®cation is clearer than that of Al-Mashary and Chen [20]. 3. Due to the simplicity of the discretion process and the sparseness of the stiness matrix, the present rigid body-spring discrete element formulation shows the superb numerical performance, especially in reducing the computational eort for the nonlinear structure analysis. The numerical studies have shown that the method is eective and agreeable with those obtained from other numerical methods and experiments. The formulation can easily be applied to a variety of nonlinear problems of framed structures involved in large displacements, large rotation but small strains and elastic± plasticity. 4. Because the elastic±plastic analysis of current rigid body-spring discrete element formulation is based on the plastic hinge concept, the formulation also belongs to an approximate method. The formulation cannot describe the graduate evolution of plastic deformation along the section and length direction of the element. However, the formulation can provide the reasonable ultimate load-carrying capacity of the structure, which is sucient to meet the requirements of engineering analysis and design.
s13
1 a K Leÿ1 sin yeÿ1 2 eÿ1
A.1. The general tangent stiness matrix of rigid bodyspring elements seT 39 2
s11 4 0 s31
0 s22 s32
s13 s23 s33
s14 0 s34
0 s25 s35
s16 s26 s36
s17 0 s37
0 s28 s38
where s11 ÿK aeÿ1
s14 K aeÿ1 K ae
s17 ÿK ae
3 s19 s29 5 s39
1 Le sin ye K aeÿ1 ÿ K ae 2
1 s19 ÿ K ae Le1 sin ye1 2 s22 ÿK seÿ1
s25 K seÿ1 K se
s28 K se
1 s23 ÿ K seÿ1 Leÿ1 cos yeÿ1 2 s26
1 Le cos ye K se ÿ K seÿ1 2
s29
1 s K Le1 cos ye1 2 e
1 s31 ÿ K aeÿ1 Le sin ye 2
s32
1 s K Le cos ye 2 eÿ1
1 s33 ÿK yeÿ1 K seÿ1 Leÿ1 Le cos yeÿ1 cos ye 4 1 K aeÿ1 Leÿ1 Le sin yeÿ1 sin ye 4
s34
1 Le sin ye K aeÿ1 ÿ K ae 2
s35
1 Le cos ye K se ÿ K seÿ1 2
s36 K yeÿ1 K ye
Appendix A. The stiness matrix of rigid body-spring elements
s16
Le a K eÿ1 cos ye ue ÿ ueÿ1 2
1 1 ÿ Le cos ye Leÿ1 cos yeÿ1 Le Leÿ1 2 2 Le 1 K ae cos ye ue1 ÿ ue ÿ Le1 cos ye1 2 2 1 Le cos ye Le1 Le 2 Le 1 K seÿ1 sin ye ne ÿ neÿ1 ÿ Le sin ye 2 2 Le Leÿ1 sin yeÿ1 K se sin ye ne1 ÿ ne 2 1 L2 ÿ Le1 sin ye1 Le sin ye e K aeÿ1 sin2 ye 2 4
L2e a 2 L2 L2 K e sin ye e K seÿ1 cos2 ye e K se cos2 ye 4 4 4
W.-X. Ren et al. / Computers and Structures 71 (1999) 105±119
s37
1 a K Le sin ye 2 e
1 s38 ÿ K se Le cos ye 2 [5]
s39
1 ÿK K ae Le1 Le sin ye1 sin ye 4 y e
[6]
1 K ae Le1 Le cos ye1 cos ye 4 A.2. The small deformation stiness matrix of rigid body-spring elements s e 39 seÿ1
se
se1
[8]
where 2
ÿK a
6 6 6 0 seÿ1 6 6 4 0 2
2K a 4 0 se 0 2
ÿK
0 2K s 0
3
0
s
1 s K Le 2
ÿK a
6 6 6 0 se1 6 6 4 0
0
1 ÿ K s Le 2 ÿK y K s L 2e =4
[9]
7 7 7 7 7 5
[10] [11]
3 0 5 0 2K y K s L 2e =2 0 ÿK s
1 ÿ K s Le 2
[7]
0
1 s K Le 2 ÿK y K s L 2e =4
[12]
[13]
3 7 7 7 7 7 5
References [1] Bathe KJ, Bolourchi S. Large displacement analysis of three-dimensional beam structures. Int J Numer Meth Engng 1979;14:961±86. [2] Kawai T, Toi Y. A new element in discrete analysis of plane strain problems. Seisan Kenkyu 1977;29(4):204±10. [3] Nakazawa S, Kawai T. A rigid element spring method with application to nonlinear problems. In: . Proceedings of 1st Conference on Numerical Methods in Fracture Mechanics, Swansea, UK, 1978. p. 38±51. [4] Kawai T. New discrete models and their application to
[14]
[15] [16] [17]
[18] [19]
[20]
119
seismic response analysis of structures. Nucl Engng Des 1978;48:207±29. Mustoe GGW. A generalized formulation of the discrete element method. Engng Comput 1992;9:181±90. Mustoe GGW, Huttelmaier HP, Chung JS. Dynamic coupled bending-axial analysis of two-dimensional deepocean pipes by the discrete element method. In: . Proceedings of the Second International Oshore and Polar Engineering Conference, San Francisco, 1992. p. 504±511. Tan XG. Discrete element method to nonlinearly static and dynamic analyses of frame structures. Masters thesis, Department of Engineering Mechanics, Tsinghua University, 1995. Argyris JH. Finite element analysis of two- and threedimensional elasto±plastic framesÐthe natural approach. Comput Meth Appl Mech Engng 1982;35:221±48. Conci A. Natural approach for thin-walled beam-column with elastic-plasticity. Int J Numer Meth Engng 1990;29:1653±79. Hodge Jr PG. Plastic analysis of structure. New York: McGraw-Hill, 1959. Chen WF, Atsuta T. Theory of beam-columns. New York: McGraw-Hill, 1977. Orbison JG, Mcguire W, Abel F. Yield surface applications in nonlinear steel frame analysis. Comput Meth Appl Mech Engng 1982;33:557±73. Kondoh K, Atluri SN. Large deformation elastic-plastic analysis of frames under non-conservative loading using explicitly derived tangent stiness based on assumed stresses. Comput Mech 1987;2(1):1±25. Shi G, Atluri SN. Elasto-plastic large deformation analysis of space frame: a plastic-hinge and stress-based explicit derivation of tangent stiness. Int J Numer Meth Engng 1988;26:589±695. Chandra R. Elastic±plastic analysis of steel space structures. J Struct Engng, ASCE 1990;116(4):939±55. Toi Y, Yang HJ. Finite element crash analysis of framed structures. Comput Struct 1991;41(1):137±49. Duan L, Chen WF. Design interaction equation for steel beam±columns. J Struct Engng, ASCE 1989;115(5):1125± 243. Bathe KJ. Static and dynamic geometric and material nonlinear analysis using ADINA, User manual 1978. Mondkard DP, Power GH. Finite element analysis of non-linear static and dynamic response. Int J Numer Meth Engng 1976;11:499±520. Al-Mashary F, Chen WF. Simpli®ed second-order inelastic analysis for steel frames. Struct Engng 1991;69(23):395±9.