The dynamic element method for analysis of frame and cable type structures

Engineering Structures 27 (2005) 1906–1915 www.elsevier.com/locate/engstruct

The dynamic element method for analysis of frame and cable type structures Anders Ansell∗ Department of Civil and Architectural Engineering, Royal Institute of Technology (KTH), SE-100 44 Stockholm, Sweden Received 25 February 2005; received in revised form 17 June 2005; accepted 20 June 2005 Available online 1 August 2005

Abstract With the dynamic element method (DEM), results more accurate than with the conventional finite element method (FEM) are obtained with the same number of degrees of freedom. This is due to the introduction of shape functions of polynomial type, introducing frequency dependence into the mass matrix expressions. It is demonstrated how this affects free vibration analysis, including the solution of nonlinear eigenvalue problems. Various numerical techniques for solving these polynomial problems are discussed. The polynomial matrix formulations for stiffness and mass matrices are given for a beam, a bar and a cable element. Numerical examples demonstrate how the DEM can be implemented for modelling of frame type structures and its efficiency is compared to that of the conventional FEM. © 2005 Elsevier Ltd. All rights reserved. Keywords: Dynamic element method; Free vibration; Nonlinear eigenvalue problem; Frame; Cable

1. Introduction The finite element method (FEM) is today the commonly used tool for static and dynamic analysis of structures. Large or complex structures are often described by large matrices, thus giving rise to a large number of operations to be performed within an analysis. Recent computer development, together with software that facilitates numerical programming, provides large computer capacity at a reasonable price and therefore gives new possibilities to perform effective and accurate FEM analyses. The implementation of the FEM consists of discretising a structure into an appropriate number of segments, elements, interconnected at nodes with a finite number of degrees of freedom describing possible motions of the structure. Results of higher accuracy are often obtained by performing further discretisation of the model, i.e. adding degrees of freedom. An alternative method in a dynamic analysis is to use more exact mathematical descriptions of the elements ∗ Tel.: +46 8 790 80 41; fax: +46 8 21 69 49.

E-mail address: [email protected]. 0141-0296/$ - see front matter © 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.engstruct.2005.06.010

in the model. Conventional FEM formulations are based on shape functions of polynomial type, approximating the real behaviour of the elements. If shape functions of more accurate forms are used, frequency dependence will be introduced into the matrix expressions. An important part of a structural analysis is the determination of free vibration properties, i.e. the natural frequencies of a structure, often called resonance frequencies, and the corresponding free vibration modes. The free vibration case is obtained when no external loads are applied, here given by Kd (ω)v = 0

(1)

which is a nonlinear eigenvalue problem. The solutions are the natural frequencies ωr with corresponding free vibration mode shapes vr = v(ωr ). The dynamic stiffness matrix Kd (ω) is assembled from element contributions Ked (ω), incorporating the stiffness and inertia properties of the structure. For the conventional FEM, the approximated shape functions lead to Kd (ω) = K0 − ω2 M1

(2)

A. Ansell / Engineering Structures 27 (2005) 1906–1915

where M1 is a consistent mass matrix, composed from element contributions Me1 , and K0 the static stiffness matrix. In this case, the approximate expression of Eq. (2) will give a linear matrix eigenvalue problem that can be solved by well known numerical procedures. The approximated displacement fields involve discretisation errors, giving exact accuracy only as the number of global degrees of freedom approaches infinity. If a more exact form is chosen, the eigenvalue problem will be of nonlinear type due to the frequency dependence of the included shape functions. The nonlinear problem has to be solved by more sophisticated numerical procedures than for linear problems. The dynamic stiffness method (DSM), often referred to as the exact method, makes use of exact shape functions for the formulation of Kd (ω). The method is thoroughly reviewed by Fergusson and Pilkey [1–3] and also described by Banerjee [4] and Banerjee and Sobey [5]. There are also methods based on approximations to the DSM, giving results still more accurate than those from a conventional finite element model with the same number of degrees of freedom. This paper demonstrates free vibration analysis using a method based on polynomial matrix approximations. 2. The dynamic element method In a finite element formulation, the dynamic displacement field u(x, t) is interpolated according to
u(x, t) = Nd (x, ω)v(t) (3) v(t) = v0 eiωt where x is a coordinate in the element, t represents time and v(t) is a column vector of nodal displacements. The static and dynamic shape functions N(x) and Nd (x, ω) are further defined as solutions to the boundary value problems given by DσT EDu N = 0

(4)

and DσT EDu Nd + ρω2 Nd = 0

(6)

and DσT σ − ρ v¨ = 0

be formulated in terms of the boundary values of the shape functions as  NT DuT EDu NdV (8) Ke0 = V

and



Me (ω) =

NT ρNd dV.

(9)

V

Replacing the dynamic shape function in Eq. (9) by a truncated Taylor series in even powers of the frequency (the odd ones vanish) Nd =

