Some characteristics of frequency-dependent structural matrices

004s7949/90

Computers & Sfrucrures Vol. 35, No. 4, pp. 413-416. 1990

$3.00 + 0.00

Pergamon Press plc

Printed in Great Britain.

SOME CHARACTERISTICS OF FREQUENCY-DEPENDENT STRUCTURAL MATRICES W. D.

FILKEY

and N. J. FERGUSWN

Department of Mechanical and Aerospace Engineering, The University of Virginia, Charlottesville, VA 22901, U.S.A. Abstract-The traditional finite element method applied to dynamic problems employs shape functions which are based on a static displacement assumption. The more exact approach uses frequency-dependent shape functions and frequency-dependent mass and stiffness matrices. These matrices are usually constructed as a series expansion in powers of the frequency of vibration. A formula is presented which relates the higher ordered terms in the mass matrix expansion to those in the stiffness matrix expansion. An additional result gives expressions for the mass and stiffness matrices in terms of the values of the shape function at the boundary of the element. The results presented hold for models based on shape functions derived as solutions to differential equations which are the Euler equations for the particular variational principle used to define the mass and stiffness matrices.

INTRODUCTION

FORMULATION

matrices used in traditional finite element analysis of the dynamic effects of vibrating structures are in actuality the first terms in the Maclauren expansion of the exact mass and stiffness matrices, which depend explicitly on the frequency of vibration. For the freely vibrating structure the assumption of harmonic motion allows the displacement field u(x, t) to be interpolated between the nodal displacements U(t) by the relationship

The free oscillations of an undamped governed by the familiar expression

The mass and stiffness

4x, f) = W, w)U(t),

Ku+Mii=O,

are

(2)

where M and K are the global mass and stiffness matrices assembled from the individual m and k matrices for each element, and u is the vector describing the displacement field. The assumption of harmonic motion allows one to write the displacement vector as

(1)

where o is the frequency of vibration [l]. Such frequency-dependent shape functions result in mass and stiffness matrices which also depend on w. If enough elements are employed the frequencyindependent shape functions N(x) and the resulting frequency-independent mass and stiffness matrices provide a satisfactory approximation. These matrices, however, are simply the first terms in the power series expansion of the full blown frequency-dependent matrices. For a more exact approach, additional dynamic correction terms can be retained in the expansions, resulting in a nonlinear eigenvalue problem. In such a case, there is evidence that better accuracy is afforded for comparable costs [l-9]. A formula is introduced here which relates the higher ordered terms in the mass matrix to those in the stiffness matrix. In addition a corollary is presented allowing both the mass and stiffness matrices to be obtained from the dynamic stiffness matrix, which can be computed solely from boundary information. The results are shown to hold true for structural members such as beams, circular plates, and shells of revolution whose shape functions are obtained as solutions to the governing differential equations for the element.

structure

u(x, t) = u(x)eiw’.

This results in a shape function form

(3)

expression of the

u(x) = N(x, w)U,

(4)

where U is the amplitude of U(t). The shape function N is then inserted into the formulae

m=

s

N’pN dI’

(5)

Y

k=

N’,DTED,N

d I’

(6)

for the element mass and stiffness matrices, where p is the generalized density matrix [lo] (usually scalar), E is the elasticity matrix, and D, is the differential operator appearing in the strain-displacement relation 6 =

413

D,u.

(7)

W. D. F’ILKEYand N. J.

414

The symbol “Dr represents the transpose of D, and acts on the matrix which preceeds it. The element matrices are combined to form global matrices K and M. These are then used to form a nonlinear eigenvalue problem [K(w) - o*M(w)]U = 0.

(8)

