Material and geometric nonlinear dynamic analysis of steel frames using computer graphics

Compurers

& Structures

Vol. 21. No. 4. pp. 825-840.

0045-7949/85 $3.00 + .OO 0 1985 Pergamon Press Ltd.

1985

Printed in Great Britain.

MATERIAL AND GEOMETRIC NONLINEAR DYNAMIC ANALYSIS OF STEEL FRAMES USING COMPUTER GRAPHICS? SAID I. HILMY and JOHN F. ABEL Department of Structural Engineering and Program of Computer Graphics, Cornell University, Ithaca, NY 14853, U.S.A.

(Received

6 July 1984)

Abstract-A promising analytical model for nonlinear dynamic analysis of steel frames is described. This approach is based on an updated Lagrangian formulation in conjunction with force-space, concentrated plasticity. Kinematic strain hardening behavior is modelled by the bounding-surface concept. The analysis is implemented in a highly interactive, adaptive fashion using computer graphics and a super-minicomputer. Several examples illustrate the effectiveness of the analysis strategy described. L length of the member

NOTATION

mi

vector of yield surface offset cross-sectional area of the member bL}, bb} total offsets of the loading and bounding surface Eli, 82) plastic parameter in the axial and bending direction at end i from the bounding surface d distance measured by the Euclidean norm tdaL1, {dab) incremental offsets of the loading and bounding surface din value of d at the most recent initiation of yield b-G1 vector of incremental displacements IhI, {dq,I vectors of incremental elastic and plastic displacements IdSI incremental forces at the member ends Id&I> Id&I perpendicular and tangent components of dS E Young’s modulus current and bounding values of the modEP, EP’ ulus yield surface equation matrix which includes the gradients of 67 the loading surface gradient of the loading surface at end i {Gi model parameter [HI 2 x 2 matrix which represents the size ratio of the loading and boundary surfaces kll, kr2r .&I. k22 coefficients to be specified by the user according to the experimental data and are used to evaluate the plastic moduli parameters kapi, kfpj axial and flexural plastic moduli at the end i ke, kp, kt elastic, plastic, and tangent moduli lKe1 elastic stiffness matrix [Kg1 geometric stiffness matrix lKh1 diagonal matrix which includes the plastic moduli coefficients [KPI plastic or hardening reduction matrix [Ktl summation of the elastic and geometric stiffness matrices

t This article is an amplification of an extended abstract presented and published at the Eighth World Conference on Earthquake Engineering, San Francisco, California, July 21-28, 1984.

magnitude of plastic deformation

at end

i

{!

Mi,

Mj

MY PL$ Pb

size of the loading and bounding surface in the bending direction bending moments at ends i and j plastic moment size of the loading and bounding surfaces in the axial direction plastic deformation matrix squash axial load vector of element end forces

conjugate force point on the bounding surface U curvature parameter W size of the yield surface

INTRODUCTION

For strong occasional loadings, a structure may experience significant nonlinear behavior. For example, a severe earthquake may induce, in a typical building, repeated excursions far beyond the elastic range. For an economical as well as safe design, such inelastic behavior is indispensable for dissipating a large portion of the energy input by strong loadings. To achieve a realistic insight into the response of structures to strong earthquakes, nonlinear dynamic analysis must be employed. For such a purpose time-history analysis is considered the mainstay for assessing transient structural behavior. This analysis may be improved not only by using models that realistically idealize the true hysteretic behavior of material but also by considering recent advances of computer technology. The main objective of the investigation reported in this article is to establish a nonlinear dynamic computational procedure effective for the design of steel frames. This objective leads to four fundamental requirements: (1) Realistic element model. An element formulation which accurately reflects structural and material behavior is of central importance. Mod-

825

826

S. I.

HILMY

elling which requires subdivision and numerical integration of a member is costly and usually met with reluctance by both researchers and designers. On the other hand, a simple model that uses nonlinear rotational springs at the member ends based on the assumption of elastic-perfectly plastic or Ramberg-Osgood force-deformation relationships may violate or oversimplify the essential characteristics of the material behavior such as the Bauschinger effect. (21 Simple and concise procedure. Besides satisfying correct mechanics basics, the procedure to evaluate the structural tangent stiffness matrix (including both material and geometric nonlinearities) must be simple and concise enough to perform a fast but accurate analysis that covers a wide range of structures from small assemblages to tall buildings. (31 Flexibility and adaptivity. The implementation of the model in a transient analysis system must be flexible and efficient. The number of analysis parameters must be minimized. However, the program should allow the user to interact with the analysis flow to control analysis sequences, redefine parameters, check the reasonableness of the analysis at different time steps, etc. (4) Graphic capability. It is tedious and time consuming to feed the computer with input data and then transform the reams of output listings, at different increments, into readable and understandable plots. In many cases, limited interpretation capabilities alone are enough to discourage the use of nonlinear analysis from the start. A friendly interaction between the user and the computer can be developed through an interactive analysis where current information regarding the analysis and results are displayed graphically in real time. A number of strategies and element models have been implemented and tested. A promising model which satisfies, in the authors’ view, the above four requirements was developed. This model is based on the common assumption of concentrated plasticity, i.e. plastic hinges, are allowed to form only at the member ends. However. a new aspect is that material nonlinearity is accounted for by a bounding-surface kinematic hardening model in force space. The tangent stiffness matrix is formulated in a concise form and implemented in a highly interactive, adaptive program which provides an environment with a considerable degree of flexibility for the performance of nonlinear analysis. The research is conducted in an advanced interactive facility centered about vector-refresh graphics displays and a virtual-memory super-minicomputer. The computer implementation is designed to allow the user to monitor continuously the analysis results, to control the analysis flow, and to intervene at any time. While the analysis is running, a real-time portrayal is provided of yield-surface movement, bending moments, deflections, or other structurai responses selected by the user. For struc-

and J. F.

ABEL

tures of large size, one may perform a batch analysis which stores the results in a readable file that may be graphically replayed at the user’s convenience.

The proposed element model is based on concentrated plasticity assumptions. Here, all the plastic deformations are confined to zero-length plastic zones at the two ends of the member. The part of a member between the two plastic zones remains elastic. This proves to be a reasonable assumption for steel structures subjected to lateral loadings. This concept, which is sometimes referred to as the ‘one-component rotational spring,’ was introduced in simple forms in the mid-1960s by a number of researchers. It was later modified by Nigam[I] to consider the interaction of forces at plastic ends based on the assumption of force-space yield surfaces. Porter and Powell[Zl provided a concise tangent matrix form for Nigam’s model. This provision gained acceptance by many investigators through the 1970s and early 1980s and was implemented and refined by, for example, GrossC31, Argyris[4], and Orbison[S]. To approximate strain hardening effects, Clough[6] introduced the parallel element model in which an elastic element with a modulus selected to represent linear hardening behavior is added to the one-component model. Although this modification can be justified in some cases, especially in the static analysis of steel frames, it cannot model a complex material hysteretic behavior characterized by a Bauschinger effect. lJzgider[7] used essentially the inverse of the Ramberg-Osgood representation to idealize the elasto-plastic force-deformation behavior at the member ends. The beam-column inflection points are assumed to occur at the member midspans. The coupling between axial forces and bending moments is considered only while checking the force interaction curves. Recently, an attempt was made by Chenl8f to improve the one-component model by introducing a discrete plastic hinge concept where each hinge is subdivided into a series of subhinges. The actiondeformation relationships for each subhinge are presented by bilinear functions which, when combined for a series of subhinges, produce a multilinear flexibility for each complete hinge. Finally, this results in a multilinear flexibility relationship for the element itself. However, Chen’s formulation does not appear sufficiently compact. and his analysis has been applied mainly to reinforced concrete or steel pipe structures. In this section, a concise element tangent stiffness formulation for a beam-column with strain hardening is developed. This is followed by a brief description of the kinematic hardening rule adapted to the analysis of framed structures. The extensions

827

Nonlinear dynamic analysis of steel frames using computer graphics

to include three-dimensional to a later article. Element tangent stiffness

structures

are deferred

matrix

The proposed stiffness formulation is in some ways similar to the Porter and Powell elastic-perfectly plastic model[2] but differs in that it includes the effect of strain hardening. This is accomplished by introducing a hardening coefficient matrix and applying a force decomposition principle for plastic flow. The main assumption is that, at each stage of plastic deformation, there is a unique yield surface (y.s.) in the force space (F(S, a) = W) so that strain hardening deformation takes place only for F(S, a) > W. {S}, {a}, and W are, respectively, the vector of element end forces, the vector of y.s. offsets, and the y.s. size. Both F and W may depend on the existing state of forces and on the deformation history. In the current formulation, it is assumed that the size and the shape of the yield surface are kept constant during its movement in the force space (kinematic hardening). For a given state of force {S}, an incremental force {dS} producing plastic flow is imposed. This increment can be decomposed into two components {dS,} and {dS,}, such that Id&} produces no plastic flow and {dS,} is proportional to the gradient of the current y.s. {dFldS}={dSt} = 0,

in which mi is the magnitude of plastic deformation at end i, and {Gi} is the gradient of the loading surface

but

{dFldS}r{dS} > 0.

(1)

Geometrically this means that the vector {dS} is decomposed into components tangent and perpendicular to the y.s.; that is {dS} = {dSt} + {dS,}.

(2)

Within each load or time increment, the relationship between the force increment and the plastic deformation increment is assumed to be linear. Since {dS,} is the only component that produces plastic flow, the following relationship is posed

(6)

{Gi} = [dF(S, a)ldSir dF(S. a)ldSiz “‘IT.

For elements with two hinges [G] contains two columns, one corresponding to each end. Equation (5) becomes

= [GlIm).

(7)

[Glr {d&l = (01.

(8)

But from eqn (1)

Using eqns (2) and (8), one obtains (9)

[Glr ({dS) - {d&J) = {Ol.

The element nodal force increments are related to the elastic displacement increment by (10)

{dSJ = ]Ktl {dqe} in which [Kt] is the summation geometric stiffness matrices

of the elastic and

[Kr] = [Kel + [Kg].

(11)

The elastic stiffness is conventional. However, the particular geometric stiffness matrix used in this investigation is based on the updated Lagrangian formulation and on the assumption of large-displacebehavior, and is given in ment, small-strain Appendix A. Details of the derivation are given by Gattass in [lo]. Equation (10) can be written as (12)

{dSJ = ]Krl {dq) - [Kfl {dq,}. {d&l = ]Khl {dq,1

(3) Substitution

in which [Khl is a diagonal matrix which includes the plastic moduli coefficients to be examined later, and {dq,} is the plastic displacement vector. Equation (3) is parallel to the compliance equation in stress space assumed by Mendelson[9]. In incremental plasticity, the displacement decomposition principle is

of eqns (3) and (12) into eqn (9) yields

]Gl= ([Ktl {dq) - []Kfl +

[Khll {c&J) = (0) (13)

and using eqn (7), one arrives at {ml = [pdml Id41

(14)

and {dqI = {dqJ + Idq,)

(4)

in which {dq,} is the elastic deformation component, and {dq} is the total deformation vector. For an element end i, the associated flow rule is expressed as (5)

[pdml = ([Glr [HI

+ ]Khll]Gl)-’

]Glr1Ktl (144

in which {m} is the vector of plastic deformation magnitudes and [pdml is called the plastic deformation matrix. Finally, using eqns (14), (12), and

828

S. I.

HILMY

(7), one obtains id51 = [[Ktl - [Ktl[Gl

([GIT[[Krl

+ [Khll[Gl)-’

[GIT[~tllkW

(15)

or

and J. F.

ABEL

the plastic range. On the other hand, kp is zero for the elastic-perfectly plastic case. The above concept can be extended to two-dimensional space by considering both axial and flexural plastic stiffnesses. The elastic axial rigidity of the element is

(16) [Kpl is called the ‘plastic or hardening reduction matrix. ’ For the limiting case of no plastic behavior, [Kh] includes infinite coefficients and eqn (15) reduces to eqn (10) which corresponds to the elastic case. Moreover, at the other extreme, elastic-perfectly plastic behavior, [Kp], becomes a null matrix and the solution is essentially similar to Porter and Powell’s[2] or Orbison’s[S]. The term ([G]T[[Ktl + [Kh]][G])-’ in eqns (14) and (15) entails no computational difftculty. As noted in Ref. [5], with a single-equation yield surface, only three non-null forms of [G] are possible: one for a plastic hinge at the ith node, one for a hinge at the jth node, and one for hinges at both ends. In the first two cases, [G] reduces to a column vector, and the matrix requiring inversion in eqns (14) and (15) is simply a scalar. In the case of hinges at both element ends, [G] contains two columns resulting in a 2 x 2 matrix simple to invert numerically. In all cases symbolic manipulation is used to avoid multiplication involving zero values. This leads to a straighforward and concise evaluation of [Kp] if a closed form of [Kh] is provided, which is considered next. Evaluation of [Kh] [Kh] is a diagonal

(17)

in which kapi and kfpi are the axial and flexural plastic moduli at the end i. First, it is helpful to clarify the basic concept used to evaluate these coefficients. Consider a typical axial force-displacement or moment-rotation relationship. From the basic assumption used in eqn (4), the total displacement increment dq can be decomposed into elastic and plastic increments, dq, and dq,. These increments may be expressed as dq = dFlkt, dq, = dFlke, and dqp = dFlkp, in which dF is the force or moment increment and kt, lie, and kp are, respectively, the tangent, the elastic, and the plastic modulus. This leads to llkt

= llke + Ilkp.

The plastic axial rigidity kup may be expressed a multiple of the elastic stiffness kapi = Bli(EAIL)

(18)

At initiation of plastic behavior kp must be infinite to achieve a smooth transition from the elastic into

(19)

as

(20)

in which Bli is a plastic parameter which varies from zero for the elastic-perfectly plastic case and infinity for the elastic case. In contrast to the axial plastic stiffness, the flexural stiffnesses of a beam-column element depend on the curvature distribution along the element which, if all loadings are applied at the joints, is primarily a function of the bending moments at the element ends. From virtual work the following formula is derived for the elastic flexural rigidity kfei

= Ui(EIIL)

(21)

in which Ui = 6/(2 - MjIMi).

(22)

Mi and Mj are the bending

moments at the element ends, and Ui varies from 2.0 to 6.0. The plastic flexural rigidity is related to the elastic one by kfpi

matrix that contains plastic moduli coefftcients for the element’s two hinges. Neglecting the shear effects, one may express [Khl of a two-dimensional beam-column element as [Khl = diagonal [kup, 0 kfpi kup, 0 kfpzl

kae = EAIL.

= B2i[ Ui( EZIL)] .

(23)

Again B2i varies from zero to infinity depending on the hardening value at the element end. As an application of the above equations, a closed-form stiffness matrix for a simple two-dimensional beam element can be easily formulated[l4] and is given in Appendix B. For a beamcolumn element the procedures to evaluate the stiffness matrix follow directly from eqn (15). PLASTICITY

MODEL

Within the general framework of the formulation outlined above there is considerable room for different implementation strategies. The construction of a particular constitutive model requires the specification of F, a, Bl, and B2. Different hardening rules have been investigated and implemented in the analysis. A detailed comparison of these rules is given in [14]. A parallel model similar to the one introduced by Clough[6] provides a bilinear strain hardening behavior. Although the model is simple in concept, it has two major disadvantages. First, compatibility is not satisfied between the two parallel elements.

829

Nonlinear dynamic analysis of steel frames using computer graphics Second, bilinear hardening behavior cannot describe the true hysteretic behavior of material under dynamic loadings. Another concentrated plasticity model that avoids the drawsbacks of the parallel model and utilizes the element stiffness formulations described above was considered. This model idealizes material response in a discrete manner based on a multi-surface Mroz approachll31, but in force space. Although this nested-surface model provides better insight into the hysteric behavior than the parallel model, a large number of surfaces must be employed to obtain the smooth behavior observed experimentally. This is prohibitive for large-scale computation especially if an automatic time step adjustment is used.

of the loading and bounding surfaces, respectively. [HI is a 2 x 2 diagonal matrix that represents the size ratio of the two surfaces and has the form (Fig. 1):

1=

Constant

(24a)

It is possible for the loading surface to contact the bounding surface but not intersect it, and the two surfaces move together while in contact. In such a case, eqn (24) reduces to

and Bounding-surface model

The third alternative considered was selected for detailed development because it provides a smooth hardening model and relatively good computational efficiency. This concept, based on the notion of bounding surfaces, was originally proposed by Dafalias and Popov[ 1l] and subsequently elaborated by these and others. The efftciency of this model is demonstrated in [ 111 by predicting almost perfectly the experimental data of a typical uniaxial stressstrain curve. In the current work this model is extended to model behavior in force space. A comparison of the behavior in force space and stress space[l4] reveals the similarity of both responses for steel structures. In this model the elastic region is represented by the interior of a loading surface F that represents force interaction at the start of yield. This surface is always enclosed by a second surface called the bounding surface which may also move in force space. The initial position of the loading surface need not be at the origin of force space; indeed, an initial offset may be used to approximate the effect of residual stress on initial yielding of the crosssection. The loading surface may translate according to an appropriate hardening assumption. In the current investigation, the direction of the translation is specified, as suggested by Mroz[l3], as the unit vector along the line connecting two conjugate points, one on each of the surfaces. These two points are characterized by identical outward unit normals to their respective surfaces. In analogy to stress space, the translation vector can be described as {da,] =

Uffl - IZll{SI - ([fZl{ad - b,H W’~WkW WdSIT (Uffl - [~ll{~~- [[fflbd - {ab}l (24) in which {daL} is the incremental loading surface translation, and {aL} and {ab} are the current offsets

{daLc) = lH]{dab] +

[[A - [W{.V

(26)

in which {dab} and {daL} are the incremental translation vectors for the bounding and the loading surfaces, respectively. If elastic unloading takes place, the loading surface detaches from the bounding surface and moves in an inward direction, while the bounding surface is assumed stationary. It is possible to allow for gradual change of surface sizes. The bounding surface may also be allowed to move but with a slower rate. However, in the present implementation the bounding surface is assumed fixed except when in contact with the loading surface, and both surfaces are kept at their original size. For a current force point {S} on the loading surface, a conjugate point {S’} on the bounding surface is found by Mroz’s criterion that the points should have parallel outward normals1131. Then the plastic moduli are taken as function of the distance d between {S} and {S’}. The current plastic modulus Ep (in stress space) is suggested by Dafaliasll21 as Ep = Ep’ + h(din) -d-&-d

(27)

in which Ep’ is the bounding value of the modulus, and din is the value of d at the most recent initiation of yield. Both d and din are measured by the Euclidean norm, h is a model parameter which is a function of din and which controls the steepness of plastic variation, and din plays the role of a discrete memory parameter of the most recent event of unloading-reloading. Because of the lack of experimental evidence about behavior in force space, more than one approach is possible in generalizing for multiple forcespace dimensions the concept illustrated above for the uniaxial case. An appropriate choice may be based on logic and a striving for simplicity of interpretation. Any form of generalization, however, must reduce to the results for the uniaxial case. The approach adopted in the current research

830

S.

I. HILMY and J. F. ABEL

follows closely the work of Dafalias on anisotropic hardening for initially orthotropic materials[l 11. First, it is significant that the behavior in forcespace is direction dependent, which means that one can obtain different responses in different directions. There are many reasons for such anisotropy, such as the different geometry of the cross-section in different directions, the different change of the plastic moduli along various axes, and the diverse dimensions of the force vector components themselves. The last factor, for instance, prevents any simplification parallel to the concepts of effective stress and effective strain commonly used in stress space. As a matter of fact, the response of the crosssection in force space is analogous to the behavior of initially orthotropic material[ 111. In the current research, the behavior is considered in normalized space in which the different axes are scaled by their equivalent yield or ultimate strengths. The orthotropic feature can be represented in a systematic way, as part of the general development of the model, by assigning different shape functions to different directions. However, the shape function parameters expressed by d and din are unique as they are measured by the Euclidean norm. The physical interpretation is that during the plastic process, the plastic modulus parameters, although conceptually different, are in genera1 interrelated. This contributes to the nonproportional variation of the force vector components, a behavior which is expected during the plastic process and which is influenced also by the force-interaction equations and the flow rule assumptions used in the formulation of the beam-column strain-hardening stiffness. In the current work eqn (27) has been converted to force space as follows: (28) (29) in which krr , k12, kz,, and k22 are coefficients to be specified by the user according to the available experimental data. The calibration of these coefficients is described in [14] and is not discussed here. The initial distance din averages in a simple way and directly measures the amount of plastic deformation in the opposite loading sense during the most recent load history. For instance, a small value of din implies a small amount of plastic loading in the opposite loading direction. Moreover, the value of din is mainly influenced by the most recent cyclic loading, thus its updating weighs automatically the recent history more than the earlier. The function d/(din - d) is selected to provide the relative steepness as well as the observed basic similarity among the hysteresis loops. As mentioned by Dafalias[ll]. it is the role of

the changing distance d to continuously bridge the gap in the material behavior between an abrupt unloading-reloading process and a process of gradually changing loading direction. A check is added to eqns (28) and (29) so that din is updated to equal d whenever d exceeds the value of din. In such a case, the resulting infinite values of Bi’s are represented by numerically large values, e.g. 106, in the analysis program. A possible restriction of the above equations may arise if one considers extremely large deformation. It is apparent from the uniaxial case that the experimental data, on which the present development is based, has an intermediate range of deformation, which means for example, for a rigid steel connection, a rotation IO to 15 times larger than the maximum elastic rotation. It is noted, however, that after this range of deformation possible failure of the beam-column element may occur due to excessive deformation or local buckling, subjects which are beyond the scope of this study. Therefore, a simple check is provided in the analysis program to set a limit on large deformations based on ductility ratios. Surface

sdection

When a multi-faceted loading surface is used, additional computations are necessary to identify the active facet and to prevent the forces at the nonlinear element ends from attempting to cross over from one facet to an adjacent one. Further, multi-faceted surfaces pose additional problems in the determination of the correct gradient matrix [G] at facet intersections. Hence, an admissible singleequation surface with full slope continuity is used in the current approach to model both the loading and the bounding surfaces. The same equation (with different sizes) is used for both surfaces to allow contact or sliding but not overlapping of the surfaces, as suggested by the Mroz mapping rule which is used in this work. The surface equation was developed by Orbison et al.[ 151 to model plastic hinge formation of a WI 2 x 31 steel cross-section, as most representative of the behavior of light- to medium-weight American shapes. The method used to develop the surface is a combination of trial and curve Iitting for the combined actions needed to cause complete plastification of the cross-section. The yield surface, in the two-dimensional space. has the form 1.15~’

+ m7 + 3.65p’m’

=

W

(30)

in which p = (P - a,)/Py is the ratio of the axial force to the squash load, m = (M - az)lMy is the ratio of the strong axis bending moment to the corresponding plastic moment, and ~1~and a, are the current offsets of the surface in the P and M directions, respectively.

Nonlinear dynamic analysis of steel frames using computer graphics Figure 1 shows the application of eqn (30) to model the loading and the bounding surfaces by scaling it to two different sizes. It is shown that, although the equation is based on the assumption of complete plastification of the cross-section, it gives an acceptable approximation to the initial yielding surface. Figure 2 shows the movement of the loading surface inside the bounding surface once the plastic process starts. According to this movement the values of the plastic parameters Bi’s can be obtained based on the initial and current distances din and d and the corresponding shape functions. An iterative procedure is typically required to control the force drift from the loading surface. This drift is caused by the assumption of a constant value of the surface gradient during a step increment. A bisection algorithm is used to keep the force on the surface within a tolerance specified by the user. The force return path is perpendicular to the surface. ANALYSIS IMPLEMENTATION

The above model has been implemented by the first author for two-dimensional dynamic nonlinear analysis program by building on the initial developments of Gattass[lO, 161. This capability includes static, dynamic, modal, and buckling analysis with

831

real-time graphical feedback. Updated Lagrangian procedures, a skyline solver, and dynamic memory allocation are employed. Symbolic manipulation is used to perform all the required matrix formulations. In the next section, a brief review of the strategies used to implement the proposed model is given. This is followed by a description of the interactive capabilities. Global solution methods

Direct integration implicit analysis is used to solve the coupled equations of motion for the system described above. The transient analysis can be performed by either the two-parameter Newmark Beta or the Wilson Theta method. These methods are well known and need not be reviewed here. However, a number of modifications have been made to the basic method to account for material nonlinearity and other features of the problems being considered. In both cases the updated Lagranian procedure is employed, where the reference configuration is updated at the end of each increment. The analysis follows a Newton-Raphson incremental-iterative procedure based on an energy tolerance convergence criterion. The time increment size is controlled during the analysis so a change of the element hardening (onset of yielding or contact with the bounding sur-

Fig. 1. Use of the same single-equation surface to idealize both the loading and the bounding surface.

832

S. I. HILMY and J. F. ABEL I

P/PY

Fig. 2. Movement of the loading surface during the plastic process.

face) does not occur within the increment itself. Along this line, at every new configuration the analysis starts with an attempted time-step size. The resulting displacement, velocity, acceleration, and element nodal force increments are obtained and saved. The total element forces at the end of the increment are checked against the loading surface in its current position. This check is performed for each element end. (Various provisions are made to accelerate this procedure by skipping unnecessary comparisons). If it is found that yielding would begin, or an overlapping of the loading and the bounding surfaces wouId occur, the revised load increment to just produce plastification or surface contact (within a tolerance provided by the user) is calculated by a bisection method. A new value for the distance between the force vector point which lies on the loading surface and the conjugate point on the bounding surface is calculated and temporarily stored, A linear interpolation is used to find the current values of the displacement, velocity, and acceleration vectors. For this purpose, circular vectors are assigned for the storage of the intermediate steps. The analysis proceeds by updating the values of the distance {d} and the offsets of the yield surfaces {ai} at every yielded member end. As the member stiffnesses are being formed, the plastic deforma-

tion matrix [pdm] of eqn (14a) is calculated and saved. After the hinge displacements have been determined, [pdm] is used to compute the plastic deformation magnitude {m}. A negative value of {m} indicates immediate unloading. In such a case, the force point is no longer to be constrained to remain on the yield surface; therefore, the result of the step is ignored, the member and structural stiffnesses are reformed, and the step is repeated. The process must be continued until all plastic deformation contributions are positive. Since for steel cross-sections the values of kij of eqns (28) and (29) are selected to provide steep variations of the plastic moduli at yield initiation, the load or time step size must be modified to reflect this characteristic behavior, i.e. to obtain a smooth and gradual transition from the elastic to the plastic range. A step control ~gorithm is implemented to adjust automatically the step size according to the slope of eqn (29). However, during unloading the step size is changed automatically to its initial, before-yielding value. This algorithm noticeably smooths computed load-deflection curves for subassemblages or structures with a small number of plastic elements. However, for a structure with a large number of plastic elements it is found that the use of such an algorithm results in little difference of the overall load-deflection although it increases

833

Nonlinear dynamic analysis of steel frames using computer graphics

the total number time.

of incremental

Adaptive interactive

analysis

steps and CPU

The analysis is incorporated in a comprehensive interactive graphic analysis and design capability which provides an environment with a high level of user flexibility for the performance of nonlinear analysis[l6]. In this environment, analysis may be performed either in a time-sharing mode or assigned to a batch queue for execution at a later time. The interactive graphics techniques to control the flow of the analysis and to monitor results are based on a hardware configuration that consists of a vector refresh display terminal, a digitizing tablet, an editing terminal, and a super-minicomputer. The graphics displays consist of various menu pages such as the one shown in Fig. 6. The structural response is monitored in 3 different viewports, and the user may switch information from one to another. The user interacts with the computer by pointing to the commands on the menu either to move to another menu or to perform a required procedure . The concepts of parallel and sequential analyses are introduced to facilitate selection of analysis parameters and understanding the nonlinear structural characteristics. Free vibration, buckling, and deformation mode analyses involving the tangent matrices of the above formulation are examples of possible parallel analysis. Database organization and graphical feedback are bases for analysis sequences. The user has the option to control the flow of the analysis and intervene when problems appear. He can also redefine parameters, masses, element properties, load histories, and earthquakes which are stored in a disk file. The options to display bending moments, displacing shapes, yield-

surfaces movement, or plastic hinge information can be selected while the analysis is progressing. Special parameter pages permit user specification and checking of the experimental data of eqns (28) and (29) and all the other parameters needed for performing nonlinear analysis. The user can interactively change surface information and assign residual stresses to the cross-sections by changing the initial offset of the loading surface. For the batch analysis mode, provision is made to save and display the result-output files with full advantages of the graphics feedback originally developed for real-time analysis. More details of the implementation are given in [ 141. EXAMPLES

In this section, four abbreviated examples are selected to indicate the effectiveness of the computational approach developed from the above ideas. Details of these and other examples are given in [14]. For the examples described here, the bounding-surface approach was used. Unless otherwise stated, the following hardening parameters were selected: k,, = 0.03, klz = 0.5,kZ, = 0.03, and k22 = 2.0. These parameters are only estimated and are not obtained by rigorous calibration process. Figure 3 shows the variations of plastic moduli along the distance between the loading and the bounding surfaces. Geometric nonlinearity is considered for all incremental analyses. Zero values are assigned to residual stresses and damping parameters. The CPU times given are for a VAX 1l/780. Beam-column

connection

Figure 4 shows an assemblage consisting of a beam-column intersection. The assemblage has 4 elements each of 125 in. length and W12 x 35 cross-

4 eo

T

Bi

B,=kilt1

+ klpd/(din

- d)

1

\

k. v-_-_-Ly__ 0.0

* dn

Fig. 3. Plastic modulus parameters Bi’s vs (din -

(din-d)

d), eqns (28) and (2%.

834

S. I. HILMYand J. F.

ABEL

Deflected (Magnif.

PLastic

Factor

-

10)

Structure

Hinge

Fig. 4. Assemblage

Shape

of beams and columns to model intersection

-0.04

Fig. 5. Moment-rotation

at an intermediate joint.

-0.02

relationship for the central joint of the assemblage of Fig. 4.

Nonlinear dynamic analysis of steel frames using computer graphics

PROGRAM

835

OF COMPUTER

STRUCTURAL

ENG.

CORNELL

GRAPHICS

DEPARTMENT

UNIVERSITY

FORCES INTERACTIONAT NODE 4

6

vi-

READ BATCH WRITE FILE -

ROTX

-

ROT

+2

ROTY

t -

ZOOM

PAN HARD

.I

“73o.d

-43606

t

RESET COPY

FULL

VIEW

’ BENDING

EXIT

MOMENT

Fig. 6. Analysis information for single-bay portal frame (Example 2).

Time

Fig. 7. Time history for the horizontal displacement of the top right node of the portal frame (Example 2).

t

836

S.

I.

HILMY and J. F. ABEL

section with a yield stress of 36.0 ksi. A lumped mass equal to 0.356 k-sec’/in. is assigned to the joint which is subjected to a concentrated moment that follows a sinusoidal load history of an amplitude of 13,860 k-in. at a period of 5 seconds. 200 time increments were used with an average time step equal to 0.0474. The complete real-time solution consumed 84 seconds in CPU time; this includes execution time as well as graphic manipulation to update display monitors and show yield-surface movements. Figure 4 shows also the magnified deflected shape of the assemblage at time = 0.89 sec. Hinges are indicated by small circles adjacent to the plastic end. The hysteretic cyclic behavior of the connection rotation is depicted in Fig. 5.

w

Single-bay portal frame

A single-bay, one-story portal frame is subjected to vertical and horizontal concentrated loads as shown in the upper left monitor of Fig. 6. These loads follow a cyclic sinusoidal load history with an amplitude of 4.0 kips and a period of 5 seconds. A light-weight steel cross-section of W14 x 22 was used for both the beam and columns. The beam and left column have a yield stress of 36 ksi. while the right column yields at 10 ksi. The fundamental mode of deformation has a period of 3.66 seconds. The analysis was performed in 200 time steps with an average time of 0.081 second and a total CPU time of 90 seconds. Figure 6 also shows the forcepoint trace of the lower end of the right column at

W 6x9

6x9

7

A 16

w 10x12 W 12x14

W 12x14

w 12x22 3 @ 115”

w 12x19

w 1?x19

500 ” I W 16x22

1

W 14x22

w 14x22

143”

w 16x36

Fig. 8. Four-story frame.

8 Y

1.0

2.0

4.0

3.0

6.0

5.0

Time

Fig. 9. Moment history at the lower end of member

16 of the four-story

frame.

Nonlinear dynamic analysis of steel frames using computer graphics .

+

WlEX36

M4x26

Wi4x26

I

360’ B 1

837

I

266’ ’ 1

200”

Fig. 10. Ten-story structure.

(4

(4

(cl

Fig. 11. Behavior of the ten-story structure subjected to vertical loading and El-Centro earthquake. (a) Deflected shape at step 52. (b) Bending moment diagram at step 159. (c) Deflected shape at step 199.

the end of the analysis. The upper monitors portray the bending moment diagram and the relationship of the applied force vs. the horizontal displacement of the upper right node. The time history of this displacement is shown in Fig. 7 for both linear and nonlinear analysis. Four-story building

A four-story two-bay building is shown in Fig. 8. A concentrated mass of 0.1 k-sec21in. is assigned

to each joint. All the members have a yield stress of 36 ksi. The fundamental mode of vibration of the structure has a period of 1.91 sec. The structure is subjected to the first five seconds of the N-S component of the 1940 El Centro earthquake. The moment history of the lower end of the fourth story central column is shown in Fig. 9 for both elastic and nonlinear cases. The CPU time for the linear analysis was 35 seconds for 62 steps of 0.1 seconds each, while the nonlinear analysis was performed

838

S. I. HILMYand J. F. ABEL

-c

Plastic

Analysis

With

Plastic

Analysis

Without

Linear

Vertical

Loads

Vertical

Loads

Analysis

-I1111111111111111114

______-_

+

2.50

1.25

0.

'TIME ____--.----

~----

Fig. 12. History of the horizontal displacement

in 102 seconds with 112 steps. Plastification observed in eleven different locations.

was

Ten-story frame

A ten-story, three-bay frame consisting of 70 elements and 44 nodes is shown in Fig. 10. The crosssections of the members are similar to those proposed in [ 171 to achieve weak-beam, strong column behavior. A concentrated mass of 0.2 kip-sec’lin. occurs at each node. All the members have a yield stress of 36 ksi. The fundamental mode of vibration of the structure has a period of 2.79 sec. The building is subjected to vertical nodal loads of 3.0 kips at the inner nodes and 2.0 kips at the outer nodes. These loads are constant during the dynamic excitation. The structure is subjected to the N-S component of the El-Centro earthquake. Three different analyses are performed. These include linear analysis, inelastic analysis with the vertical loads applied to the structure, and inelastic analysis without vertical loading. The bounding-surface model is

3.75

5.00 _--

of the upper-right node of the ten-story frame.

used in the two inelastic analyses. The model parameters are: kll = 0.03, k,;? = 0.50, kz, = 0.05, and kz2 = 3.0. The sizes of the loading and bounding surfaces are 0.89 and 2.2, respectively. 200 steps are employed in each case. The nonlinear analysis is sent to a batch queue and is performed in 10 minutes of CPU time. Figure 11 shows a few glimpes of the “playback” of the analysis. The plastic hinges start to form at the lower levels and at the supports of the structure. Next, these plastic hinges are released, and new ones form at upper levels. This slow migration of the plastic hinges is shown in Figs. 1l(a)(c). The structure exhibits excellent ductile behavior with a uniform distribution of the bending moments as shown in Fig. 11(b). Figures 12 and 13 show a comparison of the behavior of the three different analyses. In the first 3.7 seconds of the earthquake record, the vertical loading has only a slight effect on the variation of the horizontal sway of the structure.

839

Nonlinear dynamic analysis of steel frames using computer graphics

0.2 ai

Nonlinear

-0.3

4 I

I

.&

Nonlinear

8

Elastic

Analysis

With

Without

Vertical

oa

s

Analysis v

I

0.

I

I

1.00

I

I

I

I

I

I

I

IIIIIIIIIIII+

3.00

2.00

4.00

5.00

‘TIME Fig. 13. History of the vertical displacement

of the upper-right node of the ten-story frame.

CONCLUSION

in this paper incorporates, in the authors’ view, a promising tool to provide a clear insight into nonlinear dynamic behavior of steel frames. The concise formulation of the tangent stiffness matrix along with the bounding-surface approach which employs realistic and admissible single-equation yield surfaces, provides a realistic, effective, and reasonably fast nonlinear analysis. The implementation in an interactive computer graphics environment has an obvious advantage in that it permits the user to trace and understand structural nonlinear response. The approach

J

described

Acknowledgements-The research presented herein was made possible by the Grants No. PRF-7815357 and No. CEE-8117028 from the U.S. National Science Foundation. The two-dimensional program used in the analysis was initiated by Dr. M. Gattass. The collaboration of Dr. C. Pesquera during the early phase of the research is greatly appreciated. The authors also gratefully acknowledge the contributions from Prof. W. McGuire, as leader of the pro-

ject, and Prof. D. Greenberg, as Director of the Program of Computer Graphics. REFERENCES

1. N. C. Nigam, Yielding in framed structures under dynamic loads. J. Engng Mech. Div. ASCE %(EM5), 687-709 (1970). F. L. Porter and G. H. Powell, Static and dynamic analysis of inelastic frame structures, Report No. EERC 71-3, Earthquake Engineering Research Center, University of California, Berkeley, California (1971). J. L. Gross, Design for the prevention of progressive collapse using interactive computer graphics, Dissertation submitted in partial satisfaction of the requirements for the degree of Doctor of Philosophy, Cornell University, Ithaca, New York (January 1980). J. H. Argyris, B. Boni, U. Hindenlang and M. Kleiber, Finite element analysis of two and three dimensional elasto-plastic frames-the natural approach. Comput. Methods Appl. Mech. Engng 35, 221-248 (1982). J. G. Orbison, Nonlinear static analysis of three-dimensional steel frames, Dissertation submitted in partial fulfillment for the degree of Doctor of Philosophy, Cornell University, Ithaca, New York (May 1982).

840

S. I. HILMYand J. F. ABEL

6. R. W. Clough, K. L. Benuska and E. L. Wilson, Inelastic earthquake response of tall buildings, Proc. Third World Conf. on Earthquake

7. 8.

9. 10.

Engineering,

New

Zealand, Vol. II, Session II (1965). E. U. Uzgider, Inelastic response of space frames to dynamic loads. Comput. Structures 11,97-l 12 (1980). F. P. Chen, Generalized plastic hinge concepts for 3D beam-column elements, Dissertation submitted in partial fulfillment for the degree of Doctor of Philosophy, University of California, Berkeley (November 1981). A. Mendelson. Plasticity: Theory and Application. Macmillan, New York (1968). M. Gattass, Large displacement. interactive-adaptive dynamic analysis of frames, Dissertation submitted in partial fulfillment for the degree of Doctor of Philosphy. Cornell University, Ithaca, New York (May 1982).

P

-v (12a + 1.2)P

M - VL (6a + O.l)PL (4a + 2/15)PL’

Sym.

ophy, Cornell University, Ithaca, New York (May 1984). 15. J. G. Orbison, W. McGuire and J. F. Abel, Yield surface applications in nonlinear steel frame analysis. Comput. (1982).

I;EA L

0

0

12EI L’

i,,

displacement ics. Comput. 17. V. V. Bertero of multistory

Mech.

Engng

APPENDIX A Geometric

stiffness for a two-dimensional

-P V -M + VL P

V (12a + 1.2)P -(6a + O.I)PL -v (12a + 1.2)P

sectional area. APPENDIX B Tangent stiffness matrix for a two-dimensional element including material nonlinearity only?

--

EA L

6EI -L*

36(EI)’

5 - L3

)(

24(El)’

-5,,

8(E1)2 5, 0

4EI

(- L

Tucson, Arizona, U.S.A. 13. Z. Mroz, On the description of anisotropic workhardening. J. Mech. Phys. Solids 15, 163-175 (1967). 14. S. I. Hilmy, Adaptive nonlinear dynamic analysis of three-dimensional steel framed structures with interactive computer graphics, Dissertation submitted in partial fulfillment for the degree of Doctor of Philos-

beam

0

0

Sym.

Proc. Int. Conf. on Constitutive Materials, January lo-14,1983,

1

in which a = Z/AL’, P is the axial force, V the shear force, M the moment, L the member length, and A the cross-

(- -

Lawsfor Engineering

element

-M (6a + O.l)PL (2a - lI3O)PL’ M -(6u + O.l)PL (4a + 2/15)PL?

L

versus sophistication,

557-513

large analysis with real time computer graphStructures 16, 141-152 (1983). and H. Kamil, Nonlinear seismic design frames. Can. J. Civ. Engng 2, 494-516

4EI

[Kl =

33,

(1975).

- 12EI 7+

36(EZ)’

Appl.

16. M. Gattass and J. F. Abel, Interactive-adaptive,

11. Y. F. Dafalias, On cyclic and anisotropic plasticity: (i) a general model including material behavior under stress reversals, (ii) anisotropic hardening for initially orthotropic materials, Dissertation submitted in partial fulfillment for the degree of Doctor of Philosophy, Dept. of Civil Engineering, University of California, Berkeley (1975). 12. Y. F. Dafalias. Modeling cyclic plasticity: simplicity

r

Meth.

1 in which 5 = 4EI + Km Km=

i,

4( EZ)2 >

.

L

B2UiF

t In this special case, the plastitication is assumed to be caused by the bending moment at end i only and the coupling between axial and bending actions is ignored.