PREDICTION OF POST-ELASTIC SEISMIC RESPONSE OF STRUCTURES BY A MODE SUPERPOSTION TECHNIQUE

SYNOPSIS: For the determination of the dynamic response of a structural system to earthquake ground motion the mode superposition technique offers an alternative approach to the well established direct integration method. Whereas the principle of modal superposition is encountered commonly in elastic analyses, the response of a yielding structure ispredicted almost universally by direct integration of the equation of motion. In this paper the principles of modal superposition are extended into the post-elastic domain. Among the advantages which accrue from the modal based transformation are a potential reduction in the number of dynamic degrees-of-freedom considered in the solution, and a new insight into the structure response provided by the instantaneous dynamic properties of the structure.


Nakao and Takano (14) described the analysis by mode superposition of elastoplastic frame structures.
The examples considered again were very simplistichowever, they concluded the method was effective for analysing the elasto-plastic seismic response of framed structures.
Bathe (15) recognized the generality of the mode superposition formulation and suggested the approach when the response may be modelled by relatively few coordinates.
Little research effort has been directed towards the application of modal superposition techniques to the dynamic time-history response determination of EARTHQUAKE ENGINEERING, VOL. 16, NO. 3, SEPTEMBER 1983 typical multi-storey frame structures subjected to earthquake ground motion when the member deformations are extended beyond the elastic limit.
In this paper, the solution of the non-linear structural vibration problem by mode superpositon is outlined.
In the theory presented the mass matrix, M, is assumed to be constant.
A lumped mass rather than a consistent mass formulation is adopted.
The concepts developed apply to a displacement dependent mass matrix, however the lumped mass representation enables some economy in the solution of the set of dynamic equilibrium equations. The stiffness matrix, K, is assumed to be displacement dependent, however geometric non-linearity effects have been neglected. Since the eigenproblem is a function of. the stiffness matrix, as well as the mass matrix, a change in the stiffness properties implies a change also in the corresponding eigenproperties.
In the step-by-step response prediction, during periods of material yielding, the eigenproblem may require solution many times.
An efficient algorithm for the eigenproblem solution is fundamental to the successful implementations ofthe proposed solution strategy.

There is one fundamental difference between the solution of the dynamic equilibrium equations by direct integration and by mode super-position.
In direct integration, the equations of dynamic equilibrium are integrated directly in the original finite element basis whereas before solving the equilibrium equations, the mode superposition method requires a transformation to a new basis.
A large proportion of the computational effort in the mode superposition approach is directed towards the derivation of a set of vectors which together form a suitable transformation basis.
Generally, the eigenvectors are the most suitable set of vectors, however, as discussed in detail subsequently, an alternative set of vectors may sometimes be computationally more efficient.
premultiplying by T , the new form of the equilibrium statement is:

T MTy T CTy T KTy = T F(t) (3)
The displacements derived in the new system may be related to those in the original system by means of equation (2).

So far no comment has been made on the form of the matrix (T).
There is a limitless choice of coordinate transformations which could be applied as in equation (2).
The effectiveness of the mode superposition "approach lies in the choice of the column vectors for the transformation matrix, T. The eigenvectors corresponding to the n-free vibration modes of the above system provide an eminently suitable basis.
The matrix T is constructed as a square matrix of order n with the columns comprising the n individual eigenvectors.
A typical element t..

