~ompurers d Swucrwes Vol. 48, No. 3, pp. 4814W Printed in Grcat Britain.
0045.1949193s6.00 + 0.00 0 1993 Pcrgamon Press Ltd
1993
A FINITE ELEMENT PROCEDURE FOR INELASTIC DYNAMIC ANALYSIS OF CYLINDRICAL STRUCTURES b&HUNG
!?sHtAU
Institute of Aeronautics and Astronautics, National Cheng Kung University, Tainan, Taiwan, R.O.C. (Received 20 March
1992)
Al#stmet--A finite element procedure for the inelastic dynamic analysis of cylindrical structures subjected to lateral dynamic load and axial thrust is given. The plastic strains developed in the structure are treated as fictitious equivalent loads and the neutral axis at each section is determined by using a subsection concept. Equations of motion are then solved numerically by the direction integration method. The finite element procedure is verified by the exact solution using a continuum analysis. Application of the procedure is made to investigate the e&&s of axial thrust on the inelastic dynamic response of a tall cylindrical struc&m. The present method may be applied to the analysis of practical designs of cylindrical structures.
~RODU~ON Slender cylindrical structures are commonly used in industry as electric poles for transmission lines and chimneys in fossil-fuel power plants, for instance. To meet the increasing demands for air pollution control, the height of such structures needs to increase. In areas of extreme environment, the energy absorption capacity of these structures to resist severe dynamic forces, such as those induced by earthquake and wind, becomes an important factor in the design of these structures. A literature survey shows that the studies of the dynamic behavior, especially the inelastic dynamic behavior, of slender cylindrical structures subjected to moderate or severe earthquake have not received much attention. Although guidelines for estimating the dynamic forces and the corresponding displacements of such structures are provided in design codes, a more complete and detailed analysis method for inelastic dynamic response of the structures is still needed. A number of investigators[l~] studied the dynamic behavior of tall chimneys using either the mode-superposition method or the step-by-step direction integration method. However, their studies were limited to the elastic range. In the case of a strong earthquake, the stresses developed in a tall chimney may exceed its elastic limit and the dynamic behavior of the chimney may become quite different from those obtained from elastic dynamic analysis. Shiau and Yang [S] used the finite element method coupled with Wilson-8 time integration method to study the inelastic response of a tall chimney. Their results showed that plastic strains tend to resist the response motion of the chimney. When the tail, slender cylindrical structure is subjected to both Iaterat dynamic forces and axial loads its inelastic dynamic behavior may become quite different. The analysis will also become quite
complicated due to the location of the neutral axis in each section shifts with time which makes an irregular plastic strain distribution across the entire section. Toma and Chen [6] used a subsection concept and the tangent stiffness method to find the exact moment-th~st~u~ature (M-P_rb) relation of a tubular column under monotonic loading. Subsequently, they did a series of studies [7-lo} regarding the cyclic inelastic behavior of tubular columns. The purpose of this study is to develop a finite element procedure for the analysis of inelastic dynamic behavior of cylindrical structures subjected to both lateral dynamic loads and axial thrust. The plastic effect is accounted for by applying fictitious equivalent loads to the structure and the equations of motion are solved by the numerical direction integration method. FINITE ELEMENT FORMULATION
The equations of motion for an assembled system of finite elements at time t can be written as
t~lG% + IW~, f Klt~4It = ia
(1)
in which [Ml, [Cl, and [K,] are the lumped mass, damping, and tangent stiffness matrices, respectively; (q} is the vector for the nodal degrees-of-freedom; (P) is the vector for the nodal forces; t is time; and the overdot indicates a derivative with respect to time. At time t + At, eqn (1) can be written in an incremental form as
481
LE-CHUNG SHUU
482
in which {Ad},, {Ad},, and {Aq}, are the vectors of increments in acceleration, velocity, and displacement, respectively, due to the time increment Ar. Force vectors WI @I,, [Cl (QL and 1%1,(q),are known quantities in eqn (2). Equation (2) is used to find {Aq},, which, in turn, gives {Ad}, and {Ad),. The difference between elastic stiffness matrix [Kc] and tangent stiffness matrix [K,] at time t denoted as [61yl, and can be expressed as [&I - K], = The internal unbalanced divided into two parts [K,l,{q), = or in incremental
Let it
WI, *
(3)
forced at time t can be
KlbIh - wm4~r
(4)
form
Fig. 1. Tapered heam element and nodal section partition.
be de&red that U% = Wl,&IL
I@], = WUAq],
(6)
in which (F}, and (AF), are referred to as fictitious and incremental fictitious load vectors respectively at the time 2. Substituting eqns (4)-(a) into (2) gives
(7)
A circular hollow tapered beam element is used to model the cylindrical structures (Fig. la). The elastic stiffness matrix for the beam element contains three degrees-of-freedom at each nodal point and can be found in any finite element textbook. The fictitious load vector induced by the plastic strains in the element can be derived from the strain energy of the element which is obtained as u =fiqY~k,l~q~
-W&J&O,).
the initial strain stiffness matrix is assembled from many submatrices
where [kJ is the initial stiffness matrix for the strip s in Fig. l(b). It is assumed that the plastic strains within each strip of the element vary linearly in the axial direction and distribute unifo~ly in the transverse direction. Hence, the plastic strain at any point in the strip can be expressed as a function of nodal section plastic strains
The initial strain stiffness matrix [kp] can then be integrating the strains over the volume of the strip [ 111
(8)
where (fP) contains two nodal plastic strains cPi and Ed,, and [k,] is the initial strain stiffness matrix. The k is used to indicate the element stiffness matrix, while K indicates the assembled matrix for the whole system. Since the values of the element nodal degreesof-freedom remain unchanged after assemblage, {q) is used for both the individual element and the assembled system. In order to determine the plastic strain distribution at each time step, the nodal sections are divided into numerous subsections as shown in Fig. l(b). Thus,
in which [v is the usual elastic displacement relation matrix; t+ is strip in plastic range. The variable any point in the strip and can be
strain and nodal the volume of the b is the width at expressed as
483
Inelastic dynamic analysis of cylindrical structures where z,, z,, b,, and 6, are locations and widths of the top and bottom boundaries of the strip as shown in Fig. l(b). Their expressions in terms of nodal variables are the same as the one for the plastic strains or
The incremental fictitious load {Af} is calculated at each time step using the plastic strain increments determined at end of the previous time step. Solution procedure The solution
Z, = zti(, -;)
+ zti(?),
Z, = zbi( 1 _ ;) + z*j(;) (13)
of the finite element equations
using the Wilson-0 motion is determined method [12-141. Assume that {q},, {d},, and {q}, are known at time t. Also assume that acceleration varies linearly over a time interval, r = BAt, with 0 3 1.37.
bf=bt(1-f)+b(~)9 bb=bbf(l-f)+b@(~)* Thus
(~~,+.=(4z+I((1),+,+(I),)
(14) After performing the integration, stiffness matrix is derived as
4cd +fg)
[k&l= gj i
-24
0 3 -31
-61 24 -181
{qL+7= {q},+rIl},+~({~},+*+2{~~,).
+ [b(cd +fg) + c(ce + b,id +F
Mbtte+ bbih)+ b(ce + b,id +_/h + bb,g)]
51 0 -51
where [k{J is the flexural part of [k,,] and [k&J is the __ __ _ axial part of the [k,,] or
-:]+bs[-;
kl{q~ -W&J
d=2z,j+zbj-2z,i+zbi,
e =~z,~+z~~,
g=z,-+2z~-z,i-2zbi,
h=~,~+2z,,,
(18) or in an incremental
Id), - 2{4-I,. (22)
{PII,,
(23)
f=b,-bbi
(17)
~=b,~+b~~.
= {PI! + WPI,+& (q)
(19)
(16)
c = 6, - b,i
form
= {API + {Af}.
-211
eqn (22) into (20) gives
Thus, vectors Id},+1 and {q}, +T as described in eqns The equations of motion, t should also hold for time new definitions
WI),=
KI&}
-41
-WI)-22(4),-$. b = z,~- zbi,
= @I or kl{q) = (PI + 0
-61 27
-;I>,
- {4L)-;
The parameters in eqns (15) and (16) are functions of nodal variables and defined as
Fictitious louris. Performing partial differentiation of the strain energy, as defined in eqn (8) with respect to each degree-of-freedom, yields
:
Equation (21) gives
Substituting
a=zti-zbj-zti+zbi,
bd?)I
-27
-30 -51 30 -251
-J +(ar+bc+bf)[-;
[&I =
+
(21)
I 1 -3
0
+
(20)
and
the initial strain
-3
of
GI,+z -
’ ,v
{4},+r are functions of (22) and (23). eqn (7) written for time t + 7 with the following - {PM,
,A~$&r; ,
;I,, t+t
4 I’ (24)
LE-CHUNG SHLNJ
484
Therefore, one can obtain {q},+l at time t + 7 and, subsequently, obtain {4},+, and {d},+,. These vectors can be converted to those at time t + At by
0.5 NUMERICAL
RESULTS
Before using the developed program to perform the inelastic dynamic response analysis of a cylindrical structure, it is first used to analyze a static case with a cantilevered beam. An exact solution of the cantilevered beam [ 151 is available for comparison and verification. Figure 2 shows a comparison between the exact and finite element results of a uniform cantilevered beam subject to concentrated lateral and axial loads at its tip. In this figure, P is the axial load at centroid, P, is the yield axial load of the beam, the nondimensional load is the ratio of applied lateral load to the maximum ultimate load (when P = 0) of the beam, and the nondimensional tip deflection is the ratio of actual tip deflection to the tip deflection at the maximum lateral load (when P = 0) for which the beam is entirely elastic. The results obtained from the finite element analysis compare quite well up to the collapse load with the corresponding results calculated from an exact solution by using a continuum analysis. Figure 3 shows the effect of an axial compressive load on the tip deflection of the cantilevered beam. Since the beam stiffness is reduced by the axial compressive load, the nondimensional load at elastic limit is lower for higher axial load and the tip
a
1.4
0 ‘Z E:
1.3
.g a 2
1.2
c
1.1 1.0 0.6
0.7
0.8
nondimensional
0.6
0.7
nondimensional
0.9
1.0
load
Fig. 2. Comparison of finite element and exact results for a cantilevered beam.
0.6
0.9
1
load
Fig. 3. Effect of axial load on the tip deflection of a cantilevered beam.
deflection increases with the axial load. In the analysis, the lateral load was incremented to a value that resulted in the failure of the structure. This failure is indicated in Figs 2 and 3 by the near vertical slope of the load-deflection curves. It should be noted that plastic buckling is not considered as it fails at a much lower load. Inelastic dynamic analysis
To demonstrate the effect of axial thrust on the inelastic dynamic response of cylindrical structures, a self-supported chimney is chosen for the present study as shown in Fig. 4. The chimney is composed of inner and outer shells of revolution; only the outer shell is considered. The elastic limit of the outer shell is set to be 25,220 kN/m2 and the stress-strain relationship of the material is assumed to be bilinear. Because of its curve shape along the height of the chimney, the outer shell is modelled by eight pipe-type taper beam finite elements. The north-south components of the 18 May 1940, El Centro earthquake was used for this inelastic dynamic response analysis; 30 set are considered. A time increment of 0.002 set is used. Figure 5 shows the responses in tip deflection of the outer shell for two loading cases: (1) a purely elastic case and (2) an inelastic case without axial load effect. It can be seen that plastic strains developed in the chimney tend to resist the reverse motion of the chimney, the magnitude of the entire inelastic response curve is smaller than that of the purely elastic response curve. The effect of the external axial compressive load is dependent on the magnitude of initial stresses created in the chimney by the axial load. Figures 6-8 show the effect of the axial compressive load on the inelastic response of the outer shell. The magnitude of initial stresses considered are a percentage of the yielding stress (a,,) namely; lo%, 20%, and 30%. When initial stresses are less than 20% of the yielding stress, the general trends of the response curves for
Inelastic dynamic analysis of cylindrical structures El(m)
O.D.(m) t(OII) 11-74 26.67
369.65
11.89 20,32
12.33 20.32
,
485
------inelaaUc
with 0.2
(I,
!
338.33
,I
307.85
1 ’ I I Fig. 7. Effect of axial load on the responses of tip deflection of outer shell (0.20,).
14,13 31.75
246.89
I I
15.48 40.64
216.41
' I I
Fig. 4. Geometry of outer shell of the chimney.
Fig. 5. Time history responses of tip deflection of outer shell for purely elastic and inelastic cases.
-100-1
these two cases are similar to that for the case without axial load effect and the effects of the axial loads are only on the magnitudes of tip deikction of the outer shell. After 15 set response period, the maximum tip deflection for the case with 0.20, of initial stresses is about 80% of those values obtained for the purely inelastic case. When the initial stresses in the chimney further increases to 30% of the yielding stress, the tip deflection response of the chimney becomes quite different. After the first 15set period, the tip deflection response of the chimney gradually occurs on one side of its original static equilibrium position only (i.e., all negative tip deflections) and the magnitude of the tip deflection is much smaller than those values obtained from purely inelastic analysis. However, the trends of the plastic strains tend to resist the response motion of the chimney are the same for all inelastic cases. This is due to the fact that after the first yielding period the chimney is subjected to unloading and returns to elastic range. Because the permanent strains produced in the preceeding yielding period tend to resist the reverse motion, the subs~uent deflections are less than those obtained from the purely elastic analysis. When the reverse motion continues, the chimney undergoes another yielding period. New plastic strains are produced that result in reduced stiffness. However, the magnitudes of the deflection cannot be greater than those obtained from the purely elastic analysis unless the new plastic strains overcome the preceeding permanent (residual) strains. Thus, although the general trends
Y
Fig. 6. Effect of axial load on the responses of tip deflection of outer shell (0.1~~).
Fig. 8. Effect of axial load on the responses of tip deflection of outer shell (0.30,).
b&HUNG %LW
486
of the inelastic response curves are similar to that of the elastic one the former have a smaller magnitude. CONCLUSIONS A finite element formulation and procedure developed for the inelastic dynamic analysis of slender cylindrical structures using the concepts of subsection and fictitious load to account for the effects of axial load and plastic strains are found to he efficient and satisfactory. The computer program has heen verified by comparing the results of a static case with the analytical solution. The inelastic dynamic response analysis shows that the plastic strains nroduced in the structure tend to resist the response motion. Therefore, the magnitudes of the deflections obtained from the inelastic response solution are less than those obtained from the purely elastic response solution. The effect of the external axial load may have a significant affect on the structural responses. financial support by NationaI Science Council of the Republic of China through Grant No. NSC78-~i~E~ is gratefully acknowledged.
Acknowledgemew-The
3. A. C. Heidebrecht, Vibration of nonuniform simply supported beams. J. Engng Mech. Div., ASCE 93, I-15 (19671. 4. k. US:Rumman, Earthquake forces in reinforced conc&e chimneys. J. Struct. I)&., ASCE 93, X5-70 (1967). 5. L. C. Sbiau and T. Y. Yang, Elastic-plastic seismic response of chimney. J. Strucr. Div. ASCE 106,791-807
(1980). 6. S. Toma and W. F. Chen, Analysis of fabricated tubular columns. f. Struct. Div., ASCE W&2343-2366 f'lQ7Qj. ‘-- F. _‘. Chen and H. Sugimoto, Moment-curvature7. W.
axial-compression-pressure relationship of structural tubes. J. Construct. Steel Res. 5, 247-264 (1985). 8. 1. S. Sohal and W. F. Chen, Local buckling and inelastic cyclic behavior of tubular sections. Thin- Wafied Struct. 6, 63-80 (1988).
9. S: Toma and W. F. Chen, Cyclic analysis of fixed-ended steel beam-columns. J. Struct. Div., AXE 108, 1385-1399 (1982). 10. S. Toma and W. F. Chen, Inelastic cyclic analysis of n&ended tubes. J. Struct. L)iv., ASCE 19S, 2279-2294 i1982). il. H. Armen, Jr, A. Pifko and H. S. Levine, A finite element method for the plastic bending analysis of structures. Report No. RE-347J, Grumman Research Department (1968). 12. K. J. Bathe and E. L. Wilson, Stability and accuracy analysis of direct integration methods. Earthq~e ErzgngStruct. Dyn. 1, 283-291 (1973).
13. E. L. Wilson and R. W. Clouab. Dvnamic resnonse bv step-by-step matrix analysis. Proceedings, Simposium on the Use of Computers in Civil Engineering, Lisbon, Portugal, pp. 45.1%.14 (1962). 14. E. L. Witson, I. Fraboomand and K. J. Bathe, Nonlinear dynamic. analysis of complex structure. Earthquuke Engng Struck Dyn. 1,241-252 (1973). 15. W. Prager and P. G. Hodge, Jr, Theory of Perfectfy Pfastic Sofia& John Wiley (1951). I.
REFERENCES
E. G. Abu-Saba, Vibration of chimney supported on elevated concrete slab. J. Struct. Div. ASCE97,619637 (1971). 2. H. R. Busby, Jr and V. I. Weingarten, Response of nonlinear beam to random excitation. J. Engng Mech. I.
Div., ASCE 99, 55-68 (1973).
.