Collings1978.pdf

Use of trigonometric interpolation for the analytical determination of the dynamic response of linear systems to arbitrary inputs Andrew G. Collings and Leith R. Saunders Department of Civil Engineering, University of Auckland, Auckland, New Zealand (Received 20 December 1977)

Commonly the Ioadings associated with the dynamic excitation of structural systems are expressed in the form of discrete time series. By the use of trigonometric interpolation these inputs may be transformed into continuous time functions. The dynamic response of a linear structural system when subjected to such inputs may be computed by modal analysis techniques. These methods require the solution of a set of uncoupled second-order linear ordinary differential equations with constant coefficients. The trigonometric interpolation function is a finite series of sine and cosine terms and hence it becomes a straightforward matter to solve the set of uncoupled differential equations for the displacements. Solutions obtained in this manner are also seen to be continuous functions expressed as finite series in sines and cosines. The associated velocities and accelerations may be easily obtained by differentiation of these expressions. The method is also free from the errors inherent in step-by-step integration methods, and the errors in the numerical determination of Duhamel integrals. The derivation of the coefficients of the time-series inputs and of the final evaluation of the dynamic response are both executed by a simple recursive algorithm which is fast and accurate. Since the only errors occur in the initial expression of the input data as a finite series of elementary functions--and these errors can be made very small--the method is a novel and attractive alternative to existing procedures. It is illustrated here by the computation of the dynamic response of a skeletal space frame when subjected to earthquake ground motions.

Introduction The linear dynamic response of structural systems, when subjected to irregular loadings, may be adequately represented by ordinary differential equations which have constant coefficients. Although there exist general purpose equation solving 0141-0296/78/010041-12 $02.00 © 1978 IPC Business Press

procedures which may be used to compute solutions to such equation systems, their use may often require excessive storage and computational cost. There are however efficient specialized computational algorithms which successfully exploit the intrinsic properties of such equation systems, sparsity, symmetry, etc. These algorithms rely for their formulation on several distinct

Eng. Struct., 1 978, Vol 1, October

41

Dynamic response of linear systems: A. G. Collings and L, R. Saunders

approaches; in particular, the methods known as direct integration, mode superposition, and frequency domain analysis are frequently used. The texts by Bathe and Wilson I and Clough and Penzien 2 contain expositions of these techniques. In direct integration methods no transformation of the ordinary differential equation system into a different form is carried out. A step-by-step procedure which discretizes the time domain and makes assumptions as to the variation of displacement, velocity and accelerations over each interval is used to integrate directly the equations of motion. With mode superposition techniques however, a change of basis is made for the system of equations. This transformation is such that the original system of equations is replaced by an uncoupled set. These equations are then solved by the use of closed form solutions, numerical evaluation of the convolution integral, commonly referred to as the Duhamel integral, or by numerical integration. With frequency domain analyses use is made of the fact that the Fourier transform of the response is the product of the Fourier transform of the input and the Fourier transform of the pulse response of the system. This paper expresses the input as a finite trigonometric series which reproduces the discrete input exactly and gives a very smooth interpolation between the given data points. This is the trigonometric interpolation method due to Lanczos a-s. It is applicable to any input series. To find the coefficients of the trigonometric series an efficient algorithm 6'7 is used. For long time series the fast Fourier transform can be exploited to advantage. The normal mode method of analysis is applied in this paper. By means of matrix similarity transformations the problem is resolved into the solution of a set of linear second order differential equations with constant coefficients. The complementary functions give the evanescent transient responses, and the particular integrals, corresponding to the trigonometric driving functions, give the steady state responses. The earthquake record is divided into consecutive sections, the first datum of a section being the last datum of the previous section, so that the initial conditions for a section are provided by the end results of the previous section. The solution is obtained in the form of the sum of a finite cosine series, a finite sine series, and a further simple expression. It should be noted that the solution is given as a continuous function of time. However, it is convenient to evaluate the solution at the set of time values in which the original input sequence was given. This can be economically achieved by using the Goertzel algorithm again. The problem of the approximate evaluation of the Duhamel integral has thus been by-passed. An example is given later in which the input to the system is an earthquake acceleration record. The example is explained in some detail in order to indicate the utility of the method when applied to a realistic problem. When an earthquake record is processed to appear as a trigonometric series it is usually found that the higher frequency terms have low amplitudes. Moreover, the structure whose dynamic response is required may have an insignificant response to the higher frequencies. In such situations a worthwhile

42

Eng. Struct., 1978, Vol 1, October

saving in computation can be achieved by simply truncating the interpolation formula at some point N' < N. The formula then loses its interpolation property in that it is no longer exact at the data points. Nevertheless it can remain an excellent approximation, and is indeed the formula of best fit in the least squares sense to the data points, for the specified value of N'. The truncated formula in effect represents a smoothing of the data. The data is supplied as a discrete series with time interval T (commonly a!o or ~!osee.), and the sampling theorem states that this series can only properly represent those frequencies present in the actual disturbance; defined for continuous time, which do not exceed ½T. Choosing N' < N is equivalent to suppressing the higher frequencies in the band (0,½T).

Mathematical basis Trigonometric interpolation Consider. a real sequence X(j) defined for: j = -N,-(N-

l) ..... 0,1,2 ..... N - 1,N

Then it may be readily shown that for the 2N + 1 integer values o f j we may writea: N