represents the .th" component of the jtn eigenvector^"
For linear elastic studies, rewriting the transformation matrix in more familiar form as: The equation of dynamic equilibrium in the finite element coordinate system may be written in terms of a structural mass matrix M, a damping matrix C, and a stiffness matrix K as The displacements, x, and related velocities, x, and accelerations, x, may be expressed in terms of a new coordinate system by means of a transformation of the form: If the viscous damping matrix representation, C, is chosen such that the same transformation as uncouples the mass and stiffness matrices also uncouples the damping matrix, the n equations of motion become fully uncoupled and their solutions may be calculated independently.
Damping which satisfies the above criteria is commonly referred to as orthogonal, classical, or proportional damping.
The final form of the equilibrium expression is then set of n equations represented by: represents the damping ratio for mode 1 is the natural frequency of mode 1 Following the transformation to an uncoupled system, the displacement response may be computed by a variety of techniques including the direct integration approach. One immediate benefit resulting from the uncoupled system is the matrix manipulation steps involved in the solution of each of the n equations are reduced to trivial calculations.
The selection of an appropriate integrations time step length is inherent in the development of the direct integration method.
A time step formula, based on the highest participating frequency, (6) is commonly adopted: 0.01 < =^ <_ 0.1 (9) where t is the integration time-step length (seconds) T is the period of the highest frequency p mode assumed participating (seconds).
The cost of solution using direct integration may be reduced by using a varying time step related to the mode frequency.
A coarse step may be adequate for the first few modes (or low frequency modes) and the step-length further subdivided for the higher modes.
The overall response is then approximated by superimposing the solutions when the time intervals coincide.
The appropriate length of time step is now related to the accuracy with which it is desired to track the selected earthquake ground motion record in each coordinate.
When transformed to generalized form as in equation (8) the equations may be solved using an integration operator based on the Duhamel integral.
The advantage of the Duhamel-based operator is that it is precise -it has no artificial damping or phase velocity error.
A precise operator based on Duhamel 1 s integral has been presented by Dunham, Nickell and Stickler (2).
A particularly effective recursive relationship, based on the approach due to Cronin (16), has been presented by Dempsey and Irvine (17) .

PROBLEM CONDENSATION:
Benefits which accrue in rewriting the equations of dynamic equilibrium in terms of generalized coordinates result in the first instance from the simplification of the matrix operations in the solution schemes and in the second from the increase in scope in selection of the particular numerical integration scheme. If these were the only advantages, mode superposition schemes in general would . be computationally inferior to the efficient direct integration schemes when applied to non-linear analysis.
When the contribution of each uncoupled system to the gross response of the structure is examined, a major potential of the modal transformation approach -a reduction in the number-of-degrees participating In the response -is identified.
For the frame structure typical of a muIti-storey building, the response is almost fully described by the first few lateral displacement modes.
For the structures considered in this study the first p modes are assumed significant in contributing to the response. The transformation matrix,. (¥) , is a matrix of order n x p, the columns comprising the first p eigenvectors. The generalized equation of dynamic equilibrium as in equation (7) may be rewritten: The displacement response is computed by solving p equations rather than n as in the complete solution.
The efficiency of the reduced scheme is a function of the ratio of p to n -the greater the difference the more attractive will be the mode superposition approach.
In the discussion thus far it is assumed that the transformation matrix, (T) , is constructed from the first p eigenvectors of the mathematical model of the structure.
The transformation matrix is then written as ($).
The only mathematical restriction placed on the selection of the transformation matrix is that the order of the matrix must be compatible with that of the mass, damping and stiffness matrices to ensure the triple matrix products are valid.
A poor choice results in an inadequate representation of the distribution of inertia loads. A large number of vectors would be necessary to achieve an accurate preduction of the dynamic response of the structure.

Whilst it is desirable to use the eigenvectors as a basis for the transformations of the equations of motion, they do not always provide the most efficient basis.
If a set of vectors are available which, although not the exact eigenvectors for the structural system, nevertheless are a close approximation to the true eigenvectors, then these vectors may provide a satisfactory basis for the transformation matrix.
It is logical that any set of displaced shapes which approximate the correct mode shapes will provide a reasonable prediction of the distribution of the inertia loads. The set of vectors so derived retain one property of the eigenvector basis as mentioned previously.
The inertia loading may be adequately represented by a small number of vectors and thus the order of the generalized, set of equations is kept to a minimum.
However the vectors generally do not satisfy the orthonornality relationships of equations (6a) and (6b).
Since the system of transformed equations are no longer diagonalized by the selected basis, the system of equations are coupled.
The numerical effort required to solve the transformed set of equations is increased when compared with the uncoupled system. The choice of integration scheme suitable for computing the displacement response is restricted since some schemes are specifically developed for an uncoupled set of equations.
The response may still be determined efficiently since the order of the generalized set of equations, p, is less than n, the order of the structural system as initially described.
The situation where an approximate set of eigenvectors is a more efficient transformation basis than the exact eigenvectors is outlined in the following section.

NON-LINEAR SYSTEMS:
The principles presented above may be applied equally to the solution of non-linear problems.
The solution may be considered as equivalent to the solution of a series of linear problems. Two alternative formulations are discussed.

The equation of dynamic equilibrium is written in incremental form, to give:
Ay + C* Ay + K* Ay = T AF(t) MAx whe re C fc = damping matrix at time t and tangent stiffness matrix at time t