The matrix D = K - w*M is often called the global dynamic stiffness matrix, the symbol d being reserved for the element dynamic stiffness matrix. The boundary conditions applicable to the structure can be enforced, and then eqn (8) can be solved for the natural frequencies using, for example, the iterative scheme [ 111 [K(w, ) - w ?+I M(~i)lUt + I = 0.

from which the shape functions N,, N,, N,, . . are calculated. Appropriate boundary conditions for eqns (15) and (16) are supplied by the element nodal conditions on the displacements, the components of N and their derivatives assuming the values of zero or one accordingly at the various nodes. For the higher order coefficients N, , N,, N, , . . only zero boundary conditions occur. The Ni, having been computed, can be used to obtain the mass and stiffness matrices, the appropriate formulae being

(9)

Even for simple elements the transcendental expressions for m and k are complicated and the power series expansions K m= 1 mnW2” n=O

(10)

k = f k,,W2n “=O

(11)

urUDrsdI’

N,?“DrED,N,_,dI’.

(18)

u’A*sdS

=

-

urD,TsdI’

(19)

for any vectors u and s smooth enough for the integrals to converge, where AT is a matrix of direction cosines assembled from the components of the outward unit normal to the boundary aV, and is defined by the relationship

(12)

N,w2”

can be used for the shape function. Actually, expansion can be taken in the variable w rather than CB*but it will be found that the odd power terms all vanish. Such power series converge uniformly within a certain radius about o = 0, whereas the series diverges for larger frequencies, but still may serve as an approximation when truncated. The local governing equation for the element has the displacement form DfED,u - piI = 0,

(13)

where the differential operator Df, which appears in the stress form of the equilibrium equations, Drs-piI=O,

(14)

describes the divergence of the stress tensor, s being a column of generalized stresses or stress-resultants, including moments. Insertion of eqns (3), (4) and (12) into eqn (13) yields the recursive system of differential equations

DfED,N,+pN,_,=O,

(17)

The main result of this article relies on a certain formal adjoint relationship satisfied by the operators D,’ and D,

p=ATs

“=O

DiED,N,

N,TpN,_,dI’

s ”

will be used. A similar expansion N= f

FERGUSON

= 0

(15) n >I,

(16)

(20)

for the stress vector on an oblique plane. The above adjoint relationship is merely a statement of the Gauss integral theorem for tensors, the stepping stone from equilibrium to the principle of virtual work. In the case of a beam, for example, this equation reduces to the simple statement which describes integration by parts. Inherent in this discussion is the assumption that the local governing differential equation used to derive the shape functions is the Euler equation corresponding to the principle of virtual work used to derive the form of the mass and stiffness matrices. This is in contrast to the methods used by some authors [1] where, for example, the rotary inertia of a beam is taken into account by including extra terms in the mass matrix of the form

s I

NT r;pN dx,

0

but terms involving the radius of gyration ry are not included in the shape functions. It should be noted that although the results of this paper depend on this assumption, the additional accuracy obtained by including rotary inertia terms in the shape functions is not appreciable, except for higher frequencies.

415

Frequency-dependent structural matrices As an example consider a pinned-pinned beam whose radius of gyration is one-tenth the length of the beam. Let 1, and 1r be the frequency parameters corresponding to the lowest two natural frequencies. Define A2= w~~,,/(~A/EZ), where p is the mass per unit length of beam, 1 is its length, A is the cross-sectional area, and EI is the flexural rigidity. The exact values of 1 for such a rotary inertia beam are 1, = 3.0690 and 1, = 5.7817. Consider a one-element model with one extra dynamic correction term in the dynamic stiffness matrix. Using a simple beam model which does not include rotary inertia in the shape functions, but does include it in the mass matrix, yields 1, = 3.0945 and 1, = 6.0028. A more precise model which includes rotary inertia in the shape functions produces 1, = 3.0937 and 1, = 5.9800. For such a model a simple relationship between the terms m,, m,, m,, . . and k,, k,, k,, . . . , namely (n+l)k,+,=nm,,

n>O,

(21)

will now be demonstrated. For the case n = 0, the 1.h.s. of eqn (21) simplifies to k,. Equation (18) gives k, =

(N;,DrED,N, IY

+ NT,DrED,N,)

dV,

Using eqn (15) and the fact that N, vanishes on 8 V, this last expression becomes N;pN, dV+

NfpN,dV, s

sV

V

which is simply m, , thus the formula holds for n = 1. The remaining cases, n > 2, can be proved in a similar fashion. Another interesting formula which relates the mass matrix to the dynamic stiffness matrix [lo] states that

ad

m=-a(o2).

(22)

Since power series converge uniformly the expansions of eqns (17) and (18) can be differentiated term-byterm and formula (21) can be derived using eqn (22). In fact, for low frequencies (within the radius of convergence) the two formulae are equivalent and so the restrictive assumption that the shape function must solve a differential equation which is the Euler equation for the virtual work principle should also apply to eqn (22). A useful corollary is obtained by first eliminating the higher order stiffness terms in favor of the higher order mass terms, resulting in

which by eqn (19) becomes d=k,-&co*-

f

im,_,co2n

(23)

n=2

k, =

N;,D’EAN,

N,T,D’EPrN,

dS -

dV

s V

I 8V N;ArED,N,

+

and similarly, eliminating NfD,TED,N,, dV.

dS -

I PV

2k, = 2

(N;,DrED,N,

These two formulae allow the mass and stiffness matrices to be obtained from the dynamic stiffness matrix, which can be computed solely from boundary information as follows:

dV.

=

N’,D’ED,N

dV - co2

sV The first and last terms vanish for the same reasons that k, vanishes. The remaining term is split into two equal halves 2

Nf,DrED,N,

(24)

d=k-w2m

+ N[,DrED,N, + N;,DrED,N,)

1 f knco2”. n=2n - 1

d=k,-m,w2-

s V

Equation (15), in conjunction with the fact that N, satisfies homogeneous boundary conditions, implies that this expression is zero. For the case n = 1 consider the 1.h.s. of eqn (21),

the mass terms,

=

N’pN I

NrArED,N

dS -

dV

V

N’DTED,N dV s V

-CO2 N’pN dV

dV

sV =

Nf,D’ED,N,

dV +

N;D:ED,N,

s V

=

s N;,DrEAN,

4 JV

+

dV

V

dS -

N;,DrE,DrN,

dV

=

NrArED ”NdS

.

Thus

s V

N[A’ED,N,

dS -

N:D’ED,N, s V

dV.

d=

NrArED YNdS .

(25)

W. D. F?LKEYand N. J. FERGUSSON

416

This result is merely a way of stating formally that the dynamic stiffness matrix coefficients are precisely the nodal forces due to unit displacements at the various nodes. Similar boundary integrals for the terms m,, ml, m2,... and k,, k,, k,, . . . are easily obtained using eqns (23) and (24) by introducing eqn (12) into eqn (25). This results in

k0 =

N,A’ED,N,

dS

(26)

k, =O k,=

(27)

-(n

NTATED @N n_i-dS j=o i PV J

(28)

II+1 N,rArED,N, + , _j dS. C j=O i8V

(29)

- I) i

m,= -++l)

A word of caution arises when using eqn (25) in conjunction with power series. Consider the third term of the dynamic stiffness matrix which can be computed three different ways using the above results

(30), (31) and (32) depends on the accuracy desired and the tools available, SUMMARY

For the frequency-dependent mass and stiffness matrices used in dynamic analyses, the higher order terms in the mass matrix expansion are shown to be proportional to the terms in the expansion of the stiffness matrix. The constant of proportionality depends only on the order of the term being considered. Also shown is the fact that either the mass or stiffness terms can be calculated from the values of the shape functions at the nodes. These results are applicable to the class of elements whose shape functions are equilibrating in the sense that they solve an equilibrium equation consistent with the structural matrices derived from them. In particular the integrals defining the mass and stiffness matrices must follow from a variational principle whose Euler equation is precisely the differential equation used to define the shape functions. REFERENCES

(N;ATED,N,

d2 =

+ N;A’ED,N,)

+ N;ArED,N,

dS

(30)

or

dz= -;

(N;pN,

s 1’

+N:pN,)dJ

(31)

or d, = -

~N~~DrED~N* + N~*D~ED~N, sY

+ N;,DrED,N,)

dV.

(32)

The boundary method of eqn (30) has the advantage of requiring only boundary information but, on the other hand, requires the knowledge of the higher order shape function Nz not found in eqn (31). Equation (32) requires the derivative of N, throughout the whole body. When using numerical techniques this derivative may be obtained with higher accuracy than the function itself since the function itself requires two integrations whereas the derivative would require only one. The choice between eqns

I. J. S. Przemieniecki, Theory of Matrix Structural Analysis. McGraw-Hill, New York (1968). 2. J. S. Przemieniecki, Quadratic matrix equations for determining vibration modes and frequencies of continuous elastic systems. Proc. Conf: Matrix Meth. Struct. Me&., Wright-Patterson Air Force Base, OH; AFFDL TR 66-80 (1966). 3. A. J. Fricker, A new approach to the dynamic analysis of structures using fixed frequency dynamic stiffness matrices. Inr. J. Numer. Meth. Engng 19, I 11l-1129 (1983). 4. K. K. Gupta, Solution of quadratic matrix equations for free vibration analysis of structures. Int. J. il’umer. Meth. Engng 6, 129-135 (1973). 5. K. K. Gupta, On a finite dynamic element method for free vibration analysis of structures. Comput. Merh. a&. Mech. Engng 9, 105-120 (1976). 6. K. K. Gupta, Finite dynamic element fo~ulation for a plane triangular element. Inf. J. Numer. Mefh. .&ngng 14, 1431-1448 (1979). 7. K. K. Gupta, Development of a finite dynamic element for free vibration analysis of two-dimensional structures. ini. J. Numer. Meth. Engng 12, 1311-1327 (1978). matrices for 8. A. K. Gupta, Frequency-dependent tapered beams. J. Srrucf. Engng 112, 855103 (1986). 9. B. Downs, Vibration analysis of continuous systems by dynamic discretization. J. Meeh. Des., ASME 102, 391-398 (l980). 10. 0. C. Zienkiewicz, The Finite Element Method, 3rd edn. McGraw-Hill, Maidenhead (1977). 11. T. H. Richards and Y. T. Leung, An accurate method in structural vibration analysis. J. Sound Vibr. 55, 363-376 (1977).