X(J~=

~,' C(n)e i"j~/'v

(1)

n = -N

where

1 s C(n} = "2-N~=~-'uX(j) e -i"i'm

(2)

with i = v -/'-Z-II.The primes on equation (2) indicate that the two extreme members of the sums are taken with half weight. To show that the relationships in equations (1) and (2) hold we use the identity:

~-lNO + e-i~u-l~o + ... + ei(N-l)O + ½eilvo = sinNOcotO/2

(3)

which is valid for 0 :~ 2kTt. Assume the existence of equation (1) and consider the following sum: N

s. =

2 ' x(j) e- i.j./~, j=-N

=

N

N

Z'

Y,' C(m) e'j''-"~'/N

j= -Nm=

-N

N

=

Z m=

-N

N

C(m) Z' eV"-"'"/N )=

(4)

-N

If we put (m - n)~/N = 0 with m :~ n the sum over j can be identified with the left hand side of equation (3) with 0 not a multiple of 2n. It follows at once that this sum is zero. If m = n the j sum is 2N, whence S, "= 2NC(n), from which equation (2) follows. In general C(n) is a complex number, but equation (2) reveals that C(0) and C(N) are real. It can also be seen that C ( - n) = C(n), the complex conjugate, which implies that C ( - n ) e -l"~'m is the complex conjugate of C(n)e i"~'/N. Putting C(n) = E(n) + iF(n), then:

Dynamic response of linear systems: A. G. Collings and L. Ft. Saunders N-I

N-I

X()) - C(0) + ~ 2 Re{E(n) + iF(n)}

z(t) = ~ B(n)sin(mtt/NT)

n=l

x {cos njn/N + i sin njrc/N} + C(N) cos Nrq/N

(5)