The lumped mass approach is assumed in the formation of the mass matrix (M), hence (M) is considered time invariant.
The modal transformation to be outlined applied equally to the instance where the mass matrix is also time dependent. Such a case may arise where the "consistent" mass representation is adopted.
The first solution approach is to use the same basis for the transformation matrix as is used in linear elastic analysis.
The columns in matrix (T) comprise the first p eigenvectors of the undamped elastic vibration problem: The p resulting equations no longer correspond to uncoupled generalized coordinate representations.
A high level of coupling may be displayed in the matrix representation -the generalized matrices K* and C* are full rather than diagonal as in the elastic case.
The solution of the coupled set of equations is not intractable and may be efficiently computed by direct integration procedures. Although an uncoupled set of equations is not derived, the order of the system is reduced from n to p, enabling a reduction in computer effort in solving for the displacements.
The reason for considering a mathematically approximate eigenvector w basis is that the eigenproblem, which represents a high proportion of the computer demand, is solved once only. Where considerable local yielding is incurred the elastic configuration may provide a very poor representation of the distribution of inertia forces.

XtM
The repeated solution of the eigenproblem considerably increases the computer demand. Some economy should be possible from knowledge of the eigenproperties immediately prior to the current stiffness state.
In this regard, the subspace iteration algorithm (10) is a suitable choice of algorithm since the previously calculated eigenvectors may be used as a first approximation to the solution of the current eigenproblem. For small variations in stiffness, the eigensolution should converge in few cycles.

J (12)
The basis will be acceptable provided the nqn-linearities are confined to a local section of the structure such that the displaced configuration, allowing for yielding of members, resembles the elastic displaced shape; on substituting

T T MT Ay + T T C t T Ay + T^K^TAy
The transformed set of equations represented by expression (17) may be coupled or uncoupled depending on the form of the damping matrix C*.
Viscous damping effects are modelled by a Rayleigh-type representation (18, 19) which relies on three parameters for its description.
K t = current tangent stiffness matrix K Q = original elastic stiffness matrix If 3 q is specified equal to zero, the instantaneous eigenvectors will be orthogonal to C fc and thus: A series of p uncoupled equations remain for solution.
During yield excursions, (K^) may be very low in comparison to (K ). The damping represented by ° 3 Q K q type is that it retains a high value throughout yielding phases.
Damping of the former type may lead to numerical instabilities in the solution phase.
The damping matrix C* is not generally of uncoupled form when represented by equation (18). An uncoupled set of equations is desirable to simplify the solution procedure. With a low level of coupling it may be acceptable to neglect the off-diagonal terms in C*.
An estimate of the error associated with diagonalization has been presented by Warburton and Soni (20). From a study of several simple multi degree-of-freedom systems subjected to harmonic, random, or aperiodic excitation a limiting condition for the modal damping ratio of the dominant mode or modes (more than one mode must be considered for structures with closely spaced antural frequencies) was proposed in terms of natural frequencies and the elements c.. of (C*) . For the r th mode the  = 1, 2, 3 , . . . , r-1, r+1, .. . . n. For the maximum errors in the major response quantities to be of the order of 10 per cent e was chosen as 0.05.
The alternative approach is to use a direct integration procedure and include the coupled matrix (C*).
The analyst is required to select the number of modal coordinates which will adequately represent the coupled, damped response. Studies indicate that this number is similar to the number required in a standard mode superposition analysis where no modal coupling exists (21).
In the preceding discussion two alternative approaches for the solution of the dynamic equilibrium equations by mode superposition are presented.
The. first uses a transformation basis derived from the initial undamped eigenvectors of the structural system.
The second uses a time varying set of eigenvectors derived from the current stiffness properties of the system allowing for material yielding.
This second approach requires the eigenproblem to be resolved each time the stiffness characteristics of one, or more, elements are altered.
An engineering approach suggests a third alternative could offer the most computationally efficient solution.
As has been stated previously an approximate set of eigenvectors can provide a mathematically acceptable transformation basis.
The third approach establishes a transformation matrix using the initial undamped elastic structure configuration to derive the eigenvector basis. The step-by-step integration of the dynamic equilibrium equations is initiated incorporating this basis.
During the integration, depending on the component element yield strength compared with the intensity of the earthquake ground motion, at the end of a particular time increment some element forces may reach their yield limit.
The stiffness matrix is modified to account for the change in element properties.
Providing that yield is restricted to a local region the influence of the local yielding may be considered insignificant when compared with the overall stiffness characteristics of the system.
In this case the set of eigenvectors derived from the initial structure configuration still provides an accurate description of the distribution of inertia forces in the structure.