q 

Nd,i ω2i

(10)

i=1

will give a polynomial including the global Mi , corresponding to the elementary Mei . Substituting M1 in Eq. (2) results in Kd (ω) ≈ K0 − ω2 (M1 + ω2 M2 + ω4 M3 + · · · + ω2q−2 Mq ).

(11)

This formulation is, when truncated to more than one Mq matrix (q > 1), often referred to as the dynamic element method (DEM) [2]. The method is also reviewed by Mukherjee and Chattopadhyay [6] and Ansell [7]. Note that a truncation to i = 1 will result in the conventional FEM formulation of Eq. (2), making the conventional FEM a special case of the DEM. This is also the least accurate approximation of DEM type to the exact DSM formulation. 3. Polynomial eigenvalue problems The most common type of solution methods for polynomial eigenvalue problems is based on conversion to linear problems, i.e. a problem of size n and degree q is transferred to a linear problem of size nq. The quadratic problem, Eqs. (1) and (11) with q = 2, is often reformulated to the linear form using reformulations similar to (B0 + λB1 )y = 0

(5)

subjected to boundary conditions dictated by the form of v0 where ρ is density and E an elasticity matrix. The differential operator matrices, Dσ and Du , satisfy the conditions ε = Du v

1907

(7)

where ε and σ are the vectors of generalised strain and stress, respectively. In Eq. (4), N(x) give displacement fields of static type as used in the conventional FEM to approximate dynamic displacement fields. The static stiffness matrix and mass matrix of the element can then

here with   −K0 0 B0 = , 0 M2

(12)

B1 =

  M1 M2 M2 0

 and y =

v λv



which has the advantage of preserving the symmetry of the original problem. The solutions are the eigenvalues λ = ω2 with the corresponding eigenvectors v, i.e. the free vibration mode shapes. Peters and Wilkinson [15] and Ruhe [16] present a method for reformulating higher order polynomials. The method consists of reformulating the original problem of dimension n in polynomial form of degree q into a linear eigenvalue problem of dimension nq as     v v =λ (13) B va va

1908

A. Ansell / Engineering Structures 27 (2005) 1906–1915

where 

0 0   B=0  ..  .

I 0 0 .. .

0 I 0 .. .

··· ··· ··· .. .

0 0 0 .. .

      

B0 B1 B2 · · · Bq−1

with B0 = −M−1 q K0 ,

Bi = −M−1 q Mi

and I denoting the identity matrix. The solution of Eq. (13) will also include, apart from the n relevant λ that are solutions to the physical problem, (q − 1)n eigenvalues that are only solutions to Eq. (13). These irrelevant eigenvalues are in all realistic cases easily recognised as they are negative or complex [7]. The eigenvector v can be separated from the augmented part va as consisting of the n first terms of the eigenvector (vT , vaT )T . One drawback with this method, apart from the large dimension of B, is the need for evaluation of the inverse of Mq . However, it has the advantage that the well known Q R algorithm [17] can be used and Peters and Wilkinson [15] state that this method has proven to be one of the fastest and most accurate when all solutions are wanted and no approximate eigenvalues or vectors are known. Another type of methods uses improvement techniques based on solving a sequence of linear eigenvalue problems. A method of this type for the structural dynamic problem in DEM form is suggested by Sadek [10]. First the original linear eigenvalue problem of the conventional FEM formulation is solved, followed by an improvement step including one or more dynamic correction terms. Alternative corrections are suggested and the method can be described as; step 1, solve: (K0 + λM1 )v = 0

(14)

step 2, solve one of: [K0 + λ(M1 + λ0 M2 )]v = 0 [K0 + λ(M1 + λ0 M2 + λ20 M3 )]v = 0

(15) (16)

where λ0 is the old eigenvalue and λ is the improved eigenvalue. Two linear problems must be solved for each approximated eigenvalue and the method is therefore computationally expensive. The advantage is that it can, at least in principle, easily be extended to be valid also for long polynomials and it may also be used in an iterative manner, converging towards the correct λ. 4. Mass and stiffness matrices The polynomial matrix formulations for use with the DEM are here given for three types of element; a beam, a bar and a cable. The beam element is intended for frame structures while the bar is suitable for trusses, where all elements are pin-connected at the nodes. Each scalar element m i j in the frequency dependent element mass

Fig. 1. A 3D beam element with 12 degrees of freedom.

