The solution of structural dynamics problems by the generalized euler method

Cmpufm

d Smrwu,

Vd. IO,pp.50545.

Papmoa

PIUS 1979. Prided in Great Main

THE SOLUTION OF STRUCTURAL DYNAMICS PROBLEMS BY THE GENERALIZED EULER METHOD ANDREWG. CouINGst and GARRY J. TEES The Universityof Auckland,New Zealand (Receiued2.5May 1978) Mawt-StilI systems of secondarder ordinaryditlerentialequations,which describe the vibrationof structures subject to dynamic load& (for example, by earthquakes)may be solved by a variety of numericalalgorithms.A one-pammetergeneraliza& of E&r’s classical s&me for ordinarydiierential equations is investigated,and is shown to be applicableto such problems. The value of the pammeterand the size of the time&p can be chosen on the basis of detailedanalyses madeof the scheme. The Euler scheme is compared with Newmark’s “beta” method. The application to non-linear probkms is Micated. As an example a 3-stony space frame, in the form of a triangularprismwhich is subjectedto dynamicloading, is computed by several versions of the Eukr scheme, and the results are compared with the exact solution. Compkte ALGOL 60 procedures are included, for applying the generalized Euler scheme to linear structural dynamicsproblems.

operating on the velocity vector. In such cases, numerical integration schemes for solving the equations of motion have often proved successful. They are executed in a step-by-step fashion using increments of time, h, ti+C%+KX=F(t), (1) whereby the displacements, velocities, and accelerations are evaluated at the end of each step, accompanied by appropriate adjustments, (if necessary), to the connecting the displacement X(I) of a multidegree-offreedom structural system with an applied load P(t), can coefficients of the huations of equilibrium. be solved by a variety of analytical or numerical techThe dynamic analysis of a multi-degree-of-freedom niques, the choice of which will depend upon. the exact system using numerical integration techniques generally form of the mathematical model used to represent the involves the simultaneous direct integration of a set (1) of second-order differential equations, one for each physical problem. In eqn (1), M, C and K are respectively the general@@ mass, damping and stif+s matrices of ,degree of freedom. The numerical work involved will order w X, X and X are m-vectors of generalized invariably be substantial, and hence the efficiency of any accelerations, velocities and displacements in an inertial computational algorithm used will be sign&ant. The frame of reference, and p(t) is an m-vector of time- step-by-step numerical integration techniques used will dependent forcing functions. Since this equation of be subject to numerical errors, involving numerical inmotion is linear the exact analytical solution can be stability, truncation error, spurious damping, etc. In constructed, but it is often rather too complicated for general, eqn (1) represents a stiff system. A system is said to be stiff if there is a wide spread for the periods of practical computation [11. its constituent modes. The accurate integration of such a The mathematical model, which is defined by eqn (l), system may require a short time step so that the contriresults from discretizati?n of the continuum representing the structure-it has a spectrum with a minimum bution of the short-period components to the overall frequency corresponding to the fundamental mode, and a response could be adequately assessed. However, the maximum cut-off frequency which is a characteristic of use of such small time steps would invariably lead to the discretized mode1 rather than of the original excessive computation and accumulation of round-off continuum. Forced vibration of the system above this errors. Alternatively, if the time step was chosen solely cut-off frequency cannot be represented adequately by to be some fraction of the largest period in the system, then numerical integration schemes involving explicit the discretized model. The system may also have varying degrees of non- formulations would be unstable, whilst implicit schemes, linear behaviour-for example, the mass, damping and although having better stability characteristics, would fail in varying degrees to assess the contribution that the stiffness matrices may have coefficients which depend upon the stress history, or the motions of the system short-period components in the system make to the may well not be restricted to small displacements from overall response. In this paper there is investigated a fa&ly of one-point the position of static equilibrium. In addition, the damping might be non4inear. such that it cannot be numerical integration schemes, based on Euler’s classical method[2]. The Euler schemes have a number of comrepresented, even to a llrst approximation, by a matrix mendable features in addition to theii stability, not the least of which is the fact that they are simpler to impletCivil bginehg Department. SCompWationalMathematicsUnit, apartment of Matbfzma- ment than many other methods. In practice the Euler schemes have often been passed over in favour of more tics. A system of m linear secondarder ordinary differential equations:

505

506

A. G. COLLINGS and G. J. TEE

complicated methods, without full consideration being given to their advantageous features. However, some problems which have been found intractable by sophisticated numerical integration algorithms have been solved successfully by Euler schemes[3]. Several features of the Euler schemes are demonstrated in this paper, by applying them to a structural dynamical problem. The computed results are compared with the exact solutions to the governing differential equations. In addition, it is shown that the Euler schemes have similar form to the popular Newmark “beta” method-indeed, the Euler and Newmark schemes coincide for special values of their parameters [4]. Tlatmum-

In the genera&d Euler methods, the divided difference of any differentiable function U(t) over an interval (L+ 1.) is approximated by a weighted mean of the values of ir at t = t,_, and f = t,: (U(f,)-

U&-d/h

-wlj(t,)+(l

-o)ir(L,),

(2)

where h = f,, - fn_,, and o is a real parameter (0~ o 5 1). Hereafter, in this paper, a symbol such as xn will denote the discretized approximation corresponding to X at t = t,, where X(r) is the solution of the original O.D.E. An initial-value problem for a set of M-order O.D.E‘s. t = f(f, X), X(&J= x0

(3)

where both X and f are m-vectors, may be solved by the generalized Euler method, applying the approximation (2) to each component of (3)[5]. We have t, = f.-1+ h, = Of(f”, x.1 + (1 - o)fL,

(x. -x,-,)/h

X.-l).

(4)

The Explicit method, as devised by Euler[2], uses o = 0, for which X~is readily obtained as: X” = X”_l tf(t.4,

x,-,)h.

(5)

With o > 0 the eqn (4) has x,, on both sides, and hence it must be solved for x. in some manner, for instance the Newton-Raphson method (43-54). The (wholly) implicit method uses o = 1, xn = xn--I + f(L, xn)h,

(6)

the 2m vector Y as: Y=

[‘yl, I 2

where Y, = X and Y1= X; then (8) may be written as: Y = F(t, Y),

(10)

where

The Euler schemes are single-point methods, and hence are simpler to implement than multi-point methods, which require information at more than one previous point to evaluate the dependent variables at successive times. Although multipoint methods may have smaller truncation errors than single-point methods they are prone to instability; but there are some which are fairly stable and give greater accuracy than do single-point methods]& 91. (Reference [91 also contains a clear exposition of some numerical methods used in structural dynamics.) However, multi-point methods do pose additional complications, for at least two reasons. The first is that they require a separate starting procedure to compute a set of initial values, and this complicates substantially the entire computational algorithm. Secondly, difficulties may arise when dealing with systems in which the coefficients of the mass, damping and stiffness matrices change during a numerical integration time step. The Euler methods may conveniently be applied directly to a linear equation, such as (l), but for a non-linear equation of the form (8), the computation of xn involves the use of some predictor-corrector technique in which an initial approximation to the dependent variables at t. is to be computed by a suitable “predictor formula”. This initial approximation is then refined by a suitable “corrector” formula[8] (also, see 43-53). Certain limited categories of non-linear behaviour in dynamical systems may be computed by appropriate modifications of the linear second-order system (1). For example, a so-called “degrading system”, in which the coefficients of the mass, damping and stiffness matrices change after high stresses occur in the system, may be approximated in a step-by-step analysis by a sequence of linearized equations, each applied over an interval, during which the coefficients of the matrices are assumed to remain constant or to vary linearly with time.

and the Balanced Euler method uses w = l/2. X” = xn-1

+{f(fn-,, x,-d + to,, x.)W.

(7)

This Balanced Euler method is also known as the TrapeZoidal Rule, and when it is applied to heat-conduction problems (with spatial discretization) then it is called the Crank-Nicolson method [6,7]. The generalized Euler method may be applied to O.D.E’s which are not initially written in the canonical form (3). by converting them to that form. For example, a set of m second-order O.D.E’s,

FoilwRATNmoFmlI&R8cBpDlespoR JJNEM vmRATIoNEQuATIoNs The set of m linear equations (1) could be reduced (as in 11) to a set of 2m first-order linear O.D.E’s‘ i = GZtH(r),

(12)

where:

z= [;I,

G=

[&

_;-‘cl, E= [&?(J (13)

t = g(r, x; X),

(8)

may be represented as a first-order system thus: Define

Denote the velocity vector X as S, and the acceleration vector X as A.The application of the generalized Euler

The solution of structuraldynamicsproblemsby the generalikedEuler method method (4) to this system may he represented as follows:

(zn -zn-d =G(wz, +(1 - o)z,...,)+ h ._ P~tioni~

the equation:

Mao=Puo)-Cso-Kxo,

wH(t,)

+ (I- W)~(t,,_,).

(14)

507

cm

for which purpose M should be factorized as: M = IJ=u,

this eq~tion into:

(28)

by Cholesky fac~~tion. Similarly, for each u and ii, Cholesky factorization should be used to decompose:

we get that: Xn-x,-t -= h

(16)

OS”t (I- 0)8”_,,

from the first partition, and: sn -&I-1 -=M-‘(o(p(f”)-cs”-wjr,)+(l-oxP(r”-,) h -

CL1

- Kx.-4,

(17)

from the second partition. But, (1) may be rewritten as: MAtCStKx=P(t),

(18)

so that: A = M-@(t) - CS -KK),

(1%

and accordingly, (17) may be rewritten as: sn - - sn-1 = 08. -t (1 - ~)a,-,. h

(20)

sn = sn+ + hjoa, t (I- ~)a,&

(21)

x,, = xn-1 t h(os, + (1 - o)s.-,I.

(22)

Therefore,

and

Substituting (21) into (22), this reduces to: x, = x~_.~+ hs,_, + o’h*a, + ~(1 - w)h*a,-,.

(23)

It follows, from (18), that: KK=P(t)-MA-CS.

(24)

- C(s.-, + Ohan t (I- o)ha.-d. (25)

Therefore, we get an equation for the new

-K(x,-,

t hsn-, t &(l

acc&ratiOn:

+ h(l - ~)a,-,) -oh+,).

=

aan t (1 -a&l,

x,-x,-i p=s,,-,+h h

=Kx,=P(L)-Ma,-cc

fM+ ohC f o*h*K)a, = P, -C&-Z

CoMPARlWN wrrn ~~~ The widely used “beta” method of Newmark may be represented by the following formulae corresponding to (20) and (16) (see [4]). h

K(x,+ + hs,-r + o*h*a, + ~(1 - w)h2a,-,)

(29)

so that (26) may be readily solved for an (see Appendix). Note that all matrices used are of order m, and not 2m. The matrices M, C and K are real symmetric, with M being positive-definite and K being positive semi-definite (or positive-definite). It is readily shown that if C is also positive-semi-definite (or positive-definite), then Q is positivedefinite for all non-negative o and h, and hence Cholesky factorization can always be applied. However, if C is not positive semi-definite (as can happen in Rayleigh’s model of damp~g), then Q will be positive definite (and hence factorizable by Cholesky) only for sufficiently small values of oh. After an has been computed by solving (26), s,, and x,, may be computed from (21) and (22) respectively. Note, that in this application of the Euler method, if the matrices M, C and K have band form with all non-zero elements less than b rows or columns from the diagonal, then the Cholesky factorization and the back substitutions using U and V may be performed within the halfbands of width b. Thus, the computational cost of each of the two Cholesky factorization is approximately (l/2)&’ multiplications and additions, and thereafter the cost of solving for a,, involves about 2mb multipli~tions and additions. Furthermore, the evaluation of the righthand side in eqn (26) (to be solved for an) costs about 4mb multiplications and additions, plus the cost of evaluating P(L) (see Appendix). Note that a,+ in (26) is not to be re-evaluated.

Bn -%I-1 -

Pre-multiplication of (23) by K gives:

= P&)-Ma,,

Q = M t ohC t &*K = V=V,

(26)

The initial acceleration & may be computed by solving

Here a and /3 are parameters, where usually 0 rs;B 5 l/4 and (I = l/2. Note that (30) corresponds exactly to (21) witb a = or and that when a = 112 and #3= l/4 (the “unsent acceleration” method) then Newmark% equations (30) and (31) reduce to the Trapezoidal Rule, as does the Euler method with o = l/2. The application of the Newmark method to (1) may be implemented in a manner similar to that used for the Euler method. Corresponding to (2% we obtain an

508

A. G.

COLLINGS

and G. J. TEE

equation to be solved for a.: d. o. f.2

(M+ahC+@h*K)a,, =I’&)-(C(s.-,+h(l -K(x.-,

+hsn-, + (i-g)h*a.-,).

-ala.-,)

t/1

1‘

(32)

IOM

After a. has been constructed, then sn can be computed from (30) and then x, can be computed from (31). In fact, the procedure Eufetgeneralized in the Appendix would require but minor modification in order to implement the general Newmark method. .9rAMLrrrANJJC~E Detailed analyses of the stability and the convergence of the generalized Euler method are given in two other papers[l, IO]. The method is shown to be stable for all time steps if l/2<@ 51, with 0 = 1 giving the strongest stability; but for o = l/2 with a linear system (12) of O.D.E’s the method will be stable only for values of /I small enough for that circle through the origin whose centre is ll[h(l - 2011 to contain all the eigenvalues of G. Thus the time step must be a small fraction of the shortest natural period of the system, when 0 I o IS;l/2. The stability analysis (which concentrates on linear O.D.E’s)[l] shows that o = l/2 is the minimum o for which the Euler method is stable for all time steps. Thus, we should not be unduly surprised to find some form of numerical instability manifesting itself when the Balanced Euler method is applied to non-linear probkms [a]. Accordiiy, in such non-linear problems, it may be advantageous to use a more strongly stable form of Euler’s method, with o > l/2. The paper on stability includeb detailed analyses (for the case of linear O.D.I?s), of period perturbation, damping of components, the “ringi&’ of high-frequency components, etc. Also, it is shown that for sufIiciently “smooth” systems of O.D.E’s (not necessarily linear) the results of computations us& two diKerent step sizes can be extrapolated, to produce a refined estimate of the solution. That extrapolation is just&d by the convergence analysis, which[lO] shows that for sufliciently “smooth” systems of O.D.E’s (non-linear or linear) the truncation error at a tixed time is asymptotically proportional (as h-+0) to h, or to h* if o = l/2. Nonetheless, the Euler method with o > l/2 does have some properties which for some problems, make it preferable to use o > l/2 rather than 0 = l/211]. MHIBMATICAL m m mR TmJVrNAMtcANAL= The dynamic investigation was con&d to the linear elastic analysis of. a family of three-storey triangular steel-membered space frames, all having the same basic overall dimensions (Fig. l), but with different crosssectional properties for their constituent members. An equation similar in form to that of (1) represents the time-dependent response of such systems. In eqn (1) K and M, which are the generalized overall stiflness and mass matrices of the system, are constructed from the individual stiRtess and mass matrices ,of the elements into which the structure has been dissected for the purpose of analysis. The stiffness and mass matrices of the elements are given as follows: K, =

B=BB du,

(33)

I

d.af.X

1; IOM

+ IOM

Fig. 1. Steel-memberedspace frame. The system comprises 54 d.0.f.d per node. and M,=

I

P

N=N dp,

(34)

where the suflix j sign&s the jth element, du is an element of volume and dp is an element of mass. The matrices B, E and N have the properties:

l=Bd, U=K4$

(351

f=Nd, where u and E are vectors of stress and strain for the element, f is a vector which represents the displacement of an arbitrary point within the element in terms of a matrix of shape functions N, and d is a vector of generalized nodal coordinates. The shape function N used for the formulation of the mass matrices was chosen to be identical with that used for the element sti@ness matrices, and corresponds to situ& static displacement. Mass and stiKness matrices forntukted in such a way are calkd consistent with each other-they are symmetric and have respectively the proper& of positive-definiteness and positive semi-dsliniteness[ll]. Dissipative forces proportional to the rehative velocities of the system were presumed to exist, and hence the damping was represented by a matrix C operating on the velocity vector. A Baykigh mode1 was used, in which it was assumed that: C=aM+/3K

(36)

where a and fl are constants. In this particular case a and /3 are given by a model proposed by Morduchow,

The solution of structuraldynamicsproblemsby the generalizedEuler method

and are as fo~ows[lZ~:

509

so that

Ilk(t)= - Xk(t) f

A* e-‘fl cos Izt,

(39)

and R(r) = - 2X&) - (1+ a%%. where Ti and 4 are the periods of the ith and jth modes, (i < ,ii which have ratios of critical damping associated with them of Ai and Ai. (Taken as 0.02 and 0.05 in the dynamic analyses.) This simplified damping model of Bayleigh has the feature that the resulting damped vibration modes are the same as in the undamped case, but the eigenvahtes do depend upon the Bayleigh[l3] ~oeffi~nts a and 8. However, nurne~~ ~t~t~n solution techniques do not use any eigenvector analysis, and it is therefore likely that more realistic damping models could profitably be employed, e.g. the mode superposition method in which a viscous damping ratio is specified separately for each vibrational mode [ 141. ~~~~~~ In order to construct specific problems of the form (l), vectors X were chosen as arbitrary functions of time (but with some relationship to physical plausibility), and then the vector P was constructed as: PQ)=&kX+CX+XX, The computed function P(t) was then used, together with the corresponding initial values of X and X, as a test problem for various algorithms for numerically solving the eqn W-these algorithms should give solutions approximating to the function X(t) chosen originally. In particular, solutions were chosen of the form, Xk(t) = Ar e-’ sin at,

(R = 1,. . ., m),

(38)

(40)

For any t, P(t) is evahmte$ by ~0~s~~~ X(r), X(r), X(r) and then computing hfX(f) + CX(t) + XX(?). The corresponding initial values are: xk(o)

= 0,

(k = 1,. . ., m),

(41)

(k=I

(42)

and ~~(O)=Ak~,

,.a., m).

AMLYfas OF THlE coMPuTATloNs

Computations were made with two sets of crosssectional and elastic properties for the members of the space frame shown in Fig. 1. In Model A the natural periods (of the undamped system) were found to range from 0.099 see to 0.21 set, and in Model B the range was from 0.005 see to 2.2 sec. In the imposed soiution (35), G was chosen as 40 in Model A, so as to give a period of 0.157~ for the (damped) forcing function (38) and in Model B it was 1.57, giving a period of 4 sec. The values used for the Bayleigh parameters (34) were: (a,@) = (1.0591,0.00014)in Model A, and (0,0.00155) in Model B. Figure 2 displays the actual error of the results computed by Balanced Euler (o = 112)and by fully Implicit Euler (o = 1) for degreesf-freedom 20, (which is the vertical displacement at node 4 of the frame) for Model A, plotted at tiie intervals of At = 0.0125set, with Euler time-steps of h = At/S, At/l0 and Ad20. It is evident that, at each fixed time f, the error is fairly closely pro~~onal to h (when o = 1) or to h2 (when 0 = l/2).

_. -

. excct displocsmentof D.O.F. 2D

-

differeftce between exact soluticn and E&r rokrtharkx varkus schemes

Iw=QS - Euler l \_ step w*

Euler impkit scheme

balanced scheme kame as Newmarki 8=1/4)

h - time

Njwo625 scale x IO5

W

.

h.

-71

1 0.0125

i 0025

8 a0375

1 0.06

i CC625

1 0.075

t 0.0675

, 0.1

Time, seo Fii. 2. Displacementof degree-of-freedom20 with errorof Euler solutions, for model A.

I 0.1125

I 0.125

510

A. G. COLLINOS and G. J. TEE

Figure 3 displays similar results for the acceleration of degree-of-freedom 1 (which is the horizontal displacement of node 1 towardS node 3 in the frame), for Model A. The time steps used (with o = l/2 and o = 1) were h = At, h = Ad2, h = At/S and h = At/lo. Once again, it is evident that, at any fixed time t, the errors are closely proportional to h (for 0 = 1) or h2 (for 0 = l/2)-note, in particular, the errors for o = l/2 with h = 0.0125 and h = 0.00125. We observe that in Fig. 2 the errors are oscillatory, whereas in Fig. 3 they increase monotonically.

and a. is the discretized approximation to %(t,). This equation is to be solved for a,,, with s,, and x. being expressed in terms of a,, by (21) and (23) respectively. (The initial acceleration ao is to be found by solving 46, with n = 0.) The Newton-Raphson method may be applied for solving (46). Starting from an initial estimate a(O)for a. (e.g. ato)= a.-,), a sequence of estimates {aCk’)(k > 0) is constructed for a,, in the following manner. The corresponding estimates for sn and xn are defined as: S(')= sn_, t

A general set of m non-linear second-order O.D.E’s may be represented as in (8). In non-linear vibrational problems it is usually possible to define an inertia matrix M (except where geometric non-linearities dominate), such that:

X

[El=

(w

L&)l~

(1 - o)hs.-l + ohs”‘,

(48)

(*I = x+,

+ hs._, + ~(1 -

o)h2a,-, t 02h2a’“‘. (49)

(43)

where d is a vector of generalized forces. This may be converted (as in 9-11) to a set of 2m non-linear firstorder O.D.E’s;

ii

+

(47)

so that: X

Mii = aqt, x, ri),

(k)= X”_,

(1 - u)ha._, t da’“‘,

The Newton-Raphson equation: m’k’

= ,qtn,

x(k-‘),

S(k-0)

estimate a(‘) is defined by the _ jQk’(x’k)

_ x(k-O)

_ c’k’(s’k’

_ s(k-l)),

(50)

*

where the matrices Kck’and C”’ are defined as: in which we write S for 2. The application of the generalized Euler method to (44) bkcomes: ~(x”-x~-‘)=WIM4nt(1-0)1M8.-l,

6. -S.-d h

@I(t,

(51)

X, S),

I

[K’*%=2 @‘it&Xv S),

h

M

[Pk’]ij = 2

I

=

W~(tm,X.,S.)+(1--~f”-‘,X”--l,Sn--l), (45)

and thus eqns (16) and (20)-(23) apply, where now: Ma.

=

wt.,

X.,

Sn),

for i and j = 1 to m, where the partial derivatives of Cp are to be evaluated at (t., x(k-‘), s(*-‘?. Substituting (47) and (49) into (50), that equation simplifies to become: Ma(k)= @,((t x(I’-I),&I’-0) _ n,

(46)

_

02h21(‘k)(a’k’

_a(k--I))

o@k’(a’k’

_ a(k-‘)).

Time, set

Fig. 3. Accelerationof degree-of-freedom1 with errorof Eulersolutions,for model A.

(53)

The solution of structural dynamics probkms by the 8eneralized Euler method

Thus, we get a set of m iinear equations to be solved for p). (M + ohc”’ + &ZKU0)a(k) = ,I,(~ x’k-0, $k-1) “, ) + oh(C"' t whK(k’)a(k-”

(54)

We observe that the matrix operating on a(*’ is closely similar to the matrix Q operating on an in the linear case (see 26). Thus, when one time-step of the generalized Euler method is to be applied to a non-linear set of O.D.E’s (43), each iteration of the Newton-Raphson method (for soiving the non-linear equations 46, 21 and 23) involves the solution of a set of m linear algebraic equations (54), whose matrix resembles closely that matrix Q in the set of m linear algebraic equations which need to be solved for performing one time-step of the generalized Euler method for a linear set of O.D.E’s (1). Indeed, the procedure Eulergeneralized(see Appendix) could be modified slightly for use within a procedure to apply the Newton-~phson method to non-iin~ O.D.E’s. That procedure would need to supply the modified version of Eulergeneralizedwith the matrices Ctk’ and K’“’ of partial derivatives of 8. The computation of those derivatives is usually very costly, and hence in such multi-dimensional Newton-Baphson iterations it is customary[lSl to use a mod&d version of the al~~thm, in which the matrices of partial derivatives are used for several iterations, rather than being recomputed for each iteration. Whereas the NewtonRaphson algorithm gives asymptotically quadratic convergence to the solution (if the O.D.E’s are sufficiently “smooth”, and if the initial estimate is sufhciently close to the solution), the modified algorithm (with constant matrices of partial de~vatives) initially converges (if at all) approximately quadratically, but ultimately converges linearly. That linearity of the convergence may be detected by the fact that the increment vectors x(‘)x(k-I) and +kl _s(k--I) then tend to decay by a constant factor after each iteration. When the convergence is thus observed to be linear, it would then be appropriate to re-compute the matrices of partial derivatives of CBat the latest point. The matrix in (54) would then need to be re-factorized. When ack)has effectively converged to an, then sck)and x(‘) will have converged to the Euler solutions s. and xn (see 47 and 48). Note that the physical connections within the vibrating structure will be reflected in the resulting spar&y of the matrices C”’ and K”’ (see 51 and 52). as in the case of linear elasticity and damping, where those matrices are constant. Moreover, as in the linear case, judicious ordering of the elements of X wilI frequently enable M, C’*’ and K’*’ to be represented as narrow band matrices which enables considerable economies to be made in the solution of the linear algebraic eqn (54), as with linear O.D.E’s. These features justify the use of O.D.E’s in the form of (44) rather than (S), since the corresponding treatment of (8) would usually result in matrices corresponding to Cck’and K”’ which were full, and which would be much more cumbersome to handle.

The generalized Euler method has several desirable which include stability, simplicity, ease of ap plication, and susceptibility to detailed analyses. Most other numerical integration schemes are more difficult to C.&s Vd. 16. No. 3-F features,

511

anaiyse rigorously, especially in the case of non-linear equations. For example, the Newmark “beta” method has proved to be quite difficult to analyse[4]. It was published in 1959, but general stability and convergence analyses were not published until the 1970’sin papers by NickelHi, 17). Goudreau and Taylor 1181, et al.; with an analysis of the Newmark constant acceleration method (or Balanced Euler) for the integration of non-linear structural dynamic equations in a paper [6] by Hughes. The stability analysis[l] of the Euler method shows that it is unconditionally stable, if the Euler parameter w is within the range l/2 I o 11. The stability of the Euler scheme is very strong with o = 1 (the Fully Implicit Mew), and it weakens when w is reduced towards l/2 (the Balanced Euler Method). For values of w less than l/2 the Euler scheme is stable only for time steps less than some limit, depending upon the damping and the periods of the system: indeed, the step length must be a small fraction of the shortest natural period of the system. Such req~ments clearly lead to hopelessly uneco~mic ~rnputi~ times, and it is for this reason that most attention has been directed to implicit Euler methods. This analysis indicates that the Generalized Euler Method should be considered seriously, as a practicable algorithm for computing structural vibrations of linear and mildly non-linear systems. REFERENCEi

I. A. G. Callings and G. J. Tee, Stability and accuracy of the

2. 3. 4. 5. 6.

7.

generalized Euler method for ordinary dilferential equations, with reference to structural dynamics problems. Engng Srccturrs (to appear in 1979). L. Euler, l~~~~~~jo~~sCalculi Integralis, Volumen Primum. Petropoh (1768). H. H. Robertson, The so&on of a set of reaction rate equations. In Numerical Analysis: An Introduction (Edited by J. Walsh), p. 180.Academic Press, London (1966). g. M. Newmark, A method of computation for structural dvnamics. J. Enrmn Mech. Diu. A.S.C. E. 85.67-93 (1959). J.- D. Lambeg, -Comptcrntiottol Methods In ‘Ot&ary ~~e~~~f ~ua#~~, p. 240. Wiley, London (1972). T. f. R. Hughes, Stability, convergence and growth and decay of energy of the average acceleration method in nonlinear structural dynamics. Comput d Structures 6,313324 (1976). J. Crank and P. Nicolson, A practicalmethod for numerical evaluation of solutions of partial differential equations of heat-conduction type. Proc. Cumb. Phil. Sot. 43, SO-67

(1947). 8. A. Ralston, A First Course in Nume~caf Analysts. Chap. 5. McGraw-Hill, New York (1964). 9. 1. L. Humar and E. W. Wright, Numerical methods in structural dynamics. Can. 1. Ciu. Engng 1, 179-193(1974). 10. G. I. Tee, Convergence of the generalized Euler method for ordinary diRerential equations, (to appear). Il. J. S. Archer, Consistent mass matrix for dis~but~ mass systems. S. Struct.LXv.A.SC.E. 89, 161-177(1963). 12. M. Morduchow, On classical normal modes of a damped linear system. Trans. A&ME. 28,458 (1961). 13. T. K. Cam&y, Classical normal modes in damped linear dynamic sy&&. 1. Appl. Mech. 21, 269-271(1%6). 14. R. W. Clouah and S. Moitahedi, Barthquake response analysis consid&ing non-proportional damping. EaLqu. Engig Strwt. &n. 4,4&4% (1976). IS. J. C. Butcher, On the implemen~~n of implicit RungeKutta methods. BIT 16,237-240(1976). 16. R. C. Nickel, On the stability of approximation operators in problems of structural dynamics. fnt. L Solids Strvctwes 7, 301-319(1971). 17. R. E. Nickel, Direct integration methods in structural dynamics.3. Engng Mech. Diu. A.S.C.E. 99,303-317 (1973).

A. G. COLLINGS and G. J. TEE

512

CIIOLBSYY PACTCNIUATION ANDBACK-ION A sym tric positive-definite matrix A may be factorized as A.=VTV,$ here V is an upper triangular matrix with real elements, and with positive diagonal elements. Moreover, if A has half-bandwidth b then so does V (see Fig. 4). The algorithm is based on the following considerations. We require that:

18. G. L. Goudreau and R. L. Taylor, Evaluation of numerical integration methods in elastodynamics. Computer Met/r. Appl. Mech. Engng 2,6!Ml(l972).

ALGOL 60 procedures for applying the generalized Euler method to problems in stmctuml

dynamics

i-l

“f =

vi1 =

a,

a2

a,.

.

(lb+1

ab+*.

.

03

ab+2

a2b+t

.

G=

.

a2b+*

rc-I

hiuki w

Vii

Since ail = 0 for j 2 it b, it follows by induction from (60) that vi,= 0 for j P it b, and hence ah non-zero elements of V occur within the upper half-band. Thus, we may reduce the previous expressions to become:

[(a.-$V&)“’

RI =I

ifj=i

.

.

.

.

.a26

.

.

:

] )

i--l

air - x WVki _ k-4

j=itl,...,p

Vii

a2

s,v:i* i-l

aij -

Thus,

i=l

aii -

and when i < j,

The array a is caged the one-dimensional representation of G. If G id symmetric then gii = gii, and if i 5 j then g,, is found from g,,, which is represented within a according to the foregoing rule.

(57)

. ., m).

Thus, setting i = j, we get that:

(55)

(ilj
= 1,

In view of the symmetry of A, it is sufficient to consider only i 5 j. But, ski= 0 for k > i (since V is an upper triangle). and hence:

ONK.DIMD#IONAL RBpILpsENTAInmop MATIucEs In these procedures, all matrices (whether symmetric or triangular) are represented in the following manner. A matrix G of order m, with half-bandwith b, is represented by a real array a[l: m*b], in which the elements on and above the diagm are stored thus: &=a[(i-l)*btj-itl],

(v’hk(v)kj= $,okiiukj,(i, j

ait = 8,

All matrices used in these procedures are either symmetric (air = ori for all i and j), or else upper triangular (i > j j all = 0); and hence for both types of matrix it is sufficient to store only the elements on and above the diinal (i S j). Also, the mass, damping and stiffness matrices will usually be sparse (i.e. most elements are zero), and by appropriate ordering of the rows and columns it will usually be found possible to get all non-zero elements close to the diagonal. A matrix A is said to have a half-bandwidth of b when (i - jl2 b j ali = 0. Thus, b = 1 for a diagonal matrix, and b = m for a matrix of order m which has no zero elements. Ah of the matrix operations performed here preserve the zero elements outside the bands, and considerable economies in computational time and in storage space may be gained by avoiding the zero elements outside the bands. Accordingly, only the elements within the bands need be stored.

(i=l,...,m),

(61)

1

ab

. ojb

hb+l

;b

(56) a26

0

i=m_ j=l

j=m 1

The elements of a corresponding to the extension of the upper half-band beyond column m of G are not used. If the matrices being used cannot conveniently be represented in band form, with b/m * 1, then this onedimensional representation can still be used, e.g. with b = m. However, the procedures would operate more efficiently if they were mod&l in obvious ways, so as to operate on matrices which are represented in the standard manner as doubly-subscripted real arrays. Each representation of a square matrix of order m requires a working space of exactly m* real variables.

--

Fig. 4. Cholesky factorization of band matrix.

where p=min(m,i+b-l),

q=max(l,j-btl).

(62)

These equations may be used to evaluate the elements of V provided that a suitable ordering is used, e.g. vi,, viz,. . ., ulbr v22. v23..

. . . v2.b+lr

013,

UY, . . . .

Each us uses only ail and elements of V which precede vi, in this ordering, and hence C may be overwritten by the elements of V as they are computed in this order. The computation of each non-trivial v,, requires one division (or one root extraction, if i= j). In column i of V the computation of rj_b+rj. er_b+sj, . ., ujj (insofar as these exist) requires res$ectiveiy 0, 1,2,. . ., b- 1 muhiphcations and subtractions. Thus, the computation of column j of V requires (1/2)b(b - 1) mubiplka~n~ and subtractions, or fewer when C< b. Accordingly, the entire Cboksky factorizdon of A as V V requires m root extractions, less than (1/2)mb(b - 1) multiplications and

513

The solution of structural dynamics probkms by the generalized Euler method divisions, and less than m(b - 1) divisions. (In contrast, Cholesky factorization of a full non-banded matrix requirea about (1/6)m3 multiplications and subtractions, about (1/4)mZdivisions and m root extractions.) Once A has been factorized, any equation Ax= II,

(63)

may be solved for the m-vector I thus: h=Ax=V’Vx=V’z,

(64)

where vx=s.

(65)

Conversely, if z and x can be found such that V’z = It and Vx = a, then Ax=V~Vx=v~Z=b.

(66)

The equations for E and then for x may readily be solved, since V is triangular. Indeed, for i = 1,. . ., m;

(67) so that z,, z2,. . ., z,,, may be computed successively by the formula

j is an integer expression, whose value is assumed to be intherangesi-mandi-i+b-1. Result: AA supplies the value of A[i, j) from all, which is the one-dimensional representation of A. (5) pmcedurr factor (a, b, m). Input Parameters: a is the one-dimensional representation of a symmetric positive-definite matrix A, b is the half-bandwidth of A, m is the order of A. Output Parameter: a becomes the one-dimensional representation of the upper triangular matrix V, which is the Cholesky factor of A = VTV. (6) real procedure V(i,j). Similarly to AA (No. 4) it gives the value of V[i, j] from

its one-dimensional representation a. (7) real procedure AA(i, j). Similarly to AA (No. 4) it gives the value of A[i, j] from its one-dimensional representation a. (8) procedure backsub (a, r, b, m). Input Pammeters: a is the one-dimensional representation of an upper triangular matrix V, r is a vector h, b is the half-bandwidth of V, m is the order of V, and of h. Output Parameters: r becomes A-‘h, where A = V’V. (9) real procedure V(i, j).

Similar to No. 6. (IO) procedure Eulergeneralized (n, b, m, c, k, w, h, x, dx, ddx, t,

(68)

DS, first).

input Parameters :

n is the order of the matrices hi. C. K and of the vectors with q = max (i - b + 1,l). Similarly, for i = 1,. . ., m.

b

m c with p = min (i t b - 1, m). Therefore, x,, *,,,-I.. . ., XI may be computed successively as:

k

w (70) Thus, the computation of each element of x and x requires one division and at most b - 1 multiplications and subtractions, so that the entire computation of x for any h requires 2m divisions and less than 2m(b - 1) multiplications and subtractions. PAUMEWJmOPTEEpIoclpwREs (1) integer procedure min (x, y). Input Parameters: x and y are arithmetic expressions. Result: min gives the value of the smaller of x and y (rounded to integer, if real). (2) integer procedure max (x, y). Input Parameters: x and y are arithmetic expressions. Result: max gives the value of the huger of x and y (rounded to integer, if real). (3) ptucedurc allxd(al1, d, y. m, b). Input Parameters: a 11 is the one-dimensional representation of a symmetric matrix A, d is a vector, m is the order of A, and of d, b is the h&bandwidth of A. Output Pammeter: y: = Ad. (4) real pmcedure AA(i, i). Input Pammeters:

i is an integer expression, whose value is assumed to be intherange l-m,

h

x dx t ps

x. & dd.x, is the half-bandwidth of all the matrix parameters, is the one-dimensional representation of a symmetric positive-definite matrix M (the “inertia” matrix), is the one-dimensional representation of a symmetric matrix C (the “damping” matrix), is the onedimensional representation of a symmetric positive-definite (or semi-definite) matrix K (the “stitfness” matrix), is the real parameter used in the generalized Euler method, is the time-step used in the generalized Euler method, is the initial n-vector of displacements, at time t, is the initial n-vector of velocities, at time t, is the time at the start of an Euler step, is a procedure, which must have been declared with the procedure heading: procedure ps@n, n, t); e *n;&t;realarraypn;

n, t;

The procedure body (to be supplied for any specific problem), gives the n-vector pn the value of P(t), first is a Boolean variable which is to be assigned the value E before the first execution of Eulergeneralized, but for each subsequent execution it carries forward the value falser ddx is not reauired for the first execution when ihst is &, but for later executions it carries forward the old n-vector of accelerations, for the old t. Output Parameters: x is the new n-vector of displacements for the new time

t, after one application of the generalized Euler method to the system of n second-order O.D.E’s:

dx is the corresponding new n-vector of velocities, a’dx is the corresponding new n-vector of accelerations,

514

A. G. COLWNGSand G. J. TEE

t is the new time, equal to the old value plus b. first has the value false. (The procedure uses internally the own real array u[l: n X b]. During the first execution it is initially used to contain the

+

procedure min (x, y): + min: -then -- x else y

X,y; e

one-dimensional representation of the upper triangular matrix U, where UrU = M; then it is used similarly for V, where VTV= M + whC + w%*K, and this is carried forward internally for use in subsequent executions, when first is false.)

x, y ;

integer procedure max (x, y); value x. y; integer x, y; max: = -d x > y -then x else 7 procedure a 1ixd(a 11,d, y, m, b); value m, b; integerm,b;reaJarrayall,d,y;commenty:=Ad; begin integer KS 1, ~2; real yi ; real procedure AA(i, j):& i, j; integer i, j; comment i 5 j; -?tA:=all[(i-l)xb+j-itl] for ii = 1step 1 until m do he&t~=0~:=~x(i-btl,l);s2:=min(itb-l,m);

end i -et&f allxd procedure factor (a, b, m); value m, b; integerm,b;realarraya; comment A = VrV, where V overwrites A, in a ; beginintegeri,j,k,p,q;reaJe; real procedure V(i, j);-?& i, j; w i. j; i?-=aI(i-l)xbtj-itl] real procedure AA(i, j); -value i, j; integer i, j; AA:)

q9 1until i - I do e: = P - V(k, i) x V(k,j); ;;iii - 1)x b t j -ml: = $T= j then sqti(e) -else e/V(i, i)

endj end

.-

e&f

factor

procedure backsub (a, r, b, m); & b, m ; lntegerb,m;realarraya,r; comment Solves Ax = II, where A has heen factorized as VTV, and V is stored in a. x overwrites b in I: heginintegeri,j,p,q;reale; real procedure V(i,j);alue i, j; integer i, j; V:=a[(i-l)xbtj-i+l] for i: = 1 *

iG+

1 until m do

q: = max (i-b+

13); e: = r[i];

comment Solve V’z = h, where z overwrites h in r; forq*lt@i-l&e:=e-V(j,i)xr[i]; q]: = e/ V(i, i) end i; comment Now, solve Vx = z, where x overwrites z in r; fori:=--l_lldo &p:=mm(itb-l,mT,e:=r[i]; forj:=i+lstepluotilp&e:=e-V(i,j)xrU]; %I: = e/ V(i, i) end i

-en&f backsub procedure Eulergeneralized (n, b, m, c, k, w, h, x, dr, ddx, t, ps, first); wn,b,h,w; integer n, b; & w,h, t; Boolean first; procedure ps; real, array m, c, k, x, dx, ddx; hegm integer i; real wh, wl, wh 1, wh2, dxold; -

The solution of structural dynamics problems by the generalized Euler method own real array u[l: n x b]; -real array pn, q,stl, st2[1: n]; -w/r:=wxh;wl:=l-w;whl:=wlxh;wh2:=w/rlxh; jf tirst then wimment Factorize M as UTU and store in u; for w=Istep 1 until n x b do v[i]: = m[i]; factor (u, b, n); comment Initialii~fir), a&hen find initial acceleration by solving Mddx = P- Cdx - Kx; bs(pn,t); a 1lxd(c, dr, stl, n, b); a 1lxd(k,x, ~2, n, b); for i: = 1* 1 until n do ddx[i]: = pn[i] - stl[i] - sf2[i]; backsub (u, ddx,G); comment Factorize M t whCt w*h*Kas V’V; for i: = I 9 1 until n x b do ofi]: = m[i] t wh x (c[i] t wh x k[i]); factor (u, b, n); fii= falseend of initialization; comment Evaluate P(r), for the new t; t: St h; ps(pn, n, I); comment Perform one step of Euler; fori:=ls&3IuntiJ~n& be& stl(i): = dx[r] t whl X ddx[i]; st2[i]: = x[i] t h x dx[i] t wh2 x ddx[i] a lm(r. sf 1, q. n, b); a 1Ixd(k, st2, stl, n, b); for i: = 1 * 1 until n do q[i]: = pn[i] - q[i] - stl[i]; backsub (0, q, b,n);om%ent q is the new acceleration:

for i: = 1 * 1 untiii beg@ dxold: =dr[i]; dx[i]: = dxold t wh x q[i] t whl x ddx[i]; x[i]: = x[i] t whXdx[i] + wh 1 X dxold; ddx[i]: = q[i] end i -end ofulergeneralized

515