The generalized stiffness matrix (T^K t T Q ) is a full matrix of order, p, because the eigenvector basis is not orthogonal to (K t ) when t f o
The generalized damping matrix (C tQ ) is a full matrix also.

When the system of generalized equations exhibits strong coupling, signified by large off-diagonal terms in the generalized stiffness matrix, then the eigenproblem is resolved and a new basis is generated.
The generalized equation of equilibrium is then written as:

(K t ) is the tangent stiffness matrix at time t A byproduct of the scheme whereby the solution of the instantaneous eigenproblem is considered is that the analyst now has access to the instantaneous frequencies and mode shapes.
The instantaneous periods are of interest because they represent the dynamic effects resulting from fluctuations in the tangent stiffness during periods of strong shaking.

Studies (22) have indicated there does not seem to be any simple criterion by which the earthquake damage potential may be correlated to a particular structure configuration.
These studies have been restricted to displacement related comparisons associated with a range of earthquake ground motion records. It would be advantageous to attempt a correlation on a dynamic basis.
Utilizing the above calculated values, the structure instantaneous frequency spectrum can be compared with the earthquake ground motion spectra.
Resulting from these studies, the sensitivity of a structure to a particular earthquake motion may be clarified.
The solution strategies discussed in this paper have been incorporated into a general three-dimensional frame analysis computer code.
In Table 1 the program structure is demonstrated in block form. For a more complete description the reader is directed to the source document(30).

EXAMPLE PROBLEM:
To demonstrate the implementation of the mode superposition scheme, and to enable a comparison to be made with the results preducted by direct integration, the analysis of a reinforced concrete frame structure is now presented.
The test structure was derived from the proportions and strength levels of the six-storey plane frame analysed by Jury (29).
The frame has been designed to comply with the earthquake lateral loading provisions of the current New Zealand Loadings Code (28).
The frame geometry is shown in Figure 1, and other relevant parameters used in the dynamic analysis are presented in Table 2.
The lumped modal masses and equivalent gravity loads are given in Table 3.
The columns were assumed to have no limit to their elastic characteristics,' yielding being thus confined to the ends of the beams. This is in accordance with the 1 strong column-weak beam' design philosophy promoted by Code recommendations. Prior to the non-linear analyses, the structure was analysed elastically to evaluate the individual mode contributions. The displacement response predicted by the alternative direct integration scheme, which included all the dynamic degrees-offreedom, served as a benchmark for the comparisons.     to a first mode response only, whereas in curve (b) the first six modes are assumed to participate in the response. The response profile predicted by integration of the first three modes, indicated by the plotted points in Figure 4, is in close agreement with the benchmark test.
A contributing factor leading to the incongruous result whereby the inclusion of additional modes indicates a nonconvergent prediction of the system response relative to the benchmark analysis, is apparent when the modal frequency spectrum is considered.
The additional modal coordinates which participate in the response are drawn from the high frequency region of the spectrum.
The periods of these higher modes approach the length of the time-step interval, 0.01 sees., used in the step-by-step integration, resulting in an inaccurate prediction of the response in the upper modes.
The benchmark response is derived using the same time-step and thus should not be interpreted as representing the precise response time-history.
Integration of the first mode alone, curve (a), results in an overestimate of the structure displacements.
When an allowance is made for member yielding, the higher modes participate more significantly in the overall structure response than in the elastic case.
The shape of the displacement curve is similar to the first mode case, however the effect of the higher modes becomes apparent when the applied earthquake inertia loads cause a reversal in the direct of the structure displacements.
In the process of deriving the transformation basis appropriate to the current structure yield status, it was necessary to resolve the eigenproblem. The mode periods thus derived are indicative of the instantaneous free vibration frequencies of the system. The true period of the system is the time taken by a reference point in the system to complete one cycle of vibration, however this parameter can be difficult to define in response time-histories which include a significant duration of post-elastic yield excursions.
The instantaneous periods for the first and second modes are plotted with respect to earthquake duration in Figure 5.
The instantaneous period is shown to increase threefold compared to the undamped, elastic, freevibration period.
During the periods of intense ground motion the structure responds as a longer period, or less stiff, structure than is the case where the structure is constrained to remain in the linear elastic stress state.
This observation assumes some significance when a comparison is made with the elastic response spectrum philosophy. An acceleration response spectrum for a particular earthquake record is generated by determining the peak acceleration response of a single degree-of-freedom oscillator over a full range of vibration frequencies, assuming a specified level of damping in the resonator.
The response is assumed to remain in the elastic range. Design spectra generate lateral force distributions considerably less than the forces developed from spectral response curves of moderate earthquakes such as El Centro 1940.
The discrepancy is attributed to the fact that the response of the structure is assumed to be elastic, whereas code designed structures respond with deformations which may exceed the elastic limit of some members. A member ductility factor is defined as being equal to the ratio of the force developed in purely elastic response to the member yield force.
For El Centro a ductility factor of approximately four relates the code design spectrum to the earthquake generated response spectrum.
Access to the instantaneous mode period suggests an alternative rationale, in place of the ductility factor concept. The six-storey test structure has a first mode period under elastic conditions of 0.8 sees., and, in the yielded state, of the order of 3.0 sees.
An earthquake response spectrum indicates the peak response only for a particular resonator, giving no information on the elapsed time within a particular record at which this peak occurs.
Elastic time-history analyses carried out for two resonators with periods of 0.8 sees , and 3.0 sees., respectively, using the El Centro earthquake record, showed that the peak response of both resonators occurred at a similar elapsed time.
This enables some direct comparisons to be made using the response spectrum ordinates.
The response spectrum at the 5% damping level for the North-South component of El Centro is given in Figure 6. where the spectral values, Sa, are given as a proportion of g, the acceleration due to gravity. The peak accelerations vary by a factor in the range 3.5 to 4.5.
This simple example demonstrates that the effect of the increase in instantaneous mode periods during intervals of strong ground motion can be sufficient to limit the earthquake inertia loading to close to the code design spectrum level based on the elastic properties of the structure.
Adequate ductility capability must be achieved through the member detailing to permit the member hinges to form and to rotate without any secondary failures in the structure.