matrix Me (ω) is obtained as a function of the corresponding scalar element ki j in the static stiffness element matrix Ke0 , through multiplication by a dynamic function fi j = f i j (ω). The functions fi j (ω) are in the following presented as fourterm Taylor series, derived from the exact expressions given in the Appendix. 4.1. A beam element The two-node three-dimensional (3D) Euler–Bernoulli beam element in Fig. 1 is assumed to be straight and uniform, i.e. there are no changes of area A or moments of inertia I along its length L. The element is homogeneous, consisting of linear elastic material of density ρ, Young’s modulus of elasticity E and shear modulus G. The moments of inertia of the cross-section about the y and z axes are I y and Iz , respectively, the polar moment of inertia is I0 and the torsional constant It . The element is defined in a 3D x yzsystem with a full set of six degrees of freedom at each node, describing translations and rotations with respect to the three axes. Comparison to published results show that the matrices presented in the following are identical to Ked (ω) given by e.g. Williams and Wittrick [8] and to Me (ω), as given by Melosh and Smith [9]. Sadek [10] gives Ked (ω) and its four-term Taylor expansion including Ke0 , also identical to the matrices presented in the following. The static stiffness matrix for the beam element is   T K11 K21 (17) Ke0 = K21 K22 where the sub-matrices K11 , K21 and K22 are  (−1)i+ j ka 0 0 i+ j  0 (−1) k z1 0  i+ j  k y1 0 0 (−1) Ki j =   0 0 0   0 0 (−1) j k y2 1+ j 0 (−1) k z2 0  0 0 0 0 0 (−1)1+i k z2    0 (−1)i k y2 0   0 0 (−1)i+ j kt   0 k y(3+i− j ) 0 0 0 k z(3+i− j )

(18)

A. Ansell / Engineering Structures 27 (2005) 1906–1915

including the stiffness functions EA , L

ka =

kt =

G It L

{ k y1 k y2 k y3 k y4 } =

and E Iy { 12 6L 4L 2 2L 2 }. L3

The functions k z1 , k z2 , k z3 and k z4 are obtained by substitution of I y with Iz in the above expressions. The frequency dependent mass matrix is expressed as a function of the static stiffness matrix by matrix scalar multiplication (term by term)   T F1 (ω) ∗ K11 F2 (ω) ∗ K21 M (ω) = F2 (ω) ∗ K21 F1 (ω) ∗ K22 e

(19)

where ∗ denotes component-wise multiplication. Here F1 (ω) and F2 (ω) are:     Fi (ω) =    



f ia 0 0 0 0 0 0 0 0 f (i+2)z  0 f iz  0 0 f iy 0 f (i+2)y 0   (20) 0 0 0 f it 0 0   0  0 0 f (i+2)y 0 f (i+4)y 0 f (i+2)z 0 0 0 f (i+4)z

The dynamic functions are calculated as f ia = fiat (αa ) f j y = f j (β y )

f it = f iat (αt )

i = 1, 2.

f j z = f j (βz )

j = 1, 2, 3, 4, 5, 6.

with the parameters β y4 =

Aρ 4 2 L ω , E Iy

αa2

ρ L2 2 ω = E

βz4 = and

Aρ 4 2 L ω , E Iz αt2

ρ L 2 I0 2 = ω . G It

The truncated Taylor polynomials are



 α2 f1at (α) 840 = f2at (α) 2520ω2 −420
 α4 3104640 + 139708800ω2 −2716560
 α6 9686476800 + 4576860288000ω2 −9383774400
 α8 34691380486963200 + 163916772800901120000ω2 −34420354076908800 (21) + ···

1909

      f1 (β) 78  4248             f (β)        −27 −3837          2    4 8 β β 22 1784 f3 (β) + = 2 2 −13 −1681 f (β)    2520ω   4   139708800ω               6  568  f (β)          5    −9 −1097 f6 (β)    264480       −261045     12 β 113504 + −112631 4576860288000ω2       36576      −72819   18812077440       −18779817285      16 β 8093086080 + −8084877645  163916772800901120000ω2       2611657088      −5220181117 + ··· (22)

4.2. A bar element The Euler–Bernoulli bar element shown in Fig. 2 is a version of the well known simple extensible 3D-bar, in this case also allowing bending. The extension bar is usually described by 2 × 2 matrices, but the introduction of an additional four degrees of freedom gives matrices of size 6 × 6. Compared to the beam element of Section 4.1, which it is theoretically related to, the bar element has no degrees of freedom for torsional or flexural vibration. Otherwise, all features and assumptions are common with the beam element and some dynamic functions for the bar element are identical to those forming the matrices in Section 4.1. All matrices have proven to be identical to those presented in [10]. The static stiffness matrix is   ka 0 0 −ka 0 0  0 0 0 0 0 0    0 0 0 0 0 0  (23) Ke0 =  −ka 0 0 ka 0 0    0 0 0 0 0 0 0 0 0 0 0 0 using the same stiffness function ka as in Eq. (18). The quasistatic mass matrix is expressed as   F1 (ω) ∗ K −F2 (ω) ∗ K (24) Me (ω) = −F2 (ω) ∗ K F1 (ω) ∗ K where the sub-matrices K, F1 (ω) and F2 (ω) are   ka 0 0 K =  0 k z1 0  and 0 0 k y1   0 0 f ia 0 . Fi (ω) =  0 f (i+6)z 0 0 f (i+6)y In the sub-matrices, the stiffness functions k y1 and k z1 are identical to those in Eq. (18). The dynamic functions, unique

