- Email: [email protected]
Nonlinear dynamic analysis of elastic frames
NONLINEAR
DYNAMIC Kuo-MO
Department
ANALYSIS HSIAO
OF ELASTIC
and JING-YUH
FRAMES
JANG
of Mechanical Engineering, National Chiao Tung University, Hsinchu, Taiwan 30049, R.O.C.
Abstracts-A
co-rotational formulation and numerical procedure for the dynamic analysis of in-plane frames is presented. The nodal coordinates, incremental displacements and rotations, velocities, accelerations, and equations of motion of the system are defined in terms of fixed global coordinates, whereas the total strains in the beam element are measured in an element coordinate system which rotates and translates with the element but does not deform with the element. The element equations are constructed using the small deflection beam theory with the inclusion of the effect of axial force first in the element coordinate system, and then transfo~ed to the global coordinate system by standard procedure. An incremental--iterative method based on the Newmark direct integration method and the NewtonRaphson method is used here for the solution of the nonlinear dynamic equilibrium equations. To improve convergence properties of the equilibrium iterations, the previous convergent axial force is used to calculate the geometric stiffness matrix. Numerical examples show the accuracy and efficiency of the present method
I. INTRODUCTION
In recent years, the nonlinear dynamic behavior of framed structures has been the subject of considerabte research. Nonlinear dynamic response analysis requires much more computer time than does static analysis. This is because of the necessity to advance the solution at each time step over the entire desired response. Thus, several elements using many different fo~ulation strategies and procedures to accommo~ date large rotation capability during the large deformation process have been suggested to improve the accuracy and computational efficiency of nonlinear dynamic analysis [l-7]. However. the calculation of any real structure is so expensive in computer time that some methods must be eliminated as not economical. The development of efficient formulations and elements is still a serious problem. This study aims to provide an easy, e!Iicient and reliable approach for nonlinear dynamic analysis of plane frames. To this end, the formulation of beam element proposed in [8] for nonlinear static analysis is extended for nonlinear dynamic analysis of plane frames undergoing large displacements and rotations. In [S], a co-rotational formulation of beam element and numerical procedure are presented for nonlinear static analysis of elastic frames. Based on the small deflection beam theory with the inclusion of the effect of axial force, the element equations are constructed in an element coordinate system which rotates and translates with the element but does not deform with the element. and then transformed to the global coordinate system by standard procedure. This element is applied to a variety of examples of beams and frames with large displacements and rotations in [8]. and highly accurate results are obtained.
In this study. the nodal coordinates, incremental displacements and rotations, velocities, accelerations, and equations of motion of the structure are defined in terms of fixed global coordinates, whereas the stiffness, mass, and damping matrices of the individual elements are first constructed in the element coordinate system, and are then transformed to the global coordinate system by standard procedure before assemblage. The numerical algorithm used here is an incremental-iterative method based on the Newmark direct integration method and the Newton-Raphson method. To improve the convergence properties of the equilibrium iterations, the convergent axial forces of the previous increment are used to calculate the geometric stiffness matrices of the individua! elements. Numerical examples are presented and compared with the results reported in the literature to demonstrate the accuracy and efficiency of the proposed method. 2. COORDINATE SYSTEMS
One of the basic considerations in formulating nonlinear structural probIems is the selection of a description of motion. In the present study, the co-rotational (CR) formulation is adopted. In this formulation, each element is associated with a local Cartesian element coordinate system X, y (as shown in Fig. l), which rotates and translates with the element but does not deform with the element. The deformations, stiffness matrices, mass matrices and damping matrices of the individual elements are defined in terms of the element coordinates. As can be seen in Fig. I, the element coordinate system used here is a right-handed one, the origin of the element
Kuo-MO HSIAOand JING-YUHJANG
1058
Current
Similarly, the nodal velocity and acceleration vector can be transformed from the global coordinates to element coordinates using T’, the transpose of the transformation matrix T in eqn (1).
Coniiguration
3.
FINITE
ELEMENT
FORMULATION
In the following, only the basic relations of the beam element are derived. The more detailed description may be obtained from [8-lo]. 3.1. Kinematics of beam element
\
Ith
‘/
Equilibrium
L
Initial
Configuration
Configuration
I
Fig. I. Coordinate system, member deformations and associated forces.
system is located at the node 1, and the x-axis is chosen to pass through the two end nodes of the element. Also shown in Fig. 1 is a fixed global coordinate system X, Y. The nodal coordinates, incremental nodal displacements and rotations, nodal velocities and accelerations, and dynamic equilibrium equations of the system are expressed in terms of these coordinates. The element equations are first formulated in the element coordinate system, and then transformed to the global coordinate system by standard procedure before element assemblage. The element used here is the beam element proposed by Hsiao and Hou [8]. It has three degrees of freedom per node: these are the translations Au/ and AII~(j = 1,2) in the x- and y-directions, respectively, and the counterclockwise rotations AO, (j = 1,2) at nodes j. The nodal degrees of freedom for the global coordinates are translations AU, and A V, (j = 1,2) in the X- and Y-directions, respectively, and the counterclockwise rotations AtIj (j = 1,2) at nodesj. Transformation of displacement parameters at nodes j from element coordinates to global coordinates can be performed as follows:
coordinate
In this study, the deformation of an element is referred to the element coordinate system. The beam element used here is an extension of that proposed in [8], which is applicable to arbitrarily large rotations but restricted to small deformations relative to the element axis. The deformation representation is based on the kinematic assumptions that the EulerBernoulli hypothesis is valid and the membrane strain along the deformed beam axis is constant. 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. In the following, the symbol 0 is used for xderivatives and the symbol { } for a column matrix. The superscript t stands for the transpose of a matrix or a vector. If the deformed axis of the beam element is assumed to be the cubic Hermitian polynomials in the element coordinates, the transverse displacement of the beam axis may be given by u(x) = N’ub
(2)
where N is the vector of the shape functions and is given by
= ;(1-1)2(2+5),;(1 {
tcl +()2(2-i;)+-
-5?)(f
-5).
1 +52)u
-to
1 ’
(3)
where c = x2 - x, is the current chord length of the beam member; r = - 1 + 2(x - x,)/c is a nondimensional coordinate. ub is the vector of bending displacements and is given by (4)
where 0 is the angle measured counterclockwise the X-axis to the x-axis.
from
where u, and e,, j = 1, 2, are the nodal displacements and nodal rotations, respectively, at nodes j as shown in Fig. 1. Note that due to the definition of the element coordinate, the only nonzero nodal parameters are S,,j = 1,2.
1059
Nonlinear dynamic analysis of elastic frames
Making use of the second kinematic assumption we may write the membrane strain as follows:
EI Lj
k b=-
12
6L
6L -12
4L2 -6L
6L
2Lz
(5)
em= (S - SO)/% 9
where S, is the initial arc length of the beam axis and S is the current arc length of the beam axis expressed
k,=&
36
3L
-36
I
3L
4LZ
-3L
ej=e;+Aej--~,
j=
1,2.
The element matrices required for the present study are stiffness matrix k, mass matrix m, and damping matrix c. The element stiffhess matrix is formulated by superimpo~ng the bending stiffness matrix k, and geometric stiffness matrix k, of the basic beam element, and the stiffness matrix k, of the linear bar element. The element mass matrix is obtained by superimposing m, and mrr the consistent mass matrices of the basic beam element for lateral translational inertia and rotational inertia, respectively, and m,, the consistent mass matrix of the linear bar element. For this study mass proportional damping is used. The derivation of these matrices is well documented in the textbooks and thus will not be repeated here. However, these matrices are given as follows. (a) Bending stiffness matrix kn:
-3L
I * ‘lo’
4Lq
1
(ll)
-“l1’
where L is the initial length of the beam axis, &‘I is flexural rigidity, fXSis the axial nodal force at node 2, and AE is the axial rigidity. The degrees of freedom corresponding to kb and k, are I+, defined in eqn (4). The degrees of freedom corresponding to k, are
where tl, and ua are axial translations Fig. 1.
as shown in
(d) Consistent mass matrix for lateral translation m,: r
156
22L
22L 54 -13L
54
4LZ 13L -3L2
-13L7
I3L 156 -22L
-3L.z -22L 4L2
(e) Consistent mass matrix for rotation m,: 36 3L‘ 3L
(8)
3.3. Element matrices
36 -3L
L [ -1 1
k/E
(7)
where AU,, A< and At?,,j = 1,2 are the incremental nodal displacements and rotations in the global coordinates. Then the current nodal coordinates X;., & of the element are obtained by adding AUj and AP’, to .?$ and Y! respectively. Let c((Fig. 1) denote the angle of the rigid body rotation measured from the xl-axis to the current x-axis; the current deformed nodal angles are then given by
-3L -L2
3L-l -L2
(cc) Axial stiffness matrix km:
3.2. Determination of reformational rotations
AQ,=~A~~,AV~.AB,,AUz,AV,,A8~),
(9)
4Lz
-6L
r
-36 3L t
(1 + u ‘*)li2dT, 1
The procedure used here to determine the deformationat nodal rotations for individual elements is the same as that given in [8]. However, it is repeated here for completeness. Let e,‘,j = 1,2; Xj’, Y:,j =I I,2 and x’, y’shown in Fig. 1 be the total deformed nodal rotations, global nodal coordinates and the axes of the element coordinates, respectively, for a single element at the dynamic equilibrium configuration of the fth increment, and AQ, be the incremental nodal displacement vector of a single element extracted from the incremental displacement vector of the system of equations, and given by
2L2 11 -6L
(b) Geometric stiffness matrix kg:
by
where c is the current chord length of the beam member and c is given in eqn (2).
6L
-12 -6L 12
3L 4L2 -3L -L2
-36 -3L
3L -L2
36 -3L
-3L 4L” I
. WI
(f) Consistent mass matrix for axial translation
1
m-p,L2 1 m 6 [ 12’
m,:
(15)
where L is the initial length of beam axis; pA and p, are inertia constants and are defined as PA =
P
dA
and Pr =
PY’ dA,
sA
(16)
where p is the density, y is the distance from the axis of centroid to the area dA on the cross-section.
Kuo-MO HSIAOand JING-YUHJANG
1060
3.5. Equations of motion
(g) Damping matrix e: c=,um.
(17)
where /.Jis constant and m is the element mass matrix. The degrees of freedom corresponding to m, and m, are uD defined in eqn (4). The degrees of freedom corresponding to m, are u, given in eqn (12).
3.4. Ehent
nodal force sectors
The element nodal force vectors may include the inertia force vector fi, the damping force vector f,, and the internal nodal force vector due to deformation, which, for convenience, is divided into axial nodal force vector f, and bending force vector f,. The element nodal force vectors can be calculated as follows. (a) Inertia force vector f,: f,=mij,,
08)
the matrix, element mass where m is & 1 is the acceleration vector of ii,= ;iii.i:,.o;,i&,%, a single element, which is obtained by extracting the corresponding accelerations o,, 4, and bc,(j = 1,2) from the acceleration vector of the system of equations of motion first, and then transforming o,, i;;, and l?;(j = 1.2) to the element coordinate system using the transformation matrix T’, the transpose of matrix T given in eqn (1). (b) Damping force vector fd: fd= ctj,.
(1%
where c is the damping matrix given in eqn (171, & = {li, ?P,,6, * ti,, d,, t$) is the velocity vector of a single element, which is obtained by first extracting the corresponding velocities ri,, q and 4, (j = 1,2) from the velocity vector of the system of equations of motion, and then transforming ri,, r’,, and 0, u = 1,2) to the element coordinate system using T’, the transpose of matrix T given in eqn (1). (c) Bending force vector fb: fb = (kb + kg&, ,
(20)
where f, = (f,, , m, ,&, ml],& and mj V = 1,2) are shown in Fig. 1; kb, k,, and uh are given in eqns (91, (lo), and (4), respectively. (d) Axial force vector f,:
The nonlinear pressed by
equations
of motion
may be ex-
t)=Mtj+CQ+F-P=O,
(22)
where $ is the unbalanced force among the inertia force MQ, damping force CQ, internal nodal force due to deformation F, and the external nodal force P. Q and Q are the acceleration vector and velocity vector, respectively. M and C are the mass matrix and the damping matrix, respectively, for the total structure, which are assembled from element matrices. Note that the element matrices must be transformed from the element coordinate system to the global coordinate system before assemblage, by standard procedure. In this paper, a weighted Euclidean norm of the unbalanced force was used as the error measure of the dynamic equilibrium state during the equilibrium iteration. The convergence criterion for the equilibrium iteration is given by (23) where qR is a reference unbalanced force, which is chosen to be the unbalanced force obtained at the first iteration in the present study; eta, is a prescribed value of error tolerance. Unless it is stated otherwise. the error tolerance is set to 10e4 in this paper. 4. NUMERICALALGORITHM An inc~mental-iterative method based on the Newmark direct integration method [I 1] and the Newton-Raphson method is used here. The basic steps involved in the numerical solution of eqn (22) are outlined as follows. dynamic equilibrium Assume that the configuration at time t, is known. Let Q, and (5, denote the velocity and acceleration at time t,, respectively. Q,, , and Q,,+, , the velocity and acceleration at time t,+, = t, -+ At, respectively, may be obtained by the following incremental-iterative procedure. The initial incremental displacement AQ for the next time step may be chosen to be AQ = 0. Then from the Newmark method, #ai, and Q,+, can be expressed as B
_AQ w-m-pat-
(2,
1
gi
l~ ”
(24)
i
and
= AEc,{ - 1, 1),
Q,+,=Q,+AW
-~Yk+~iln+,l,
cm
(21)
where fX,, and fKz are shown in Fig. 1; AE is axial . . I_. rigidity; 6, is the membrane strain given m eqn (3).
where At is the time step size; fi and y are parameters of the Newmark method. In this study, /I = 0.25 and Y , = 0.5 are employed. _ -
Nonlinear
1
I
L/2 E = v=o
30,000
L =
20 in
b = h = P =
1 in l/8 in 2.536X10b4
I
ksi
F(t) 640
lb
I
t
lb-sec2/in4
----
of elastic frames
1061
Mondkar Present
STUDIES
Example 1. Clamped beam subjected to step load
L/2
1.00 ‘;;‘ :
analysis
5. NUMERICAL
F(t)
f ” I
dynamic
[ 3 ]
A clamped beam subjected to a central concentrated step load is analysed. The geometry of the beam and material properties are shown in Fig. 2. Because of symmetry, only half of the beam is analysed and 10 elements are used for this. The time step size is chosen to be 50 p sec. The average number of iterations per time increment is about three. The response results for central deflection are also shown in Fig. 2, together with the solution given in [3]. Very good agreement between these two solutions is observed. Example 2. Finite vibration of a cantilever beam
0
1.0
2.0
3.0
4.0
5.0
Time (msec)
Fig. 2. Dynamic response of clamped beam.
Using the displacement increment AQ and the method described in Sec. 3.2, the current configuration and deformation can be determined. Then the element nodal force vectors can be calculated in the element coordinates using the current deformation, the current acceleration and velocity vector Q,+, and Qn+,, and the method described in Sec. 3.4. Then from the assemblage of the element nodal force vectors and eqn (22), the unbalanced force $ can be calculated. If the convergence criterion in eqn (23) is not satisfied, a displacement correction R is added to the previous AQ to obtain a new incremental displacement for the next iteration. The displacement correction R may be determined by using the Newton-Raphson method as R = -K-‘$,
A cantilever beam is initially subjected to a concentrated end load. Subsequently, the load is removed, and the beam undergoes free vibration. Figure 3 shows the geometry of the beam and material properties. The beam is idealized using four elements and the time step size is set to 0.05. The average number of iterations per increment is about four. The time history of the vertical tip displacement is shown in Fig. 3. Also shown in Fig. 3 is the result reported in [7]. It is observed that the present results are in agreement with that in [7]. Example 3. Two-member frame
The two-member frame shown in Fig. 4(a) consists of two members connected by a hinge at B. Each member is discretized by six elements. The static and
EA = El =
PA = 1 S = 10 L = 10 v=o u : damping
(26)
where K is the so-called effective stiffness matrix, and may be expressed by
R-L
/IAt’
lo8 1,000
0
0.75
1.50
t
coefficient
7.0 5.6
M+&C+K,
(27)
where M, C, and K are mass matrix, damping matrix, and stiffness matrix, respectively; fl and y are the parameters of the Newmark method. This procedure is repeated until the convergence criterion is satisfied. It should be mentioned that to improve the convergence properties of the equilibrium iterations, the convergent axial forces of the previous increment are used to calculate the geometric stiffness matrices of the individual elements.
)
4.2
& E
2.8 1.4
8 0.0 = -1.4 d -2.8 -4.2 -5.6 -7.0 0
10
20
Fig. 3. Dynamic
30
40
50 60 Time
response
70
of cantilever
80
90 100
beam.
1062
Kuo-MO HSIAO and JING-YUH
JANG
E = 206.8
GPo L = 1.0 m b = 0.05 m = 0.02 m :: = 7.86X30’ kg/m3
w I--
0.6L
t
4
(b)
Time (set)
70 60
Static
Analysis
50 z? 2 40 ;
30
1
20
5
3
0 0.10
1-1
% 0
10
0.05 0.01
0.00
Cl.02
Displacement
(4
I-
f E
U(m)
0.20 0.15 0.10 0.05 0.00
0
0.01
0.02 0.03 0.04 Time (set)
0.05
0 36
F = 50 kN
d f=;:po” 2 E
-0.15 -0.20 F = 55 kN
tf I? 0” z z
0.80 0.60 0.40 0.20 0.00
I- ;;;; -0:SO -0.80
Fig. 4. (a) Two-member frame. (b) Static response of two-member frame. (c) Horizontal displacement at the loading point. (d) Transverse deformation at the midpoint of member BC. (e) Compressive axial force of element a.
ied.
in @I. The results for horizontal displacement at the loading point vs applied load are shown in Fig. 4(b).
Static analysis. The linear buckling load of the two-member frame is Fc, = 54.43 kN. The nonlinear static analysis is carried out using the method given
Dynamic analysis. The response of the frame subjetted to a step load is studied. The magnitudes of the step load considered here are 50, 55, and 60 kN. The
dynamic behavior of the two-member frame is stud-
1063
Nonlinear dynamic analysis of elastic frames time step size is set to 1O-4 set, and the average number of iterations per increment is about three.
The time history of the horizontal displacement at the loading point is shown in Fig. 4(c). It is observed that the amplitude of the horizontal displacement is increased with time for cases F = 55 and 60 kN. This phenomenon is well known as dynamic buckling. The transverse deflection at the center of the member BC is shown in Fig. 4(d). The variation of compressive axial force of element a (see Fig. 4a) with time is shown in Fig. 4(e). 6. CONCLUSIONS
A co-rotational formulation of beam element and numerical procedure for the nonlinear dynamic analysis of plane frames have been described. The incremental displacements and rotations, velocities, accelerations, and equations of motion of the system are defined in terms of fixed global coordinates, whereas the element matrices are constructed first in an element coordinate system which rotates and translates with the element but does not deform with the element, and are then transformed to the global coordinate system by standard procedure. The beam element used here is an extension of that proposed by Hsiao and Hou [8]. In the proposed approach the major geometric nonlinearities are shown to be embodied in the coordinate transformation when the element assemblage is formed. An incremental-iterative method based on the Newmark direct integration method and the Newton-Raphson method is used here for numerical studies. Three examples show that the present formulation of beam element and procedure are applicable to dynamic analysis of plane frames undergoing large
overall rotations and translations. Despite the fact that the formulation is very simple, highly accurate solutions are obtained. It is believed that the co-rotational formulation and simple element used here may represent a valuable engineering tool for the dynamic analysis of plane frames. REFERENCES
1. C. Oran and A. Kassimali,
Large defo~ations of framed structures under static and dynamic loads. Compur. Struct. 6, 539-547 (1976). 2. R. E. Nickell, Nonlinear dynamics by mode superposition. Cornput. Meth. appl. Mech. Engng 7(f), 107-129 (1976). 3. D. P. Mondkar and G. H. Powell, Finite element analysis of nonlinear static and dynamic response. In!. J. Numer. Meth. Engng 11, 499-520 (1977). 4. T. Belytschko and L. Schwer, Large displacement transient analysis of space frames. Int. J. Numer. Meth. Engng 11, 65-84 (1977).
5. S. N. Remseth, Nonlinear static and dynamic analysis of framed structures. Comput. Struct. 10, 879-897 (1979). 6. T. Y. Yang and S. Saigal, A simple element for static and dynamic response of beams with material and geometric nonlinea~ties. fnr. J. Namer. Meth. Engng 20, 851-867 (1984).
7. J. C. Simo and L. Vu-Quoc, On the dynamics of flexible beams under large overall motion-the plane case: Part I and Part 11. J. aotd. Mech. 53. 849-863 (1986). 8. K. M. Hsiao and iF’.Y. Hou, N&linear fiiite element
analysis of elastic frames. Comput. Struct. 26, 693-701 (1987).
9. K. M. Hsiao, F. Y. Hou and K. V. Spiliopoulos, Large displacement analysis of elasto-plastic frames. Comput. Struct. 28, 627-633 (1988).
10. W. Weaver, Jr. and P. R. Johnston, Structural Dynamics by Finife Elements. Prentice-Hall, Englewood Cliffs, NJ (1987). 11. K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982).
DYNAMIC Kuo-MO
Department
ANALYSIS HSIAO
OF ELASTIC
and JING-YUH
FRAMES
JANG
of Mechanical Engineering, National Chiao Tung University, Hsinchu, Taiwan 30049, R.O.C.
Abstracts-A
co-rotational formulation and numerical procedure for the dynamic analysis of in-plane frames is presented. The nodal coordinates, incremental displacements and rotations, velocities, accelerations, and equations of motion of the system are defined in terms of fixed global coordinates, whereas the total strains in the beam element are measured in an element coordinate system which rotates and translates with the element but does not deform with the element. The element equations are constructed using the small deflection beam theory with the inclusion of the effect of axial force first in the element coordinate system, and then transfo~ed to the global coordinate system by standard procedure. An incremental--iterative method based on the Newmark direct integration method and the NewtonRaphson method is used here for the solution of the nonlinear dynamic equilibrium equations. To improve convergence properties of the equilibrium iterations, the previous convergent axial force is used to calculate the geometric stiffness matrix. Numerical examples show the accuracy and efficiency of the present method
I. INTRODUCTION
In recent years, the nonlinear dynamic behavior of framed structures has been the subject of considerabte research. Nonlinear dynamic response analysis requires much more computer time than does static analysis. This is because of the necessity to advance the solution at each time step over the entire desired response. Thus, several elements using many different fo~ulation strategies and procedures to accommo~ date large rotation capability during the large deformation process have been suggested to improve the accuracy and computational efficiency of nonlinear dynamic analysis [l-7]. However. the calculation of any real structure is so expensive in computer time that some methods must be eliminated as not economical. The development of efficient formulations and elements is still a serious problem. This study aims to provide an easy, e!Iicient and reliable approach for nonlinear dynamic analysis of plane frames. To this end, the formulation of beam element proposed in [8] for nonlinear static analysis is extended for nonlinear dynamic analysis of plane frames undergoing large displacements and rotations. In [S], a co-rotational formulation of beam element and numerical procedure are presented for nonlinear static analysis of elastic frames. Based on the small deflection beam theory with the inclusion of the effect of axial force, the element equations are constructed in an element coordinate system which rotates and translates with the element but does not deform with the element. and then transformed to the global coordinate system by standard procedure. This element is applied to a variety of examples of beams and frames with large displacements and rotations in [8]. and highly accurate results are obtained.
In this study. the nodal coordinates, incremental displacements and rotations, velocities, accelerations, and equations of motion of the structure are defined in terms of fixed global coordinates, whereas the stiffness, mass, and damping matrices of the individual elements are first constructed in the element coordinate system, and are then transformed to the global coordinate system by standard procedure before assemblage. The numerical algorithm used here is an incremental-iterative method based on the Newmark direct integration method and the Newton-Raphson method. To improve the convergence properties of the equilibrium iterations, the convergent axial forces of the previous increment are used to calculate the geometric stiffness matrices of the individua! elements. Numerical examples are presented and compared with the results reported in the literature to demonstrate the accuracy and efficiency of the proposed method. 2. COORDINATE SYSTEMS
One of the basic considerations in formulating nonlinear structural probIems is the selection of a description of motion. In the present study, the co-rotational (CR) formulation is adopted. In this formulation, each element is associated with a local Cartesian element coordinate system X, y (as shown in Fig. l), which rotates and translates with the element but does not deform with the element. The deformations, stiffness matrices, mass matrices and damping matrices of the individual elements are defined in terms of the element coordinates. As can be seen in Fig. I, the element coordinate system used here is a right-handed one, the origin of the element
Kuo-MO HSIAOand JING-YUHJANG
1058
Current
Similarly, the nodal velocity and acceleration vector can be transformed from the global coordinates to element coordinates using T’, the transpose of the transformation matrix T in eqn (1).
Coniiguration
3.
FINITE
ELEMENT
FORMULATION
In the following, only the basic relations of the beam element are derived. The more detailed description may be obtained from [8-lo]. 3.1. Kinematics of beam element
\
Ith
‘/
Equilibrium
L
Initial
Configuration
Configuration
I
Fig. I. Coordinate system, member deformations and associated forces.
system is located at the node 1, and the x-axis is chosen to pass through the two end nodes of the element. Also shown in Fig. 1 is a fixed global coordinate system X, Y. The nodal coordinates, incremental nodal displacements and rotations, nodal velocities and accelerations, and dynamic equilibrium equations of the system are expressed in terms of these coordinates. The element equations are first formulated in the element coordinate system, and then transformed to the global coordinate system by standard procedure before element assemblage. The element used here is the beam element proposed by Hsiao and Hou [8]. It has three degrees of freedom per node: these are the translations Au/ and AII~(j = 1,2) in the x- and y-directions, respectively, and the counterclockwise rotations AO, (j = 1,2) at nodes j. The nodal degrees of freedom for the global coordinates are translations AU, and A V, (j = 1,2) in the X- and Y-directions, respectively, and the counterclockwise rotations AtIj (j = 1,2) at nodesj. Transformation of displacement parameters at nodes j from element coordinates to global coordinates can be performed as follows:
coordinate
In this study, the deformation of an element is referred to the element coordinate system. The beam element used here is an extension of that proposed in [8], which is applicable to arbitrarily large rotations but restricted to small deformations relative to the element axis. The deformation representation is based on the kinematic assumptions that the EulerBernoulli hypothesis is valid and the membrane strain along the deformed beam axis is constant. 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. In the following, the symbol 0 is used for xderivatives and the symbol { } for a column matrix. The superscript t stands for the transpose of a matrix or a vector. If the deformed axis of the beam element is assumed to be the cubic Hermitian polynomials in the element coordinates, the transverse displacement of the beam axis may be given by u(x) = N’ub
(2)
where N is the vector of the shape functions and is given by
= ;(1-1)2(2+5),;(1 {
tcl +()2(2-i;)+-
-5?)(f
-5).
1 +52)u
-to
1 ’
(3)
where c = x2 - x, is the current chord length of the beam member; r = - 1 + 2(x - x,)/c is a nondimensional coordinate. ub is the vector of bending displacements and is given by (4)
where 0 is the angle measured counterclockwise the X-axis to the x-axis.
from
where u, and e,, j = 1, 2, are the nodal displacements and nodal rotations, respectively, at nodes j as shown in Fig. 1. Note that due to the definition of the element coordinate, the only nonzero nodal parameters are S,,j = 1,2.
1059
Nonlinear dynamic analysis of elastic frames
Making use of the second kinematic assumption we may write the membrane strain as follows:
EI Lj
k b=-
12
6L
6L -12
4L2 -6L
6L
2Lz
(5)
em= (S - SO)/% 9
where S, is the initial arc length of the beam axis and S is the current arc length of the beam axis expressed
k,=&
36
3L
-36
I
3L
4LZ
-3L
ej=e;+Aej--~,
j=
1,2.
The element matrices required for the present study are stiffness matrix k, mass matrix m, and damping matrix c. The element stiffhess matrix is formulated by superimpo~ng the bending stiffness matrix k, and geometric stiffness matrix k, of the basic beam element, and the stiffness matrix k, of the linear bar element. The element mass matrix is obtained by superimposing m, and mrr the consistent mass matrices of the basic beam element for lateral translational inertia and rotational inertia, respectively, and m,, the consistent mass matrix of the linear bar element. For this study mass proportional damping is used. The derivation of these matrices is well documented in the textbooks and thus will not be repeated here. However, these matrices are given as follows. (a) Bending stiffness matrix kn:
-3L
I * ‘lo’
4Lq
1
(ll)
-“l1’
where L is the initial length of the beam axis, &‘I is flexural rigidity, fXSis the axial nodal force at node 2, and AE is the axial rigidity. The degrees of freedom corresponding to kb and k, are I+, defined in eqn (4). The degrees of freedom corresponding to k, are
where tl, and ua are axial translations Fig. 1.
as shown in
(d) Consistent mass matrix for lateral translation m,: r
156
22L
22L 54 -13L
54
4LZ 13L -3L2
-13L7
I3L 156 -22L
-3L.z -22L 4L2
(e) Consistent mass matrix for rotation m,: 36 3L‘ 3L
(8)
3.3. Element matrices
36 -3L
L [ -1 1
k/E
(7)
where AU,, A< and At?,,j = 1,2 are the incremental nodal displacements and rotations in the global coordinates. Then the current nodal coordinates X;., & of the element are obtained by adding AUj and AP’, to .?$ and Y! respectively. Let c((Fig. 1) denote the angle of the rigid body rotation measured from the xl-axis to the current x-axis; the current deformed nodal angles are then given by
-3L -L2
3L-l -L2
(cc) Axial stiffness matrix km:
3.2. Determination of reformational rotations
AQ,=~A~~,AV~.AB,,AUz,AV,,A8~),
(9)
4Lz
-6L
r
-36 3L t
(1 + u ‘*)li2dT, 1
The procedure used here to determine the deformationat nodal rotations for individual elements is the same as that given in [8]. However, it is repeated here for completeness. Let e,‘,j = 1,2; Xj’, Y:,j =I I,2 and x’, y’shown in Fig. 1 be the total deformed nodal rotations, global nodal coordinates and the axes of the element coordinates, respectively, for a single element at the dynamic equilibrium configuration of the fth increment, and AQ, be the incremental nodal displacement vector of a single element extracted from the incremental displacement vector of the system of equations, and given by
2L2 11 -6L
(b) Geometric stiffness matrix kg:
by
where c is the current chord length of the beam member and c is given in eqn (2).
6L
-12 -6L 12
3L 4L2 -3L -L2
-36 -3L
3L -L2
36 -3L
-3L 4L” I
. WI
(f) Consistent mass matrix for axial translation
1
m-p,L2 1 m 6 [ 12’
m,:
(15)
where L is the initial length of beam axis; pA and p, are inertia constants and are defined as PA =
P
dA
and Pr =
PY’ dA,
sA
(16)
where p is the density, y is the distance from the axis of centroid to the area dA on the cross-section.
Kuo-MO HSIAOand JING-YUHJANG
1060
3.5. Equations of motion
(g) Damping matrix e: c=,um.
(17)
where /.Jis constant and m is the element mass matrix. The degrees of freedom corresponding to m, and m, are uD defined in eqn (4). The degrees of freedom corresponding to m, are u, given in eqn (12).
3.4. Ehent
nodal force sectors
The element nodal force vectors may include the inertia force vector fi, the damping force vector f,, and the internal nodal force vector due to deformation, which, for convenience, is divided into axial nodal force vector f, and bending force vector f,. The element nodal force vectors can be calculated as follows. (a) Inertia force vector f,: f,=mij,,
08)
the matrix, element mass where m is & 1 is the acceleration vector of ii,= ;iii.i:,.o;,i&,%, a single element, which is obtained by extracting the corresponding accelerations o,, 4, and bc,(j = 1,2) from the acceleration vector of the system of equations of motion first, and then transforming o,, i;;, and l?;(j = 1.2) to the element coordinate system using the transformation matrix T’, the transpose of matrix T given in eqn (1). (b) Damping force vector fd: fd= ctj,.
(1%
where c is the damping matrix given in eqn (171, & = {li, ?P,,6, * ti,, d,, t$) is the velocity vector of a single element, which is obtained by first extracting the corresponding velocities ri,, q and 4, (j = 1,2) from the velocity vector of the system of equations of motion, and then transforming ri,, r’,, and 0, u = 1,2) to the element coordinate system using T’, the transpose of matrix T given in eqn (1). (c) Bending force vector fb: fb = (kb + kg&, ,
(20)
where f, = (f,, , m, ,&, ml],& and mj V = 1,2) are shown in Fig. 1; kb, k,, and uh are given in eqns (91, (lo), and (4), respectively. (d) Axial force vector f,:
The nonlinear pressed by
equations
of motion
may be ex-
t)=Mtj+CQ+F-P=O,
(22)
where $ is the unbalanced force among the inertia force MQ, damping force CQ, internal nodal force due to deformation F, and the external nodal force P. Q and Q are the acceleration vector and velocity vector, respectively. M and C are the mass matrix and the damping matrix, respectively, for the total structure, which are assembled from element matrices. Note that the element matrices must be transformed from the element coordinate system to the global coordinate system before assemblage, by standard procedure. In this paper, a weighted Euclidean norm of the unbalanced force was used as the error measure of the dynamic equilibrium state during the equilibrium iteration. The convergence criterion for the equilibrium iteration is given by (23) where qR is a reference unbalanced force, which is chosen to be the unbalanced force obtained at the first iteration in the present study; eta, is a prescribed value of error tolerance. Unless it is stated otherwise. the error tolerance is set to 10e4 in this paper. 4. NUMERICALALGORITHM An inc~mental-iterative method based on the Newmark direct integration method [I 1] and the Newton-Raphson method is used here. The basic steps involved in the numerical solution of eqn (22) are outlined as follows. dynamic equilibrium Assume that the configuration at time t, is known. Let Q, and (5, denote the velocity and acceleration at time t,, respectively. Q,, , and Q,,+, , the velocity and acceleration at time t,+, = t, -+ At, respectively, may be obtained by the following incremental-iterative procedure. The initial incremental displacement AQ for the next time step may be chosen to be AQ = 0. Then from the Newmark method, #ai, and Q,+, can be expressed as B
_AQ w-m-pat-
(2,
1
gi
l~ ”
(24)
i
and
= AEc,{ - 1, 1),
Q,+,=Q,+AW
-~Yk+~iln+,l,
cm
(21)
where fX,, and fKz are shown in Fig. 1; AE is axial . . I_. rigidity; 6, is the membrane strain given m eqn (3).
where At is the time step size; fi and y are parameters of the Newmark method. In this study, /I = 0.25 and Y , = 0.5 are employed. _ -
Nonlinear
1
I
L/2 E = v=o
30,000
L =
20 in
b = h = P =
1 in l/8 in 2.536X10b4
I
ksi
F(t) 640
lb
I
t
lb-sec2/in4
----
of elastic frames
1061
Mondkar Present
STUDIES
Example 1. Clamped beam subjected to step load
L/2
1.00 ‘;;‘ :
analysis
5. NUMERICAL
F(t)
f ” I
dynamic
[ 3 ]
A clamped beam subjected to a central concentrated step load is analysed. The geometry of the beam and material properties are shown in Fig. 2. Because of symmetry, only half of the beam is analysed and 10 elements are used for this. The time step size is chosen to be 50 p sec. The average number of iterations per time increment is about three. The response results for central deflection are also shown in Fig. 2, together with the solution given in [3]. Very good agreement between these two solutions is observed. Example 2. Finite vibration of a cantilever beam
0
1.0
2.0
3.0
4.0
5.0
Time (msec)
Fig. 2. Dynamic response of clamped beam.
Using the displacement increment AQ and the method described in Sec. 3.2, the current configuration and deformation can be determined. Then the element nodal force vectors can be calculated in the element coordinates using the current deformation, the current acceleration and velocity vector Q,+, and Qn+,, and the method described in Sec. 3.4. Then from the assemblage of the element nodal force vectors and eqn (22), the unbalanced force $ can be calculated. If the convergence criterion in eqn (23) is not satisfied, a displacement correction R is added to the previous AQ to obtain a new incremental displacement for the next iteration. The displacement correction R may be determined by using the Newton-Raphson method as R = -K-‘$,
A cantilever beam is initially subjected to a concentrated end load. Subsequently, the load is removed, and the beam undergoes free vibration. Figure 3 shows the geometry of the beam and material properties. The beam is idealized using four elements and the time step size is set to 0.05. The average number of iterations per increment is about four. The time history of the vertical tip displacement is shown in Fig. 3. Also shown in Fig. 3 is the result reported in [7]. It is observed that the present results are in agreement with that in [7]. Example 3. Two-member frame
The two-member frame shown in Fig. 4(a) consists of two members connected by a hinge at B. Each member is discretized by six elements. The static and
EA = El =
PA = 1 S = 10 L = 10 v=o u : damping
(26)
where K is the so-called effective stiffness matrix, and may be expressed by
R-L
/IAt’
lo8 1,000
0
0.75
1.50
t
coefficient
7.0 5.6
M+&C+K,
(27)
where M, C, and K are mass matrix, damping matrix, and stiffness matrix, respectively; fl and y are the parameters of the Newmark method. This procedure is repeated until the convergence criterion is satisfied. It should be mentioned that to improve the convergence properties of the equilibrium iterations, the convergent axial forces of the previous increment are used to calculate the geometric stiffness matrices of the individual elements.
)
4.2
& E
2.8 1.4
8 0.0 = -1.4 d -2.8 -4.2 -5.6 -7.0 0
10
20
Fig. 3. Dynamic
30
40
50 60 Time
response
70
of cantilever
80
90 100
beam.
1062
Kuo-MO HSIAO and JING-YUH
JANG
E = 206.8
GPo L = 1.0 m b = 0.05 m = 0.02 m :: = 7.86X30’ kg/m3
w I--
0.6L
t
4
(b)
Time (set)
70 60
Static
Analysis
50 z? 2 40 ;
30
1
20
5
3
0 0.10
1-1
% 0
10
0.05 0.01
0.00
Cl.02
Displacement
(4
I-
f E
U(m)
0.20 0.15 0.10 0.05 0.00
0
0.01
0.02 0.03 0.04 Time (set)
0.05
0 36
F = 50 kN
d f=;:po” 2 E
-0.15 -0.20 F = 55 kN
tf I? 0” z z
0.80 0.60 0.40 0.20 0.00
I- ;;;; -0:SO -0.80
Fig. 4. (a) Two-member frame. (b) Static response of two-member frame. (c) Horizontal displacement at the loading point. (d) Transverse deformation at the midpoint of member BC. (e) Compressive axial force of element a.
ied.
in @I. The results for horizontal displacement at the loading point vs applied load are shown in Fig. 4(b).
Static analysis. The linear buckling load of the two-member frame is Fc, = 54.43 kN. The nonlinear static analysis is carried out using the method given
Dynamic analysis. The response of the frame subjetted to a step load is studied. The magnitudes of the step load considered here are 50, 55, and 60 kN. The
dynamic behavior of the two-member frame is stud-
1063
Nonlinear dynamic analysis of elastic frames time step size is set to 1O-4 set, and the average number of iterations per increment is about three.
The time history of the horizontal displacement at the loading point is shown in Fig. 4(c). It is observed that the amplitude of the horizontal displacement is increased with time for cases F = 55 and 60 kN. This phenomenon is well known as dynamic buckling. The transverse deflection at the center of the member BC is shown in Fig. 4(d). The variation of compressive axial force of element a (see Fig. 4a) with time is shown in Fig. 4(e). 6. CONCLUSIONS
A co-rotational formulation of beam element and numerical procedure for the nonlinear dynamic analysis of plane frames have been described. The incremental displacements and rotations, velocities, accelerations, and equations of motion of the system are defined in terms of fixed global coordinates, whereas the element matrices are constructed first in an element coordinate system which rotates and translates with the element but does not deform with the element, and are then transformed to the global coordinate system by standard procedure. The beam element used here is an extension of that proposed by Hsiao and Hou [8]. In the proposed approach the major geometric nonlinearities are shown to be embodied in the coordinate transformation when the element assemblage is formed. An incremental-iterative method based on the Newmark direct integration method and the Newton-Raphson method is used here for numerical studies. Three examples show that the present formulation of beam element and procedure are applicable to dynamic analysis of plane frames undergoing large
overall rotations and translations. Despite the fact that the formulation is very simple, highly accurate solutions are obtained. It is believed that the co-rotational formulation and simple element used here may represent a valuable engineering tool for the dynamic analysis of plane frames. REFERENCES
1. C. Oran and A. Kassimali,
Large defo~ations of framed structures under static and dynamic loads. Compur. Struct. 6, 539-547 (1976). 2. R. E. Nickell, Nonlinear dynamics by mode superposition. Cornput. Meth. appl. Mech. Engng 7(f), 107-129 (1976). 3. D. P. Mondkar and G. H. Powell, Finite element analysis of nonlinear static and dynamic response. In!. J. Numer. Meth. Engng 11, 499-520 (1977). 4. T. Belytschko and L. Schwer, Large displacement transient analysis of space frames. Int. J. Numer. Meth. Engng 11, 65-84 (1977).
5. S. N. Remseth, Nonlinear static and dynamic analysis of framed structures. Comput. Struct. 10, 879-897 (1979). 6. T. Y. Yang and S. Saigal, A simple element for static and dynamic response of beams with material and geometric nonlinea~ties. fnr. J. Namer. Meth. Engng 20, 851-867 (1984).
7. J. C. Simo and L. Vu-Quoc, On the dynamics of flexible beams under large overall motion-the plane case: Part I and Part 11. J. aotd. Mech. 53. 849-863 (1986). 8. K. M. Hsiao and iF’.Y. Hou, N&linear fiiite element
analysis of elastic frames. Comput. Struct. 26, 693-701 (1987).
9. K. M. Hsiao, F. Y. Hou and K. V. Spiliopoulos, Large displacement analysis of elasto-plastic frames. Comput. Struct. 28, 627-633 (1988).
10. W. Weaver, Jr. and P. R. Johnston, Structural Dynamics by Finife Elements. Prentice-Hall, Englewood Cliffs, NJ (1987). 11. K. J. Bathe, Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ (1982).