CONCLUSIONS:
In this paper the principle of mode superposition, which is encountered commonly in elastic response time-history determination, is extended into the post-elastic domain.
Previous investigators have tended to dismiss the approach as being impractical, however, it is submitted that with some reconstitution of the coordinate transformations fundamental to the approach, modal superposition can be successfully extended beyond the elastic limit and the resultant scheme is competitive with the direct integration Providing the transformation basis is rederived for each change in element yield status and accepting some bounds on the viscous damping specifications, the equilibrium statement in the generalized coordinate system is given by a set of uncoupled equations equal in order to the number of generalized coordinates assumed to participate in the response.
One advantage of this formulation is that the response can be derived by the integration of a selected number of generalized coordinates.
The integration stability limits are related to each individual mode rather than to the single highest frequency mode in the system, with the consequence that the linear acceleration integration scheme becomes a viable algorithm. It is possible to vary the time-step used in the integration of each coordinate, superimposing the contributions when the time-steps coincide.
The time-step used for the lowest frequency, or longest period, mode can be more coarse than that used in the higher modes, with negligible loss in accuracy.
An additional benefit of the modal approach is provided by the time-history of the dynamic properties of the system which is produced in this scheme. Access to the instantaneous frequencies and mode shapes of the structure provides a record of its dynamic sensitivity to a particular earthquake ground motion record.
From the analysis of a six-storey planar frame it was possible to relate its observed peak non-linear displacement response to the displacement envelope predicted by the elastic response spectrum of the earthquake, The traditional approach in this comparison has been to relate the peak elastic spectral response, based on the initial elastic stiffness properties of the system to the maximum observed nonlinear response through a structure 'displacement ductility factor 1 which is generally of the order of 3 to 4.
As a result of the single analysis undertaken in this study a less subjective explanation is suggested.
Because of the influence of element yield on the dynamic properties of the structure, the effective modal periods are lengthened during cycles of yield.
A structure designed for ductile yielding during a moderate earthquake response does not experience the peak elastic spectral response because the influence of yielding displaces the effective period band to the reduced ordinate zone of the spectrum. The structure 1 s non-linear spectrum tends to a flat curve, rather than a peaked curve which is characteristic of fully elastic response.
When the shift in the effective modal period during yield excursions is taken into consideration, the non-linear displacement envelope can be directly related to the earthquake response spectrum without recourse to the concept of a displacement ductility factor.