1910

A. Ansell / Engineering Structures 27 (2005) 1906–1915

and a static stiffness matrix Ke0 . Here, the latter was formed using the first, frequency independent terms in a Taylor series expansion of Ked (ω), depending on the horizontal component of the cable tension H through the relation Tϕ = H / cos ϕ. This represents the static cable tension where the cable is parallel to the chord and corresponds approximately to the average cable tension. This gives 

Fig. 2. A 3D bar element with six degrees of freedom.

k 0 0 0 1 0  Tϕ  0 0 1 Ke0 = l  −k 0 0  0 −1 0 0 0 −1 where k = c1 = Fig. 3. A 3D cable element with six degrees of freedom. The angle ϕ is always related to the plane orthogonal to the direction of gravity.

for the bar element, are f j y = f j (β y ) and f j z = f j (βz ) with j = 7, 8. The truncated Taylor polynomials, derived from the exact functions given in the Appendix, are


  β4 β8 f 7 (β) 70 24640 = + f 8 (β) 2520ω2 35 139708800ω2 23870
 12 β 8153600 + 4576860288000ω2 8137675
 β 16 2995025510400 + 163916772800901120000ω2 2994659906700 + ···. (25) 4.3. A cable element A dynamic stiffness matrix for an extensible, flexible, sagging cable suited for modelling cable-supported structures such as bridges and guyed masts is presented by Starossek [11,12]. This work is also referred to by Kim and Chang [13] and Au et al. [14] who further discuss dynamic stiffness formulations for cable structures. The cable element in Fig. 3 is considered as a continuum and its unloaded geometrical shape is given by the chord length l and the sag d perpendicular to the chord. The cable is assumed to be uniform with cross-section area A and of linear elastic material of density ρ and with Young’s modulus of elasticity E. The element is derived in a local x yz-coordinate system where the x-axis is parallel to the chord. Due to the vertical direction of the gravity g, the angle of inclination ϕ of the chord versus the horizon must be known in the local system. The dynamic element matrix Ked (ω) has in the following been rewritten as a frequency dependent mass matrix Me (ω)

12c22 (12+c22 )c12

 −k 0 0 0 −1 0   0 0 −1  k 0 0  0 1 0 0 0 1

(26)

and the cable parameters

E Al ρ Agl d c22 = c12 cos ϕ = 8 , Tϕ l Tϕ L     2 2 c d L ≈l 1+8 2 =l 1+ 1 . 8 l

and

The corresponding mass matrix is 

m 1 m 2 0 −m 1 m 2  m 2 m 3 0 −m 2 m 4  Tϕ  0 m5 0 0 e  0 M (ω) = −m −m 0 m −m l  1 2 1 2   m 2 m 4 0 −m 2 m 3 0 0 m6 0 0

 0 0  m 6 . 0  0 m5

(27)

The element functions for  a truncated Taylor series can, by introduction of Ω = ωl ρ A/Tϕ , be written as   6c22   c22 Ω 2 m 1 (Ω ) 5 = 2 m 2 (Ω ) c (12 + c22 )  c1 (12 + c22 )2 ω2     − 1  2  c2 (1020 + c22 )     2  c22 Ω 4 700 + 2 3c (12 + c22 )  c1 (12 + c22 )3 ω2    − 1  5   2 2 4   c2 (111600 + 240c2 + c2 )       c22 Ω 6 63000 + 2  c1 (1020 + c2 )(12 + c2 )   c1 (12 + c22 )4 ω2   2 2  −  1400 c22 Ω 8 + 2 c1 (12 + c22 )5 ω2   2 c (417916800 + 1400400c22 + 9120c24 + 37c26)       2 194040000 ×   c1 (111600 + 240c22 + c24 )(12 + c22 )     − 126000 + ··· (28)




   

A. Ansell / Engineering Structures 27 (2005) 1906–1915

1911

solved by using linear reformulation, Eq. (13). In addition, the improving methods of Eqs. (14)–(16) were also tested. The numerical methods were implemented using the Matlab numeric software [18], as far as possible using the built-in functions in this package. Table 1 Natural frequencies of the four-storey frame Freq. no.

Fig. 4. The four-storey plane frame.






  Ω2 2 Ω4 m 3 (Ω ) 192 + c22 = + 2 m 4 (Ω ) 6ω2 1 720(12 + c22 )ω2 168 − c2
 Ω6 9216 + 24c22 + c24 + 2 4 2 2 2 30240(12 + c2 ) ω 8928 − 24c2 − c2
 Ω8 442368 + 864c22 + 36c24 + c26 + 438912 − 36c24 − c26 1209600(12 + c22 )3 ω2 + ··· (29) 