where T is the time interval between adjacent members of the time series X(j~. At once we have:

z(jT) = Z(j~

This expression may be rewritten: N

X(j) = ~.. A(n)cos(njrc/N) n~O N-I

+ ~ B(n)sin(njrqN)

(6)

n=l

where A(n) = 2E(n) and B(n) - -2F(n). We now have X(fl expressed as the sum of an even function and an odd function, defined for integer j. It is easy to express X(j) as the sum of an even sequence and an odd sequence. Defining:

Y(j) = ( X(J3 + X ( - j ) ) / 2

(7)

Z(33 = ( X(3) - X(-j~)/2 so that Y(j~ = Y(-j~ and Z(j~ = - Z ( - f l . X(f) = Y(j) + Z(j3

( - N <<.j<<.N)

(14)

Provided that the boundary conditions Z(O) = Z(N) = 0 hold, z(t) is the interpolation of a hypothetical function specified at the finite points iT. A general comment on trigonometrical interpolation can be made. Consider the situation where a function f(t) is given as a continuous function, but the integrals required for a standard Fourier analysis cannot be conveniently determined. We may then use the sequence f ( j T ) and trigonometrical interpolation. A lucid exposition of the closely related topics of trigonometric interpolation at equidistant points, finite Fourier series, band-limited functions, aliasing and the sampling theorem is given in Cooley et al. 18 Parallel to equation (13) we may write:

Then:

y(t) = ~. A(n)cos(nnt/NT)

( - N <~j <~ N)

n=O

(8) with

In addition, equation (2) may be expanded: 1

(13)

n=l

(15)

y(jT) = Y(j)

ht

( - N <<.j<~N)

C(n) = ~--~j =~-'NX(j~{cos(njn/N) - i sin(njrc/N)} = 2

x(o) +

{x(j? + x ( - j 3 } ./=1

Equation (8) is now replaced by the interpolation formula:

cos(nj~/N)

x(t) --- y(t) + z(t)

+ ½{X(N) + X( - N)} cos nn

The question can be raised: is it possible to adapt the discrete Fourier transform to produce an interpolation formula? If this were the case then the powerful fast Fourier transform method s'9 could be applied directly. It is now shown that the discrete Fourier transform cannot be used directly. The discrete Fourier transform formulae are:

- iN-1 ~ {X(j) - X ( - f l } sin(nflt/N) 1 ./=1

=

~ NL[½X(O) + ~-' j=l

Y(j3cos(njn/N)

+ ½Y(N)cosnn - i N~,- I Z(./)sin(njn/N) 1

(9)

j=l

1 N-I

A(n) = -~ j~-o F(j)e- 2,a,,j/N

Equations (6)-(9) may be collated by presenting the following formulae:

A(O) = ½A'(O)

A(N) -~ ½A'(N)

(18)

n=O

N~=o

(1 ~< n ~< N - 1)

(17)

N-I

F(j) = ~ A(n)e 2"~"j/N

A'(n) = 2 ~, y(fl cos(nj~/N) A(n) = A'(n)

(16)

(10)

where F(j) is a real sequence, W e consider a continuous function f(t) such that: N-I

f(t) = ~, A(n)e 2"t=mr

(19)

n=.O

and: 2N-1

B(n) = -~ ~ffil Z(J~ sin(njn/N) (1 ~< n ~< N -

1)

(11)

So far we have relationships between the discrete sequences X(fl, Y(j), Z(j), A(n) and B(n). We now investigate the possibility of extending these ideas to the continuous domain. We have:

N

N-1

Z(fl = ~ B(n)sin(njn/N)

with f ( j T ) = F(j). It may easily be shown that for F(fl real, A(N - n) equals the complex conjugate ,,](n). But since e 2 f l ( M - n)IINT is not the complex conjugate of e 2"~'mr, as t not an integer multiple of T, the n and N - n terms of equation (19) are not complex conjugates, from which it follows that their sum is, in general, not real. This rules out a straightforward use of the discrete Fourier transform. But the interpolation formula derived from equation (I), namely:

(12)

n==l

where the B(n) are given by equation (11). Consider the continuous function defined by:

x(t)=

~

C(n)e i='mr

n = --N

is such that x(t) is always real. Lanezos 3'4 derived the trigonometrical interpolation from fitting the data

Eng. Struct., 1978, Vol 1, October 43

Dynamic response of linear systems." A. G. Collings and L. R. Saunders points with particular orthogonal functions a n d a least-squares criterion, and the interpolation formula represents a limiting case. At this stage we consider a point that has been neglected so far. Equation (13) necessarily implies that Z(N) = 0. However, the application of equation (7) will, in general, yield a non-zero Z(N). This difficulty is overcome by subtracting a straight line through the origin from Z(j). Explicitly:

Z~(J3 = Z(j) - (j/N)Z(N)

(20)

We thus obtain the condition Z~(N) = 0, and consequently we apply sine interpolation to ZI(j~, and roplac~ equation (20) by: x(t) = y(t) + zt(t) + (t/N73Z(N)

(21)

approximately by the law n-3 where n is the nth coefficient. This rate of decrease also applies to the coefficients of a sine series of an odd function which vanishes at both ends of the interval. It follows that truncation will only produce small errors in these situations and our separation of the input function into the sum of an even function and an odd function accordingly ensures accuracy in the interpolation. Equation (24) is used in this paper to represent the acceleration record. In the example the physical conditions are such that the input is a linear combination of an earthquake acceleration record and the corresponding velocity function. The latter is obtained by integration of the interpolation equation (24). The input function of the example therefore has the form:

Using Equations (13), (15) and (21) we write:

N

a(r) = ~,, ~l(n)cos(mt3/NT)

N

x(t) = ~, A(n)cos(nm/NT)

n=l

N

n=O

+ ~ ~2(n)sin(nnr/NT) + (~ + Jr + / ' I t 2 (26)

N--I

+ ~,, B(n)sin(nrtt/NT) + (t/NT)Z(N)

This formula applies for - N T ~ t <~ NT, so that the point t = 0 is at the centre of the range. It is more convenient in practice to have the point t = 0 at the beginning of the interval, especially when initial conditions are involved. We make the transformation:

t = r - NT

(23)

and since

cos(nm/NT) = cos(nnr/NT) cos nn sin(nnt/NT) = sin(nnr/NT) cos nn then x(t) can be replaced by: N

X(r)

=

Y.

A(n)cos(nnr/NT)

n=0 N-1

+ ~ B(n)sin(nrcr/NT) n=l

+ (r/NT-

1)Z(N)

n=[

(22)

Solution of differential equation The normal mode method of deriving structure response involves the solution of a set of scalar differential equations typified by: w+ 2tiff + co2ff = P(Q

(27)

where P(3) is the driving function, or the input to the particular normal mode. The dots signify derivates with respect to 3. The interpolation formula for P(r), (cf. equation (26)), involves terms of the form cos ~ , sin ~r and p + q3 + r32, where = = nn/NT. It is easy to show that these expressions have the following three particular integrals respectively: (o92_ ~t2)cos ar + 2fl~tsin r,3 (o92 _ ~2)2 + (2fl~)2

(28)

(ogZ _ ~2) sin ~r - 2fin cos ~3 (o92 _ ~2)2 + (2fin)2

(29)

(24)

p 0)2

where:

2flq + 2r + 8f12r ( q o94 7 + 0) 2

4fir I 034 ]

r 2 r + -0 )-2 3

(30)

/

A(n) = ( - 1)"A(n) /~(n) = ( - l)"B(n)

(25)

with 0 ~< ~ ~< 2NT. Thus the transition from equation (22) to (25) requires little more than changing the sign of every second coefficient in the finite series. This is an appropriate point at which to discuss the accuracy of trigonometric interpolation. Lanezos s has shown that the error bound of a trigonometric interpolation is the same as that for a finite Fourier series of the same order. Though the coefficients may differ a little the approximations are of similar accuracy. Incidentally, equation (2) may be regarded as the trapezoidal rule applied to the integrals with would be used to obtain Fourier coefficients. Because of the similarity between the two approaches it follows that conditions on the time functions which make for accuracy in a Fourier series also apply to an interpolation series. For an even function the coefficients of a cosine series decrease

44

Eng. Struct., 1978, Vol 1, October

The complementary function of equation (27) has the form: e - n ( K l cos c3 + K2 sin cr)

(31)

where c 2 = o)2 _ f12. It is assumed that the system is under-damped so that fl < to. Equations (28)-(31) show that the solution takes the form: N

if(r) = ~ t~t(n)cos(nn3/NT) n=t N

+ E 0,(n)sin(n, r/Nr)+ : + Or + a3' nml

+ e-n'(Kt cos cr + K2 sin or)

(32)

The constants KI and K2 are determined by initial conditions on displacement and velocity. Equation (32) gives the displacements as eontinuous functions of 3, and is therefore a general analytical representation of the response. It is convenient in practice to evaluate

Dynamic response of linear systems: A. G. Collings and L. R. Saunders the response at the same discrete time intervals j T for which the earthquake acceleration record is given. Prior to this it is desirable to separate out the complementary function part of equation (32). Thus, we define: if(z) = ~p(~) + ~,(~)

(33)

where ~p(~) and ~,(~) are the particular integral and complementary function. Also: I~(./) = i f j T ) = l,Vp(]') + I,V=(j')

(34)

where

(35)

The terms e-l~rcos(jcT) and e -~#r sin(jcT) are the real and imaginary parts of e~-~+~r, and this permits an •flicient recursive algorithm to be adopted for the computation of I,V=(j). Define: ~j = e~-~+~lir = e-i~r{cos(jcT) + isin(jcT)}

(36)

Then ~+ ~ = ~j + r / ~ where

Fl = e - # r c o s c T - 1 + i e - # r s i n c T

n=O

N- l

(42)

Tp(fl = ~ Oz(n)sin(njg/N) + Q Tj n=l

The functions S~f) and Tp(fl are respectively even and odd functions of the integer j. The response at the points z --- kT, 0 ~< k ~< 2N is given by: # d k ) = S , ( N - k) - T , ( N - k)

l~p(2N

-

k)

=

Sp(N

-

)

k) + Tp(N

-

k) ;

(43)

for 0 ~< k ~< N. On adding l~p(k) and l~'dk) we have the total response ifT) at the points "c = k T, 0 <~k <~ 2N.

l~'~(j) = e - ~ r { r ~ cos(jcT) + Kz sin(jcT)} (0~<j~<2N)

N

s,(h = Y~ O~(n)cos(nj,~/N)+ RT'j ~

(37)

and ~o = 1

The recursion is carried out for j = 0, 1, 2. . . . . 2N. This method involves the use of the computer library functions only three times, to compute e -#r, cos cT, and sin cT. Despite the recursion the accumulation of rounding errors is negligible. Thus, we obtain ¢¢,(3),0 <<.j<~2N. Since fl can be small, e -#r can approach unity, and then this algorithm is advantageous t°. The particular integral is:

Evaluation of finite trigonometric series To evaluate the sums in equations (10), (11) and (42) it is natural to use one of the economical algorithms which have appeared in recent years. For long series it is possible to adapt the fast Fourier transform method used for discrete Fourier transforms. This will be elaborated on below. In situations where the given input series has about 100 terms or less the splitting into even and odd series of about 50 terms allows a simpler approach to be adopted. Lanczos 4 has shown that the symmetry and periodic properties of trigonometric functions can be exploited to reduce the number of multiplications by a factor of four, against direct evaluation of the sum. However, it is still necessary to generate a table of cosines and sines. We consider that in most applications the method of evaluation of finite trigonometric series due to Goertzel 6 has much in its favour. This method uses recursion and does not require the prior setting up of tables. It is a straightforward algorithm and is very easily implemented by an efficient computer program. In our application we used the Goertzel method in the form given by Hamming 7, with the difference that our N replaces Hamming's 2N - 1. Define: N

C=

N

~p(z) = ~ Ol(n)cos(nnl:/NT)

N

~ a~coskx

and

S=

k=0

n=l

~ aksinkx k=l

Also let: N

+ ~ t~,(n)sin(nn,/NT)

Uo ----0

n=l

+ P + 0.~ + 1~ 2

0 <~• <<.2 N T

(38)

Ut = a~v

Um= 2U,,_1 c o s x -

U,,-2 + ate-,,+ 1

(44)

for m = 2, 3.... , N. Then At this stage we make the change of variable r = N T + t so that - N T <~t <~NT. Consequently:

C = U N C O S X - - U N + 1 "~- ao

wp(t) = ~, Ol(n) cos(nnt/NT) n=l. N

+ ~ Oz(n)sin(nm/NT) + P + Qt + Rt 2

(39)

where 0t(n) = (-l)'/~l(n)

and

02(n) = ( - 1)"/~2(n) (40)

Putting t = jT, - N ~<j ~< N, P = 0~(0), and grouping even and odd functions, yields: w,(~ =

where:

(45)

S = UN sin x

N

sd~ +T,(~

In our application x is of the form nn/N. To obtain the sum ~ ' of equation (1(3) we first put Y(0) = Y(0)/2, Y(N) = Y(N)/2. For the sine summations the algorithm is applied with U1 = 0. A single computer procedure was developed for all of the summations in equations (10), (11) and (42), (an A L G O L procedure is supplied in the Appendix). An alternative approach, attractive for long series, is to relate the trigonometric series to the discrete Fourier transform, and apply a fast Fourier transform computer routine. Consider equation (2):

(41)

C(n) = 2 ~

N=~_'tcX(J)e-'~'m

( - N <~n <~N)

1

Eng. Struct., 1978, Vol 1, October

45

Dynamic response of linear systems: A. G. Collings and L. R. Saunders Put k = j + N and define )~'(k) = X(j~, (0 ~< k ~ 2N). Also let M = 2N then we may write: 1

M

=

C(n) = Vi y"' X( k )ei . . e-2i'"~/M a,, k = O

Vl M-1 : ( - - 1 ) I'M ,~o "Y(k)e- 2"':'/"

the end of the data of section q - 1. By the time that the overlap region has been transversed the transient terms have become negligible, and may therefore be omitted from the solution at the outset. This technique has not been adopted in the example. To be precise, there is the overlap of a single term, in that the first term of the section is taken as the last term of the previous section, in order to permit the specification of initial conditions for each section.

(46)

Integration of interpolation formula The first term inside the square brackets is the discrete Fourier transform of X(k), which we may refer to as D(n). We now let n lie in the interval 0 ~< n ~< M - 1 and interpret C(M - n) as identical to C( - n). Thus:

C(n)=(-1)"[D(n)+~---~{2((M)-.~'(0)} 1

(47)

A further simplification is entailed when a straight line is subtracted in accordance with equation (20). This has the effect of rendering X(M) equal to X(0) so that the term 2 1 { , y ( M ) - .X(0)} vanishes. The original (7(- N) is not covered by the above, but it is a real quantity equal to C(N) and is therefore determined. Equation (47) caters for equations (10) and (11), and it is easy to also express the sums (equation (42)) in terms of a discrete Fourier transform. The discrete transforms are those of a real series, so that additional economy can be obtained by using a fast algorithm devised especially for such series s-~°.

Sectioning and overlapping The input record is usually a long series, and the possible advantages of dividing the record into consecutive sections must be considered. To be specific, suppose the earthquake record to be supplied as 500 values at intervals of a~ or ~ s e c . It would be possible to find a single interpolation for the entire record. However, such a procedure would ignore the frequency spectrum of the record ~x. For example, suppose that the amplitudes are small for frequencies below 0.5 Hertz (cycles per second). Then, if the data are spaced at intervals of ~ sec. a section of 80 items would deal faithfully with the low frequency part of the spectrum. In general, the lowest frequency representable by a section of M data with interval T is I/MTHz. In the example given below the input record is CALTECH C2 ~2 for which the choice of M - 80 is quite appropriate. In fact we chose M = 83 so that 6 sections would adequately cover the record. For trigonometric interpolation, M = 2N + 1 implies N = 41. With such a small value of N the accuracy of the Goertzel algorithm is quite adequate. Transient terms appear in the solution, and these will occur after the start of each section. It is necessary to include these terms since they inevitably arise in the interpolation method. In situations where the damping is large the transient terms will decay rapidly. The possibility then arises of using overlapping: the start of section q is made on data a few time intervals before

46

Eng. Struct., 1978, Vol 1, October

In some engineering problems, such as the example given in this paper, it is necessary to have available the input velocity record as well as the acceleration record. Although the velocity record may possibly be supplied along with the acceleration record it is easily produced by integration of the trigonometric interpolation formula. The accuracy of the resulting formula can be discussed by noting that the error of the original interpolation is of the form u(t)sin(Nnt/T), since this function vanishes at the points jT, for j an integer. If we let f ( t ) be the true function and f*(t) be the approximation provided by trigonometric interpolation then:

f(t) = J*(t) + u(t)sin(rtt/T)

(48)

Integrating this expression by parts, and by passing the use of a definite integral by the inclusion of a constant of integration K we obtain:

F(t) = F*(t) - (T/n)u(t)cos(rttiT) + 0(T 2) + K

(49)

in which capital letters represent integrals. Hence, for N reasonably large the error bound is given by:

le(t)l ~< (Z/rOlu(t)lm=x Since T/n is small the formula resulting from integration of the interpolation has greater accuracy than the interpolation itself, and is therefore an excellent method of obtaining the integral of a function supplied at discrete points.

Illustrative

example

Specification of problem The time-dependent dynamic response of linear complex structures satisfies ordinary differential equation systems of the following form:

M)~ + Cf~ + K X = P(~)

(50)

In equation (50), M, C and K are respectivoly the generalized inertia, damping and stiffness matrices of order m where m specifies the number of degrees of freedom in the discretized system. X is an m-vector of generalized displacements, measured in an inertial frame of reference, with ~" and ,~" as its first- and second-time derivatives. The m-vector P(z) is a set of generalized time-dependent forcing functions. If we apply a binary partitioning to equation (50), the resulting substructured system is one which is frequently encountered in structural dynamics problems. The partitioned system may be written as follows:

Dynamic response of linear systems: A. G. Collings and L. R. Saunders

[M,I

M21

MlqFX, l F ll

c, IFX, l

M22JLX2j + LC21

[ KI'

+ LK21

"qF

J= -Kil-lKl2

c22jLX2j

'l

K22_JLX2j = LP2(oj

i,,l

where the vectors of generalized displacements X1, X2 and their associated time-derivatives have orders of n and p respectively, with n + p = m. The importance of equation (51) is apparent from the fact that a judicious ordering of the rows and columns of equation (50) will allow a wide range of dynamical problems to be handled. For example, suppose we let one partition of equation (51) contain those elements for which the displacement, velocity and acceleration vectors have known prescribed time-dependent values, while the other partition contains those elements for which the displacement responses are required. This type of problem often occurs in earthquake engineering situations. In particular, we may let X2, X2 and ~'2 be a set of displacements, velocities and accelerations that have prescribed values, such as those which could be computed from the accelerations recorded at the foundations of a structural system when subjected to seismic disturbances. PI(O and P2(r) will typically be vectors of dynamic loading associated with the nonprescribed and prescribed degrees of freedom respectively, with PI(O prescribed. In the absence of dynamic loading associated with the non-prescribed degrees of freedom, such as that which may be due to wind or waves etc., PI(O will normally be a null vector. P2(r) might sometimes be regarded as a set of time-dependent dynamic reactions. The prediction of the time-dependent dynamic reactions. The prediction of the time-dependent response of the non-prescribed degrees of freedom depends upon the solution of the norder differential equation system which results from the expansion of the first line of equation (51). Hence: MrllXl + CII,~'I + KlIXl = PI(z) - M12)/'2 - C I 2 X 2 - K I 2 X 2

(52)

where the right-hand side is a known function of 3. The system which equation (52) represents does not cater for the complicated matter of structure-foundation interaction. Rather, it is concerned with a system which assumes that the structure is connected to a rigid localized region (here called the 'foundation'). This foundation under the action of seismic disturbances may typically be subjected to prescribed displacements, velocities and accelerations. W e observe that the inertial frame, in which equation (52) has been written provides a convenient system for measuring the prescribed displacements, velocities and accelerations to which the system may be subjected. However, a local non-inertial frame of reference Y, rigidly attached to the foundation, from which the time-dependent response of the nonprescribed degrees of freedom can be measured, is of practical importance. A transformation between the non-inertial and inertial coordinate systems may be written as follows 1":

Y = Xl - JX2 (53) where J, a transformation matrix, may be shown to have the following construction I a:

Using this relationship for J and substituting equation (53) into (52), with Pl(~) null and rearranging yields:

M I I ~" + CII ~" + KII Y = R)(2 + S,~2

(54)

where R = -(MI2

+ Mll J)~

and

t

(55)

S = -(C12 + Cll,J) J The vectors of prescribed velocity and acceleration ,~2, )/'2 will invariably be highly irregular in nature and are often presented in the form of digitized records where the digitization interval Tmight typically be r~-o, ~ , or a~ see., etc.

Solution using trigonometric interpolation We divide the given acceleration record into segments as explained previously. We assume that the damping in the system is low and consequently transients will form a significant part of the solution. Accordingly, we do not use the overlapping technique described earlier in the paper. Each section has 2N + 1 terms and the 2Nth point of section q - 1 is assumed to coincide with the 0th term of section q. In this way boundary conditions can be handled conveniently. Our first step is to consider section q and write a trigonometric interpolation expression for the discretized record of accelerations. Accordingly, we have by equation (24): N

.~q(~)= ~ ,:l(n)cos(nn~/NT) n=0 N-I

+ ~ B(n)sin(mtr/NT) + ( r / N T -

1)(D/2) (56)

n=l

where ,~(n) and ~(n), 0 ~< n ~< N, are the frequency coefficients of the trigonometric interpolation for the qth set of real time series acceleration data. The time variable T, which operates over the qth set of real time data, has the range 0 ~< ~ ~< 2NT. D = Z(N) is the final member of the odd time series components. If we integrate equation (56) then an expression for the associated prescribed velocities for the qth set of time series data is obtained: 7r / , = I

n

-Y. n=l

+ n

+ :i(o)-

n~l

~+4--~.~ + ¢

(57)

The constant of integration (~ may be evaluated by assuming continuity of velocity over the boundaries of the sets of real time data. In particular, we may assume that: ~q- t(2NT) = .¢q(0)

(58)

Thus, if we set z = 0 in equation (57) we then have:

C = Jq(0)= .~-'(2NT)

(59)

Eng. Struct., 1978, Vol 1, October 47

Dynamic response of linear systems: A. G. Collings and L. R. Saunders Returning to equation (54) we may now substitute on the right hand side for the vectors of prescribed accelerations and velocities by putting: ,~2 = ~a(r)e

(60)

terms involving powers of 3. Standard ordinary differential equation theory affords convenient solutions in such cases. We may rearrange the right hand side of equation (64) so that it has a similar form to that of (26). It may be shown that:

(b,(n) = r,A(n)

where E is a p-vector which has values 1 or 0 in accordance with the decision as to which degrees of freedom have the prescribed acceleration. Hence on substituting equations (56), (57) and (60) into (54) we obtain:

#

zJ(n)cosn/~z +

Also

+ -s, - ii(n)

=

It

B(n)sinng~ n=l

,4(n) sin n#z n

n

(o)

cos n#~ +

~ n= I

?!

o2 +

+ /](0)-~- ~+4--~,

(67)

K l l $ - Mll4~fl = 0

/'1

E

]

(61)

(62)

where M l l , K l l are positive definite symmetric and is a diagonal matrix of real eigenvalues. Now, if we choose CI 1 to be constructed from a linear combination of KI 1 and M11 as in Rayleigh's model la then the similarity transformation that diagonalizes M l l and/(11 will do the same to C l l . Let: Y = ~bW

s,(# ,----1 n

err

-

mrr

k,,

(68)

J = r,~--~ + s, .x](0) -

(69)

/~

(70)

O = S,4N T

0t(n) = ~l(n)Ls(n) - ~2(n)L2(n) ~ LI(n) (02(n)L3(n) - (Ot(n)L2(n) (1 ~< n ~< N) 02(.) =

(71)

Ll(n)

where Lx(n) = (co,2 - (n/~)2)2 + (2fl, n#) 2" L2(n) = 2fl, nu

k, = r, 5~(z) + s,~(r)

(64)

mrr

where m,,. c,, and k,, are the rth coefficients of the diagonalized mass, damping and stiffness matrices. The terms ~q(,) and ~(~) are given by equations (56) and (57) whereas the coefficients r, and s, are determined from the following relationships: 1 rntr it= 1

+ C}

It is now possible to write the solution to equation (64) by using the standard forms for the complementary function and particular integrals given by equations (28)-(31). This procedure results in an expression similar in form to equation (32). In particular:

(63)

where W is measured with respect to a set of normalized coordinates. Then the rth uncoupled scalar ordinary differential equation is obtained from equation (61) as follows. Substitution of equation (63) into (61) with premultiplication by q~r, and rearrangement gives: -

se A(N) /~ N

--

In addition n=l

where g = n/NT. Under suitable conditions, equation (61) may be uncoupled by a similarity transformation that uses the matrix of eigenvectors $ of the homogeneous system:

~i,, + - - ~,, +

(l~n~
n

with 2(N) =

+ S

(66)

~I(N) = r,.4(N)

¢;2(.) R

(1 ~< n ~< N - 1)

n

with

X l 1 ~b - M 1 l q ~ n = 0

=

s, B(n)

j= l

(65)

(72)

and

La(n) = o~,2 - (n#) 2 In equations (72), to,, the rth natural frequency of the system, and fl, are given by k,,/m,, and c,,/2rn,,

respectively. The expressions for/5, (~ and/~ in equation (32) may be shown to be given by:

p= d

1 rtlrr k = l

2=1

~here e~ is the flh element of the vector E and rki and ski are the kj elements of the matrices R and S given by equation (55). A solution to equation (64) may be obtained if one notes that the expressions for ~(¢) and ~q(~) are composed entirely of sine and cosine terms, plus a few

48

Eng. Struct., 1978, Vol 1, O c t o b e r

o,,2

(73)

Equation (32) also contains two constants Kt and K,, these may be determined from the initial

Dynamic response of linear systems: A. G. Co~tings and L. R. Saunders the purpose of analysis. The element stiffness and inertia matrices have the following forms:

conditions associated with the solution at the beginning of the qth block of data. In particular:

~$- t(2NT) = g$(O) (74)

and

"

g,,~- t(2NT ) = ~b,q(0)

If we apply these initial conditions then we obtain the following expressions for K t and K2: N

Kt = g'~q-I)(2NT}- ~ dr(n} - / ~ ii = |

f BrEBdv

gi =

(75)

We may take the time derivatives of equation (32) to obtain expressions for the velocities and accelerations. A transformation from the normal coordinate system, used for the purpose of the solution, to a set of local coordinates may be achieved by use of equation (63). Test structure The numerical investigation was confined to tests on a steel membered skeletal space frame in the form of a prismatic hexagon. (cf. (Figure 1)). Equation (50) gives the time dependent linear dynamic response of such a system when it is subjected to dynamic loads. In equation (50) K and M, the generalized overall stiffness and inertia matrices of the system, are constructed from the individual stiffness and inertia matrices of the elements into which the structure has been divided for

M~= fp~Ndv v

where the suffix j denotes the jth element and p is taken as the mass density per unit volume of the element. The matrices B, E and N are defined by the following relationships:

e=Bd I Ee

a

=

(77)

f=Nd

where ~r and e are vectors of stress and strain for the e l e m e n t , f is a vector which represents the displacement of an arbitrary point within the element in terms of a matrix of shape function N, and d is a vector of generalized nodal coordinates. A consistent inertia matrix t4 formulation was adopted by which the same matrix of shape functions N used to generate the element stiffness matrices was also used to generate the inertia matrices. Dissipative forces proportional to the relative velocities of the system were presumed to exist, and hence the damping was represented by a matrix C operating on the velocity vector. A Rayleigh model for damping was used and it was assumed that15: C = :~M + flK

dof3

(76)

(78)

where ~t and fl are constants. The values of :t and fl were determined from a modal proposed by Morduchow 16 in which: ~t = 4n(Ti2 j - TjA,)/(Tf - T~) fl = T~Tj{Tfit,- T~2~)/n(T 2 - r~)

(79)

where Ti and Tj are the periods of the ith and jth modes, (i < j) which have ratios of critical damping associated with them of ,;-i and 2/(each taken as 0.004 in the dynamic analyses). For the purpose of the investigation i and j were set to values of l and 2 respectively, corresponding to the first and second modes of the undamped system. The periods T~ and Tj of the system were determined to be 0.6483 and 0.5946 sec. respectively. The structural system was stiff in the sense that the ratio of the maximum to minimum periods was 264.3.

Discussion o f errors

Figure 1 Steel membered skeletal space with seismic inputs at foundation nodes. System comprises 72 d.o.f.; 6 per node

The error function for the trigonometrical interpolation method does not vary much over each block of input data. It follows that when the input function is small the relative error is large. This is in contrast to the situation with the balanced Euler (trapezoidal rule) numerical integration method where the global error is some function of time and may be shown to be proportional to the numerical integration step length squared ~7. In practice it is of interest to evaluate the dynamic response of the structure when

Eng. Struct., 1978, Vol 1, October

49

Dynamic response of linear systems." A. G. Collings and L. R. Saunders the input is of large magnitude, and it is then that the error of the trigonometric interpolation method is negligible. Figures 3 and 4 show the time dependent displacements of degtee of freedom I and degree of freedom 20, as computed by the balanced Euler and trigonometric interpolation solutions. These solutions are for the case when the structure shown in Figure 1 is subjected, to 2 sec. of ground motion as recorded by the Caltech C2 acceleration record (cf. Figure 2a). In Figure 2a straight lines connect the individual data points on the graph. In Figure 2b the trigonometric interpolation of the earthquake record is given as a smooth continuous curve that passes through the individual data points of the acceleration record. Both Figures 3 and 4 show that the two solutions are almost coincident. Table 1, which records the values of the displacements of d.o.f. 1 and 20 for the Euler and trigonometric solutions, clearly shows that when the input function is relatively small the differences between the two solutions are most pronounced but when the input function increases in magnitude these differences become much smaller, In the balanced Euler method an integration time step h of T/15 ( T = 0.025 sec.) was chosen. This selection was made on the basis that any algorithm used for computing a numerical solution to equation 154) using a time step h should be checked by repeating the entire computation with a time step of hi2 and then comparing the results of the two computations. The Euler solution also required the integration of the acceleration record to obtain velocities. The method adopted was to apply linear interpolation between the data points of the acceleration record and then integrate once to obtain a piecewise-quadratic function for the velocities. This simple approach is obviously liable to errors; for example, there will be a finite velocity at the termination of the acceleration record. In the example given in this paper the trigonometric interpolation method also required velocities, but these are readily found by integration of the interpolation formula of accelerations. Since the interpolation formula is a smooth function in the form of a trigonometric series it can be integrated analytically to obtain the velocities very accurately. The coefficients of a trigonometric interpolation are of the same order as the corresponding Fourier series coefficients. Lanczos 4 has shown that when a continuous even function is expanded into a cosine series the coefficients decrease with speed n-2. If the function is odd and vanishes at the end points of the interval the coefficients of a sine series decrease with speed n-3. These considerations governed our decision to separate the input into an even and an odd function.

Conclusions Trigonometrical interpolation represents the input to the system as a series of elementary functions. It is therefore a straightforward matter to obtain the response of a linear system in the form of a finite series of elementary functions. The solution is achieved without recourse to numerical methods, from which it follows that the accuracy depends only on the accuracy of the interpolation process. We have ensured a high

50

Eng. Struct., 1978, Vol 1, October

06 o5

a

o4 03 ,g" o z

~ o.~ ,-:E

o

,/ -o.I ~ -oz ~ -o.3 ~ -o4 < -0.5 -0,6 -o.7 -o.s -o.9~-----~ io or 0.6

20

L

30 40 50 Time,(x O.O25sec)

60

70

60

70

80

b

05 0,4 03 ,~ o2 :~ o1 ~~_ o '/ la

,c-0.1

~ -o.2 _~ -03 ~ -o. 4 <

-0.5 / - 0 . 6 ~ -0.7 -0.8

-0"90

I0

20

30 40 50 Time, {x 0.025 see)

80

Figure 2

(a) Caltech C2 earthquake ground acceleration record; (b) Trigonometric interpolation of Caltech C2 earthquake ground acceleration record

degree of accuracy by separating the input record into an even and an odd function and using a cosine series to represent the former and a sine series to represent the latter. The error function of the trigonometric interpolation is an oscillating function whose amplitude does not change much during the entire time interval. Consequently relative errors are large when the function to be approximated is small, but are small when that function is large. It follows that the error is small in the domain of interest, where the earthquake has achieved its 'build-up'. The problem which we chose for demonstration, which is in fact a typical problem, required a

Dynamic response of linear systems: A. G. Collings and L. R. Saunders 1.0

0.8

09

0.7

0.8

0.6

0.5 0.4 0.3 0.2

0.7 0.6 0.5

0.4 ?o

'

0.3

~

?o

0.2 ~/ 0.1

I

0.1

,,~

N

o

-0.1

K

K

-0.2

C c

c

i -0.1 -~ -0.2 i5 -03 -0.4 -0.!5

~

-O.3

d

-O.4

~ -0.5 • -0.6

-0.7

-0.6

-0.8

-0.7

-0.9

-0.8

-I .0

-0.9

,~

z'o

~o

4*0

5o~--~o

70

80

i

0 " - ~ 1 0 " ~ J ' - -20 ----L

Figure 3

Displacement of degree of freedom 1. ( O ) , trigonometric interpolation solution; ( x ), balanced Euler solution

Table 1

30

L 50

~o

do

¢o

.80

T i m e , ( x O . 0 2 5 se¢ )

T i m e , ( x 0 . 0 ; ) 5 sec)

Figure 4

Displacement of degree of freedom 20. ( O ) , trigonometric interpolation solution; ( x ), balanced Euler solution

Displacement of degree of freedom 1 for balanced Euler and trigonometric interpolation

solutions Displacement of degree of freedom 1

Trigonometric Caltech C2

interpolation

Time

acceleration record

Euler solution

solution

(unit, 0 . 0 2 5 sec)

M/sec =

(M)

(M)

1

- 1.4694.10 -2

1.5608.10 -6

1.7341 . 10 - s

2

- 4.4686.10-

=

1.5420.10-

6.1353.10 -=

=

6.0446.1

3

- 6.3918.10-

4

-5.5234.10

-2

5

-4.5185.10

-2

10

2.9601 . 10 - =

20

- 3.2641.10- '

30

5.0161.10-'

40

- 1.6877.10- '

50

- 1.2145.

s o- s

1.2340.10- =

1.4557.10 -4

2 . 2 9 4 1 . 1 0 - '=

2.6102.10 .4 -1.4469.10

-`=

- 2.0873.1 o- =

3.5887.10- 4 8.8425.1 o- =

-

-2.2091.10

3.4793.10 -= - 1.3537.10- =

-=

3.6312.10 -3 1.4491 . 1 o - =

-

9.1261 . 10 - =

9.1631.10 -= - 1.2265.1 o- =

60

2.2252.

- 1.2326.10 -=

70

- 1.5348.

- 2.8745.10- =

2.9046.1 o- =

2.9445.1 o- =

2.9694.10- = 3 . 9 5 5 1 . 1 0 -=

75

1.1 2 3 2 . 1 o - 1

76

3.7180.10 -I

3.9177.10 -=

77

1.3493.

4.5697.10 -2

4.6127.10 -=

78

7.6528.1 o- I

4.8281 . 10 - 2

4.8731.1 o- =

79 80

-5.1954.10

-1

1.3723.10-'

4.6913.10 -=

4.7237.1 o- =

4.2511 . 10 - 2

4.2838.1

Eng. Struct.,

1978,

o-

=

Vol

1, October

51

Dynamic response of linear systems: A. G. Collings and L. R. Saunders knowledge of the velocities of the input in addition to the accelerations. The method of trigonometric interpolation lends itself well to the term-by-term integration of the acceleration interpolation series. Since the final result of the method is a formula for displacement as a finite series of cosine and sine terms together with a quadratic term, it is a simple matter to differentiate the result to obtain velocities and accelerations. The method entails an eigenvalue analysis and the determination of eigenvectors. By far the greater computational effort is involved in this analysis. This is so because of the use of a very efficient algorithm for obtaining the interpolation function and the final evaluation of the displacements. A reduction in the overall computational effort may be made by reducing the order of the uncoupled system by selecting only the significant modal shapes for the eigenvector transformation. In addition we may also simplify the solution by truncating the high frequencies in the input series, in which case the input function ceases to be an interpolation function and becomes a least squares fit to the data points.

Acknowledgements The authors thank G. J. Tee and W. J. Walker of the Department of Mathematics, University of Auckland, for their interest and helpful contributions. They also thank the Winstone Centenary Educational Trust for financial support for this research project.

References 1 2 3 4 5 6 7 8 9 l0 II 12 13 14 15 16

52

Bathe,K. J. and Wilson, E. L. 'Numerical Methods in Finite Element Analysis', Prentice-Hall, Englewood Cliffs, USA, 1976 Clough, R. W. and Penzien, J. 'Dynamics of Structures', McGraw-Hill, 1975 Lanczos,C. J. Math. Phrs. 1938, 17, 123 Lanczos,C. "Applied Analysis', Pitman, London, 1957 Lanczos,C. "Linear Differential Operators', Van Nostrand, London, 1961 Goertzel, G. Am. Math. Monthly 1958, 65, 34 Hamming, R. W. "Numerical Methods for Scientists and Engineers', McGraw-Hill, New York, 1962 Cooley, J. W. et al. J. Sound. Vib. 1970, 12, 315 Bergland, G. D. Commun. Assoc. Comput. Much., 1968, it, 703 Singleton, R. C. I.E.E.E. Trans. Audio Electroacoust. 1969, AU-17, 93 B~ith,M. "Spectral Analysis in Geophysics', Elsevier. Amsterdam, 1974 Strong motion earthquake aceelerograms digitized and plotted data, from Calif. Inst. Technol. Earthquake Eng. Res. Lab., Vol. I., Parts A-F, Pasadena, California, 1970 Caughey, T. K. J. Appl. Mech., 1960, 27, 269 Archer, J. S. J. Struct. Div. A.S.C.E., 1963. 89, 161 Rayleigh, J. W. S. "The Theory of Sound', Vol. I., Dover, New York. 1945 Morduchow, M. Trans., A.S.M.E., 1961, 28, 458

Eng. Struct., 1978, Vol 1, October

17

18

Collings, A. G. and Tee, G. J.. "An analysis of Euler and Implicit Runge-Kutta numerical integration schemes for structural dynamic problems', Sixth Australasian conj', mech. struct, mat.. Christchurch, New Zealand, August, 1977 Cooley, J. W. et al. IEEE Trans. Audio Electroacoust. 1967. AU-15, 79

Appendix A L G O L 60 procedure for summing trigonometric series Parameters. Procedure Goertzel (n, a, b, sine, Jbrward) Input parameters: n is the upper limit of the cosine series, with n - 1 as the upper limit of the sine series. a is the array of coefficients of the sum. sine is a Boolean variable, when it is true the sine sums are evaluated, otherwise cosine sums are evaluated. forward is a Boolean variable, when it is true a forward transformation is evaluated and otherwise an inverse transformation is made. Output parameter: b is an array which will contain the interpolation coefficients when the forward transform is computed or the results of the summation of the series when the inverse transform is performed. T H E A L G O L 60 procedure procedure Goertzel (n, a, b, sine, forward): value n, sine, forward: integer n: real array a, b: Boolean sine, forward: begin integer i, il, i2, j: real gamma, pi: array u[0:n] : pi: = 4 x arctan(l); il: = if sine then l else 0; i2; = n - il; for i : = i l step l until i2 do begin gamma: = 2 x cos(i x pi/n): u[0]:= 0: u[1] : = if sine then 0 else if forward then a[n]/2 else a[n] : for j : = 2 step 1 until n do u[j]: = gamma

x u[j -

l ] -

u ~ -

2]

+ a [ n - - j + 1]: bill := if sine then sin(i x pi/n) x u[n] else if forward then gamma/2 x u[n] - u[n - 1] + a[0]/2 else gamma~2 x u[n] - u[n - 1] + a[0] : if forward then b[i]: = 2 x b[i]/n end i: if forward ^ --1 sine then begin b[0] : = b[0]/2: bin]:= bin]~2 end: comment extreme coefficients for forward cosine are given half weight end Goertzel.