Ω2 2 Ω4 Ω6 8 32 m 5 (Ω ) = + + 2 2 m 6 (Ω ) 6ω 1 360ω 7 15120ω2 31
 Ω8 128 + + ···. (30) 604800ω2 127

From numerical investigations, Starossek [12] states that Eqs. (26)–(30) are valid for inclined cables only if the limiting conditions  2 c2 ≤ 24 or ϕ ≤ 60◦  c1 ≤ 0.10 (i.e. d/l ≤ 1/80)  2 c2 ≤ 24 ϕ ≤ 30◦  c1 ≤ 0.24 (i.e. d/l ≤ 1/33) are satisfied. Otherwise, the presented matrices are only applicable for horizontal cables, i.e. with ϕ = 0◦ . 5. Examples Three examples are given, aiming at demonstrating possible improvements in the resulting frequencies when compared to results from conventional FEM calculations. It is not demonstrated how the resulting natural frequencies and vibration modes are used in e.g. a mode superposition analysis aiming at finding the structural response to dynamic loading. The studied ranges of frequencies were chosen to demonstrate the improvements possible for frequencies corresponding to vibration modes of increasing number. The polynomial forms of Eq. (11) were used with 1–5 mass matrices and the nonlinear eigenvalue problems were

1 2 3 4 5 6 7 8 9 10 Relative CPU-time

No. of M-matricesa q =1 q =2 q =3

q =4

q =5

Exactb

12.87 40.69 72.38 102.64 210.11 243.52 251.87 295.28 314.61 318.47

12.87 40.63 72.16 102.37 191.67 217.51 223.58 255.31 268.48 271.15

12.87 40.63 72.15 102.36 188.97 213.07 218.61 247.23 258.74 261.08

12.87 40.63 72.15 102.36 188.47 212.06 217.44 244.96 255.84 258.04

12.87 40.63 72.15 102.36 188.37 211.82 217.14 244.24 254.86 257.00

12.87 40.63 72.15 102.36 188.34 211.73 217.03 243.89 254.31 256.41

1.0

2.1

6.3

13.5

22.9

(>300)

In rad/s, computed with different M(ω). a Mass matrix defined as: M(ω) = q ω2(i−1) M . i i=1 b Results from DSM analysis [7].

5.1. A four-storey plane frame In this example, a four-storey frame was modelled in 2D, using one beam element per structural member. All the members are identical and the geometry of the frame is shown in Fig. 4. The length of each element is 2.0 m and the cross-sections are 50 × 100 mm2, with the weak directions in the plane of motion. The material of the beams and columns is steel with E = 210 GPa and ρ = 7850 kg/m3 . For this frame with 48 active degrees of freedom, the ten lowest natural frequencies were determined both through linearization and by improving methods. Calculations using the conventional FEM were also performed for models with two and four beam elements per structural member, giving 132 and 300 active degrees of freedom, respectively. The ten first natural frequencies obtained for the 48 degrees of freedom model, using various numbers of Mq matrices, are tabulated in Table 1. The exact results presented for comparison have been determined directly from the DSM formulation [7]. The results obtained by using the improving methods of Eqs. (14)–(16) are shown in Table 2, together with FEM results from the 132 and 300 degrees of freedom models. The relative CPU-times to obtain each result are also given in the tables. The CPU-time for the case with 48 degrees of freedom and q = 1 was 0.015 s using a PC with 2.66 GHz Intel Xeon processor and 1 Gb RAM. The assembly of the global stiffness matrix K0 and the mass matrices Mi from the polynomial element functions require approximately 0.005 s for each matrix. The example demonstrates that introducing more M-matrices gives results that converge towards the exact

1912

A. Ansell / Engineering Structures 27 (2005) 1906–1915

Table 2 Natural frequencies of the four-storey frame Freq. no.

48 d.o.f.a 132 d.o.f.a 300 d.o.f.a Impr.m. 1b Impr.m. 2c

1 2 3 4 5 6 7 8 9 10 Relative CPU-time

12.87 40.69 72.38 102.64 210.11 243.52 251.87 295.28 314.61 318.47

12.87 40.63 72.17 102.41 189.12 212.82 218.20 245.53 256.18 258.32

12.87 40.63 72.15 102.36 188.40 211.80 217.10 244.00 254.44 256.54

12.87 40.63 72.16 102.37 188.50 212.15 217.47 245.11 255.82 258.02

12.87 40.63 72.15 102.36 184.02 204.01 208.09 228.21 234.45 235.74

1.0

8.3

109.2

23.9

26.0

In rad/s, from conventional FEM analysis and computed with improving methods. a From FEM calculations. b Improving method no. 1: Eqs. (14) and (15). c Improving method no. 2: Eqs. (14) and (16).

results, as seen in Table 4. For higher frequencies, most vibration modes will converge towards the modes corresponding to the exact natural frequencies when using increasingly longer polynomials. This is due to the increasing flexibility of the elements gained when the number of correcting terms increases. A better description of the elements will thus give modes that better approximate the exact free vibration modes. It should be noted, however, that the small number of degrees of freedom gives limited information on the exact mode shapes compared to results from the conventional FEM with many degrees of freedom. The frequencies calculated with the improving method of Eq. (16) in Table 2 were all smaller than the exact value, and should therefore be regarded as less accurate than the results from Eq. (15). Adding further M-matrices to Eq. (16) did not improve the results. The number of additional Mmatrices should therefore be limited to only one, as in Eq. (15). The disadvantage of the improving technique is that an eigenvalue problem of similar size to the original must be solved for every frequency of interest. For this example, this means solving 11 times a linear eigenvalue problem of size 48 × 48 to obtain the improved values for the ten first frequencies. This should be compared to solving a linear eigenvalue problem of size (4 × 48) × (4 × 48) = 192 × 192 to produce results of similar accuracy as by Eq. (13). In Table 2, the results showed that a conventional finite element model with 132 degrees of freedom model give results less accurate than using q = 4 and that the 300 degrees of freedom model gave results more accurate than using q = 5. Comparison of the relative CPU-times in Tables 1 and 2 shows that solving the 300 degrees of freedom model requires almost six times as long calculation as for the DEM case with q = 5 while the 132 degrees of freedom model is close to that of the q = 4 formulation. It is also seen that the improving methods are less effective than using a DEM formulation. Note that the iterative method for DSM

Fig. 5. The bar-tower truss. The dimensions are given by L = 1.90 m.

analysis [7] requires more than 300 times as long CPUtime to determine the frequencies in Table 1 compared to the FEM with 48 identical degrees of freedom. Table 3 Cross-section of members in the bar-tower truss Member between node i– j

1–2, 3–4, 4–6, 5–6

1–4, 1–6, 2–3, 2–5, 3–7, 4–8, 5–9, 6–10

1–3, 1–5, 2–4, 2–6

3–8, 3–9, 4–7, 4–10, 5–7, 5–10, 6–8, 6–9

A (m2 )

3.2 × 10−4

3.2 × 10−3

1.3 × 10−3

1.0 × 10−2

I (m4 )

6.7 × 10−7

6.7 × 10−5

1.0 × 10−5

6.7 × 10−4

A and I are the area and moment of inertia, respectively, for a bar between nodes i and j.

5.2. A bar-tower truss This example is chosen to demonstrate the improvements possible for space trusses. The example was originally presented by Sadek [10] who compared a FEM model with a model of the (analytically) exact DSM form, also demonstrating the improving techniques for the DEM formulation. The example is chosen to demonstrate the convergence of the DEM, towards the DSM solution, as Mq matrices are added. The structure in Fig. 5 is a 3D truss tower composed of 25 bar elements. The properties of the material are E = 68.95 GPa and ρ = 2779 kg/m3 . The crosssections of the bars are given in Table 3. For each bar, the moments of inertia in the x y and the x z-planes were equal, I y = Iz = I . The five first natural frequencies, obtained with between one and five Mq matrices and solved with Eq. (13), are shown in Table 4 together with the improved and exact results. The polynomial form was much more efficient for this example than for structures composed of beam elements, as in the previous example. The large improvements for truss structures are due to the assumption that the bar elements have finite bending rigidities, i.e. that the elements do not have to remain straight during vibration. This is opposed

A. Ansell / Engineering Structures 27 (2005) 1906–1915

1913

Table 4 Natural frequencies of the bar-tower truss Freq. no. 1 2 3 4 5

No. of M-matricesa q=1 q=2

q=3

q =4

q=5

473.94 489.92 679.57 767.13 917.30

409.25 417.74 523.85 573.17 640.89

404.45 412.26 507.24 547.44 603.23

402.60 410.10 498.99 533.65 582.47

423.13 433.29 562.47 627.86 719.43

Impr. met.b

Exactc

412.54 435.24 549.04 625.33 706.94

401.28 408.48 487.70 509.95 574.85

In rad/s, computed with different M(ω). a Mass matrix defined as: M(ω) = q ω2(i−1) M . i i=1 b Improving method: Eqs. (14) and (15). c Results from DSM analysis [10].

for the lower cables H = 72 kN. The mast was modelled in 2D, using only three beam-elements and two cable-elements, giving 10 active degrees of freedom. Results provided by the polynomial form and solved with Eq. (13) are presented in Table 6. Conventional FEM models with 1, 2, 5, 10 and 30 beam elements per structural member resulted in the frequencies given for comparison in Table 7. Table 5 Measured natural frequencies of the guyed mast

Fig. 6. Overview and 2D view of the guyed mast.

to the conventional FEM, where infinite bending rigidity is assumed. This is clearly seen in Eq. (25) where the E I term for the M1 matrix is missing in the polynomial matrix functions. As in the previous example, the best results by the improving methods were obtained when using Eq. (14) and Eq. (15). It is also obvious that the results from using Eq. (13) are much better than the best improved results obtained by Eqs. (14)–(16). 5.3. A guyed mast In this example, all the natural frequencies in the interval 0 ≤ ω ≤ 2 Hz for the guyed mast shown in Fig. 6 were determined. Natural frequencies from measurements on the mast are put forth in [19], presenting the results from acceleration measurements on the tip of the mast and from a strain measurement at the root of the mast, here shown in Table 5. Note that there are some frequencies that were not measurable at these two points of the mast. The mast had a total height of 182 m with cables attached at 70 m and 143 m and consisted of a steel shell with a diameter of 1.6 m with the bending stiffness E I = 1800 MPa and the mass per unit length that was Aρ = 418 kg/m. The upper cables had areas and distributed masses that were A = 1000 mm2 and Aρ = 8.97 kg/m, respectively. For the lower cables, these data were A = 600 mm2 and Aρ = 5.35 kg/m. The cable tensions for the upper cables were H = 120 kN and

Freq. no.

Tip of masta

Root of mastb

1 2 3 4 5 6 7 8

0.28 0.35 0.44 – 0.72 0.94 1.16 –

0.28 0.35 0.44 0.55 0.73 0.94 – 1.55

In Hz, according to [19]. a From acceleration spectrum, measured at the tip of the mast. b From strain spectrum, measured at the root of the mast.

Table 6 Natural frequencies of the guyed mast Freq. no. 1 2 3 4 5 6 7 8

No. of M-matricesa q=1 q=2

q=3

q=4

q=5

0.28 0.43 0.56 – 1.03 2.06 – –

0.27 0.38 0.49 – 0.79 1.24 – –

0.27 0.38 0.48 – 0.77 1.18 – –

0.27 0.38 0.47 – 0.77 1.16 – –

0.27 0.39 0.51 – 0.84 1.42 – –

In Hz, computed with different M(ω). a Mass matrix defined as: M(ω) = q ω2(i−1) M . i i=1

The results in Table 6, from using the polynomial form, showed good convergence towards some of the measured frequencies in Table 5, for increasing number of included M-matrices. The results from calculations using 1, 2, 5, 10

1914

A. Ansell / Engineering Structures 27 (2005) 1906–1915

and 30 beam elements per structural member in Table 7 clearly demonstrate that using the conventional FEM gave convergence already for 5–10 elements per member. Table 7 Natural frequencies of the guyed mast Freq. no. 1 2 3 4 5 6 7 8

No. of elements per structural member 1 el/m. 2 el/m. 5 el/m. 10 el/m.

30 el/m.

0.28 0.43 0.56 – 1.03 2.06 – –

0.28 0.42 0.54 – 0.93 1.50 – –

0.28 0.42 0.54 – 0.94 1.51 – –

0.28 0.42 0.54 – 0.93 1.50 – –

0.28 0.42 0.54 – 0.93 1.50 – –

In Hz, from conventional FEM analysis.

However, the convergence was towards frequencies less accurate than the ones shown in Table 6, as the models, for increasing numbers of q or elements per structural members, converged towards the same vibration modes. The reason is that only including the M1 matrix, as in the FE formulation, gives a cable behaviour that is less flexible, resulting in higher vibration frequencies. 6. Conclusions The numerical examples demonstrate that the results from using the DEM were more accurate than with the conventional FEM. The method is more efficient for structures composed of bar elements than for structures of beam elements. This is due to the assumptions of infinite bending rigidities of the bars made in the conventional FEM. The results indicate that for these kinds of problems further discretisation is the best method if detailed descriptions of the vibration modes are important. If that is not the case and it is possible to use an efficient solution method for linear eigenvalue problems, the choice of longer DEM polynomials is the most efficient method with respect to CPU-time. If only a very small number of good approximative natural frequencies are sought, the improving techniques of Eqs. (14)–(16) can be competitive. The latter is also very easy to implement numerically. For the further research, it must be studied how the cable elements in DEM form behave in large cable supported structures, e.g. suspension and cable-stayed bridges. It is also of interest to investigate how the DEM can be efficiently implemented for solving dynamic response problems for frame and cable type structures subjected to e.g. wind and earthquake loads. As, in a mode superposition analysis, the frequency range of interest is given by the frequency content of the applied dynamic load, it should be investigated how various choices of frequency range affects the results. This should be done for each of the presented methods and DEM formulations.

Appendix. Exact element stiffness functions The exact stiffness functions for the beam element in Section 4.1 are: 1 α cos α 1 α f 2at (α) = 2 − 2 f1at (α) = 2 − 2 ω ω sin α ω ω sin α 1 β 3 (cosh β sin β + cos β sinh β) f1 (β) = 2 − ω 12(β) 1 β 3 (sinh β + sin β) f2 (β) = 2 − ω 12(β) 1 β 2 sin β sinh β f3 (β) = 2 − ω 6(β) 2 1 β (cosh − cos β) f4 (β) = 2 − 6(β) ω 1 β(cosh β sin β − cos β sinh β) f5 (β) = 2 − ω 4(β) 1 β(sinh β − sin β) f6 (β) = 2 − ω 2(β) with the parameter (β) = ω2 (1 − cos β cosh β). For the bar element in Section 4.2, the exact functions are: β 3 (sinh β cos β − cosh β sin β) 24ω2 sin β sinh β β 3 (sin β − sinh β) f8 (β) = − . 24ω2 sin β sinh β f7 (β) = −

For the cable element in Section 4.3 the exact functions are: m 1 (Ω ) =

12c22 (12 + c22 )c12 ω2

m 2 (Ω ) = − m 3 (Ω ) =

−

c2  2  c12 ω2 1 + c22 (κ−1) 2 Ω

c22 c1 (κ − 1)   2c12ω2 1 + c22 (κ−1) 2 Ω

κΩ 2 + c22 (κ − 1) 1 1  − + 2 (κ−1) 2 ω2 κω 2 4ω 1 + c2 Ω 2

m 4 (Ω ) = −

κΩ 2 + c22 (κ − 1) 1 1  + + 2 (κ−1) ω κω2 4ω2 1 + c22 Ω 2

1 κΩ 2 1 + − 2 2 ω 4ω κω2 2 1 κΩ 1 m 6 (Ω ) = − 2 + + 2 ω 4ω κω2 with the auxiliary function κ = tan(Ω /2)/(Ω /2). m 5 (Ω ) =

References [1] Fergusson NJ, Pilkey WD. Frequency-dependent element mass matrices. J Appl Mech 1993;59:136–9. [2] Fergusson NJ, Pilkey WD. Literature review of variants of the dynamic stiffness method, part 1: The dynamic element method. Shock Vib Dig 1993;25(2):3–12.

A. Ansell / Engineering Structures 27 (2005) 1906–1915 [3] Fergusson NJ, Pilkey WD. Literature review of variants of the dynamic stiffness method, part 2: Frequency-dependent matrix and other corrective methods. Shock Vib Dig 1993;25(4):3–10. [4] Banerjee JR. Dynamic stiffness formulation for structural elements: A general approach. Comp Struct 1997;63(1):101–3. [5] Banerjee JR, Sobey AJ. Dynamic stiffness formulation and free vibration analysis of a three-layered sandwich beam. Int J Solids Struct 2005;42(8):2181–97. [6] Mukherjee N, Chattopadhyay T. Improved free vibration analysis of stiffened plates by dynamic element method. Comp Struct 1994;52(2): 259–64. [7] Ansell A. Frequency dependent matrices for frame type structures. Lic thesis. Stockholm (Sweden): Dep Struct Eng, Royal Inst Techn (KTH); 1996. [8] Williams FW, Wittrick WH. Exact buckling and frequency calculations surveyed. J Struct Eng 1983;109:169–87. [9] Melosh RJ, Smith HA. New formulation for vibration analysis. J Eng Mech 1989;115:543–54. [10] Sadek EA. On the dynamics of framed structures. Comp Struct 1985; 20:1013–9.

1915

[11] Starossek U. Brückendynamik - Winderregte Schwingungen von Seilbrücken [Bridge Dynamics - Wind-Induced Vibration of CableSupported Bridges]. Ph.D. thesis. Germany: Univ of Stuttgart; 1991. [12] Starossek U. Dynamic stiffness matrix of sagging cable. J Eng Mech 1991;117(12):2815–29. [13] Kim J, Chang SP. Dynamic stiffness matrix of an inclined cable. Eng Struct 2001;23(12):1614–21. [14] Au FTK, Cheng YS, Cheung YK, Zheng DY. On the determination of natural frequencies and mode shapes of cable-stayed bridges. Appl Math Mod 2001;25(12):1099–115. [15] Peters G, Wilkinson JH. Ax = λB x and the generalized eigenproblem. SIAM J Numer Anal 1970;7:479–92. [16] Ruhe A. Algorithms for the nonlinear eigenvalue problem. SIAM J Numer Anal 1973;10:674–89. [17] Jennings A, McKeown JJ. Matrix computation. 2nd ed. Chichester (UK): John Wiley & Sons; 1993. [18] MATLAB: the language of technical computing. MATLAB notebook user’s guide. Natick (MA): MathWorks Inc; 1998. [19] Hirsch G, Nonhoff G, Volkmar H. Wind induced vibrations of guyed masts. Rotterdam (Netherlands): Balkema; 1990.