Vibration and elastic instability of thin-walled domes under uniform external pressure

Thin-Walled Structures Vol. 26, No. 3, pp. 159-177, 1996

Copyright © 1996 Elsevier Science Ltd Printed in Great Britain. All rights reserved

ELSEVIER

0263-8231/96 $15.00 PII:S0263-8231(96)00027-4

Vibration and Elastic Instability of Thin-Walled Domes under Uniform External Pressure

C. T. F. Ross Department of Mechanical and Manufacturing Engineering, University of Portsmouth, Portsmouth, U.K.

ABSTRACT Two thin-walled varying meridional curvature axisymmetric shell elements are presented for the vibration and elastic instability of thin-walled hemiellipsoidal domes under uniform external pressure. The theoretical analysis is an extension of previous work carried out by the author, where for the two elements presented in the present report, a cubic and a quadratic variation was assumed for the meridional and the circumferential displacements along the meridian of these elements. In the previous study, only linear variations were assumed for the meridional and circumferential displacements along the meridian of these elements. Comparisons were made between experiment and theory for both buckling and vibration of hemi-ellipsoidal shell domes, which varied from very flat oblate vessels to very long prolate vessels. In general, agreement between experiment and theory was good for the hemi-spherical dome and the prolate vessels, but not very good for the flat oblate vessels. Additionally, the two new elements gave poorer results than the original simpler element for the cantilever mode of vibration, but better results for the lobar modes of vibration. Copyright © 1996 Elsevier Science Ltd.

INTRODUCTION Thin-walled domes are usually used to block off the ends o f submarine pressure hulls; 1' 2 these d o m e ends can vary in shape from flat oblate forms to long prolate forms. 159

160

C. T. F. Ross

Fig. 1. Axisymmetric buckling of a flat oblate dome.

Under uniform external pressure, flat oblate domes can buckle axisymmetrically, with their noses denting inwards, as shown in Fig. 1. 3,4 In the cases of hemi-spherical domes and long prolate domes, the mode of buckling due to uniform external pressure, is of lobar form, as shown in Figs 2 and 3. Similarly, under periodic exciting forces, these vessels vibrate in numerous modes, but the modes of interest in the present study are those similar to the static buckling modes; that is, the vibration mode of interest for fiat oblate domes, is similar to the axisymmetric mode of buckling shown in Fig. 4, and for hemi-spherical domes and prolate hemi-ellipsoidal domes, the modes of interest, are similar to the lobar buckling modes of Fig. 5, where n = the total number of circumferential waves in the flank of the vessel; m = the number of half-waves in the meridional direction. In 1981, Ross & Mackney 3 presented a constant meridional curvature element (CMC), for the buckling of hemi-ellipsoidal domes under uniform external pressure, and in 1983, Ross & Johns 5 used this element to study the vibration of the same vessels under water. In 1990, Ross 4 presented a varying meridional curvature element (VMC), to extend this study; this element is shown in Fig. 6. Ross assumed the following displacement functions for the three local displacement, namely u, v and w.

Vibration and elastic instability of thin-walled domes

Fig. 2. Lobar buckling of a hemi-spherical dome.

Fig. 3. Lobar buckling of a tall prolate dome.

161

162

C. T. F. Ross

(o)

(b)

Fig. 4. Eigenmodes o f oblate hemi-ellipsoids. (a) First eigenmode, n = 0 a n d m = 2; (b) second eigenmode, n = 0 a n d m = 3. In the elevations, the d a s h e d curves show displacem e n t shapes a n d in the plans, n o d a l curves.

-,i

(o)

....

(b)

(c)

Fig. 5. L o b a r elgenmodes for hemi-ellipsoidal a n d prolate domes. (a) n = 1, m = 1; (b) n = 2, m = 2; (c) n = 4, m = 2. In the plans, the d a s h e d curves show displacement shapes a n d in the elevations n o d a l curves.

A ij-I -Bj Axisof;ymmetr~

v~

PartsectiononA-A

Fig. 6. V M C element.

Vibration and elastic instability of thin-walled domes (1 - 0 U

- -

-

u; cos n4~ + ~(1 + ~) u cos n~b

-

2 (1 - 0 V -- -

-

(~3

vi sin n4) +

2 _

163

3~ + 2)

W --

(_(3

-t-

+

(1 + 0 2 vj sin n4)

wi cos n~b +

(1 + 0 ( 1

(1) -

()2

Rl ~Oj cos n~)

4

3( + 2) (1 - ~) (1 + 0 2 4 wj cos nq~ -4 Rl~Ojcosnflp

where ( = s/(Ra) or, in matrix form,

(u} v

= [U]

W

{ Ui} = a matrix o f nodal displacements and [N] = a matrix o f shape functions = [N']* (cos n4~ or sin n4~). These displacement functions, which varied along the meridian o f the shell, assumed a linear variation for u and v, and a parabolic variation o f w. The assumptions for the displacement functions in the circumferential direction, 3~ was of a sinusoidal nature; this simplified computation in the circumferential direction, as integration in this direction was carried out explicitly. Integration in the meridional direction was carried out numerically, using four Gauss points per element. The stiffness matrix was given by: [k] ----

[B]T[D][B] dx dy =

= nRle

I'

[B]T[D][Blr ddp R10~ d~

-1

(2)

r[Bll T [Ol[B 1] d~

-1

where [B1]=[B]/(either constants; r = radius at

cosn~b

or

sinnq~);

[Dl=matrix

{e} = a vector of strains = [e,e4)e,4~KsK4~Ks4)]"r = [B]{Ui} where

of elastic (3)

C. T. F. Ross

164

1

Ou

gs----

w

q-

RI~ O~ Rl

e~ = -

+ u sin/3 + w cos 13

r

es~ = r

v sin/3 +

~ 0~

-1

02w

02/3

Ks - Rl2a 2 0 ~ 2

1 Ou

UOfis2+ R12a 0~

cos /3 0v

l[!omw

-

Ks~=-

2[-1

02w

sin/30w r 0~b

---,,)sin/31

1 (~ Ow

cos ,60v Rla 0(

sin/3 cos/3

l Ou) •

The relationship between local { Ui} and global { Ui ° } was given by: {Ui}

= [DC]

{u?} =

[UlU1WIOlU2U2W202]T.

(4)

The elemental stiffness matrix in global coordinates, namely [k°], was given by: [k °] = [DC]T[k] [DC]

(5)]

[coSll]o o4

0

[~1] =

[~2]

1

1

0

el

0

c2

0 0

0

1

--S 2

0

C2

0

0

0

(6)

=

and cl = c o s i l l ; c2 = c o s / 3 2 ; s1 ~-- sin/31; s2 = sin/32" The wall thickness was assumed to vary linearly along the meridian of the shell, as follows: (1 + ~ ) (1 ~) tl t2 t - ---2+-5--

where t~ = wall thickness at node 1; t2 = wall thickness at node 2.

(7)

Vibration and elastic instability of thin-walled domes

165

To obtain the geometrical stiffness for this element, it was necessary to consider the additional strains due to large deflections, 7 as follows:

6~,

~

(ow

b ~ - v cos/~ ~

~

)]

•

(8)

- o cos t~

N o w the additional strain energy, due to large deflections, is given by

½{u,}~[6]T[~{ u,} so that the geometrical stiffness matrix is given by [kl] =

I

t[a]T[a][O--]Rled~rd~ = R l c ~

t[G]T[a][G]rd~dq5 -1

0

(9)

= TrRlO~II_l t[G1]T[~l[Gllrd~ where [G 1] = [G]/(either cos n~b or sin nq~) [o] =

o~

•

F o r hemi-ellipsoidal thin-walled domes, the prebuckling m e m b r a n e stresses 8 were approximated by:

-pa 2 {7 s

t{2(a2 cos2/~ + b2 sin 2 fl)l/2} -pa2{b 2 - (a 2 - b 2) cos 2 fl}

(l 1)

6¢~ = t{2b2( a2 cos 2 fl + b 2 sin2 fl)l/2} where p = uniform pressure (external positive); a = radius of hemi-ellipsoidal d o m e at base; b = height of hemi-elliptical dome. The mass matrix was given by:

[m] = I I [N]T p[N] tr dqb R, ~cd~ (12)

= nRlep

tr[Nl] T [N l] d~ -1

where [N 1] = [N]/(either cos n~b or sin n4~).

C. T. F. Ross

166

In global co-ordinates,

[kl o] =

[DC]

(13)

[m °] = [DC]T[m] [DC].

(14)

[DC]T[kl]

and

For varying meridional curvature, it can be assumed that: (15)

fl =- ao + al s + a2 s 2.

Cook 9 shows that the three constants in eqn (15) namely, ao, al and a2 can be determined from the following three conditions: fl = fll at s = O

fl = flE at s = L

and sin (fl - flo)as ,~

(fl - tic)ds = 0

so that a0

=

al = (6tic -- 4fll -- 2fl2)

fll

a2 =

(3fll + 3fl2 - 6tic ) L2

where 1/R1 = -Ofl/Os; L = arc length of element; tic = chord angle. Cook shows that 2-~I~(fll--flc)~T[ 4 flc)j [ - 1

L"~l+vv[(fl2

--1]~(fll --flc)~ 4 t(fl2 flc)J

(16)

and r = R1 +

I s s i n fl d s o

and fl ~ tic + d y / d x

where

y

=

[x(1

-

x / l ) ( x 2 / l ) ( - 1 q- x/l)]

- flc)J

with l = chord length. For hemi-ellipsoids it is a simple matter to carry out the necessary integrations, and for axisymmetric elements of more complex shape, inte-

Vibration and elastic instability of thin-walled domes

167

gration can be carried out numerically with the aid of an additional midside node. It should be noted that for domes, Ow 0 -

(17)

u/R1. Os

One deficiency o f the V M C element was that it was assumed that u and v varied linearly in the meridional direction. Rajagopalan 1° pointed out this deficiency, and found its inclusion for circular cylinders was significant. However, in Rajagopalan's Chapter 6, not all the continuity requirements o f his polynomial were satisfied. Hence, the present study is reported.

THEORETICAL ANALYSIS Two elements are presented, one which assumes a cubic variation for u, v and w, along the meridian, (the CCC element), and one which assumes a quadratic variation for u and v, and a cubic variation for w, along the meridian, (the QQC element). The derivation of these elements is now described. The CCC element For this element, the displacement functions for u, v and w were all assumed to be o f cubic form. N o w the displacement function for w is of the Hermitian form shown by eqn (1), but the displacement functions for u and v will be obtained, as follows. For the node element o f Fig. 7, the assumed displacement functions for u and v are given by the following cubic forms: u = a + b ( + c~2 + d( 3 v = e +f(

+ g(2 + h(3.

To determine the shape functions for u and v, consider the b o u n d a r y conditions for u, as follows: at ~ = - 1 , u = ul .'. Ul = at~=-l,u=u3 at

1

~ = g, u = u4..

a-b+c-d;

" u3 = a - b / 3 + c / 9 - d / 2 7 ; u4 = a + b / 3 + c / 9 +

at(=l,u=u2.'.uz=a+b+c+d.

d/27;

168

C. T. F. Ross

¢=1/3 ~=-1/3 g=-I

/

J

i ¢=1

2 4

1 Fig. 7. N o d a l displacement positions for u a n d v of cubic form.

Solving the above four simultaneous equations, the following expressions are obtained for the constants: a=

- -]-~ (u2 -+- Ul) -]-

,,1

(u3 --]- u4)

9

b = - 1---2(u2 - Ul) +

(3u4 - 3u3)

(18)

t3 ¢ = ~ 6 ( u l -~ u2 -- u3 - - u4)

d

(/22 -- b/l -~- 3u3

= 9

- 3u4).

Similar expressions can be obtained for the constants from e to h, and substituting these constants into the displacement functions, the shape functions for u and v are given by: N1 = 1 (_ 1 + ( + 9( 2 - 9(3) N2 = 1 ( - 1 _ ( + 9( 2 + 9( 3) lO

(19)

fl

N3 = 1-~(1 - 3( - (2 +

3 ( 3)

N4 = 9 ( 1 + 3(

3(3),

- (2 _

i.e. u = (Nl.ul + N2.u2 + N3.u3 + N4.u4)cosn~p

(20)

and v = (Nl.ul + N2.u2 + N3.u3 + N4.u4) sinndp.

(21)

Vibration and elastic instability of thin-walled domes

169

These displacement functions can be seen to be very different to those used by Rajagopalan 1° in his Chapter 6. The [B] matrix was obtained from the strain-displacement relationship of eqn (3) where, DN1 -

DN2 --

ON1

o; ON2

a( (22)

DN3 D N 4 --

ON3 ON4

(23)

= L/(2R1)

L = arc length es = DN1/(RI~)Ul + (~3 _ 3( + 2)/(4R1)Wl

+ (1 + ~)(1 - ~)2~/4.01 + DN2/(RI~)u2 + (_(3 + 3( + 2)/(4R1)w2 - (1 - () (1 + ~2)~/4.02 + DN3/(Rl~)U3 + DN4/(RI~)U4 e¢ = N l . s i n f l / r . u i + n N 1 / ( r ) v l + (~3 _ 3~ + 2).cosfl/(4.r)wl

+ (1 + ~) (1

-

~2)2R 1 ~ COSfl/(4r).Ol

+ N2.sin fl/r.u2 + nN2/(r) U2 + (_(3 + 3~ + 2)cos fl/(4r).w2 - (1 - ~)(1 + ~)2RI c¢cos fl/(4r)02 + N3 sin fl/r.u3 + nN3/(r) v3 + N4 sin fl/r.u4 + nN4/(r) v4 es~ = - nN1/r.ul + [(-N1 sin fl + DNI.r/(R1 oO]/r.v~ -

nN2/r.u2 + [ ( - N 2 sin fl + DN2.r/(Rl~)]/r.vz

- nN3/r.u3 + [ ( - N 3 sin fl + DN3.r/(Rl~)]/r.v3 - nN4/r.u4 + [ ( - N 4 sin fl + DN4.r/(RI~)]/r.v4

170

C. T. F. Ross

Ks = [ D N 1 / ( R ~ ) - 2Nl.az] ul - 3~/(2R12o~2).wl - ( - 1 + 3~)/(2R1~)01 + [DN2/(R12~) - 2N2.a2] u2

+ 3~/(2R12~ 2) w2 - (1 + 3~)/(2R1~)02 + [DN3/(R12~) - 2N3.a2] u3 + [DN4/(RI2oO u4 - 2N4.a2]

KO = N1 .sin fl/(Rlr).ul + n.Nl.cos fl/(r 2) vl + [n2(~ 3 - 3~ + 2)/(4r z) - 3(~ 2 - 1)sin fl/(4Rl~r)] wl + [n2(1 + ~)(1 - ~)2Rl~/(4r2) + (1 + 2~ + 3ff2) sinfl/(4r)] 01 + N2.sin fl/(Rlr).u2 + n.N2.cos fl/(rZ).v2 + [n2(-~ 3 + 3~ + 2)/(4r 2) - 3 ( - ~ 2 ÷ 1)sin fl/(4Rl~r)] w2 + [-n2(1 - ~)(1 + ~)2.Rl~/(4r2) + (1 - 2ff - 3~ 2) sin fl/(4r)] 02 + N3.sin fl/(Rlr).u3 + n.N3.cos fl/(r 2) v3 + N4.sin fl/(Rlr).u4 + n.N4.cos fl/(r 2) 04 Ks4~ = - 2n.N1/(RI r).ul

+ [ - 2 sin fl cos fl.N1/r 2 + 2 . D N l . c o s fl/(R1~r)].vl + [ - n sin fl (~3 _ 3~ + 2)/(2r 2) + 3n.(~ 2 - 1)/(2Rlo~r)wl + [ - n sin fl (1 + ~) (1 - ~)2.Rl~/(Zr2) + n ( - 1 - 2~ + 3~2)/(2r)] 01 - 2n.N2/(Rlr) u2

+ [ - 2 sin fl cos fl N2/(r 2) + 2 DN2.cos fl/(Rl ~r)] v2 + [ - n sin fl (_~3 + 3~ + 2)/(2r 2) + 3 n ( - ~ 2 + 1)/(2Rl~r)] w2 + In sin fl (1 -- if) (1 + ~)2. R1 ~/(2r 2) - n(1 - 2~ - 3 ~2)/(2r)] 02 - 2n.N3/(Rlr) u3

+ [ - 2 sin fl.cos fl.N3/r 2 + 2 DN3.cos fl/(Rl~r)] v3 -

2n.N4/(Rlr) u4

+ [ - 2 sin ft. cos ~.N4/r 2 + 2 DN4.cos fl/(Rl~r)] v4. T h e g e o m e t r i c a l stiffness m a t r i x [kl] was o b t a i n e d f r o m eqn (9) a n d Table 1, a n d the m a s s m a t r i x [m] was o b t a i n e d f r o m eqn (12) a n d Table 2. It is evident t h a t all these matrices are of order 12 × 12, but to simplify c o m p u t a t i o n , G u y a n r e d u c t i o n 11 was used to eliminate the u3, v3, u4 a n d

Vibration and elastic instability o f thin-walled domes

171

o~

o

Y +

I

-q- ~s

~

~,:z

o

°

%

o~ ~o I

t

i

o I

172

C. T. F. Ross

U4 displacements in a m a n n e r similar to that o f Rajagopalan l°. Thus, these matrices will be o f order 8 × 8, so that the original programs for the V M C element could be used with minimal changes.

The QQC element F o r this element, the displacement functions for u and v were assumed to be in quadratic form along the meridian o f the element, as shown in Fig. 8. The positions of the nodal displacements for u and v are shown in Fig. 8. The assumed boundary conditions were at ( = - 1 , u

=

u 1

and v = vl;

at ( = 1, u = u2 and v = v2; at(=0,

(24)

u=u3 andv=v3.

This resulted in the following displacement functions for u and v u = N l . u l + N2.u2 + N3.u3

(25)

v = N l . v l + N2.v2 + N3.v3

(26)

where N1 = -½(1 - ( ) ( N2 = ½(1 + () (

(27)

N3 ----(1 - ~2) and their derivatives with respect to ~ are ¢=1

\

¢ = - 1 ~ /

3

1 Fig. 8. Three node varying meridional curvature element.

Vibration and elastic instability of thin-walled domes

173

ON1 DN1

-

DN2

-

DN3

--

ON2

(28)

ON3

Substitution of the above displacement functions, together with their derivatives into eqns (2), (9) and (12), will lead to the stiffness, geometrical and mass matrices for the QQC element. Here, again, Guyana reduction was used to eliminate the u3 and v3 displacements, so that the 10 x 10 matrices became of order 8 x 8. Hence it was a simple matter to modify the computer program for the CCC element to incorporate the QQC element.

C O M P U T A T I O N A L ANALYSIS In this section comparisons will be made for both buckling and vibration of the experimental results of 10 hemi-ellipsoidal domes (see Fig. 9), with the theoretical predictions of three varying meridional curvature elements. The models were made in solid urethane plastic (SUP), where the SUP liquid was poured between machined male and female aluminium alloy moulds and then thermoset. The models had an internal base diameter of 0.2m and a wall thickness of 0.002m. SUP was found to have the following properties: Young's modulus -- 2.89 x 109 N/m2; Poisson's ratio -- 0.3 assumed; density = 1230 kg/m 3.

Instability The three elements of Ross, namely the VMC element (linear-linear cubic), 4 the CCC element (all cubic) and the QQC element (quadraticquadratic-cubic), are compared with the experimentally obtained buckling pressures for the 10 hemi-ellipsoidal domes of Ross and Mackney 3 in Table 3, where AR

= aspect ratio -

dome height base radius "

F r o m Table 3 it can be seen that the CCC and the QQC elements give better results than the V M C element for all hemi-eUipsoidal domes, except for the oblate dome of aspect ratio 0.7. However, it must be emphasised that this dome is very prone to even the slightest variations in meridional

174

C. T. F. Ross

Fig. 9. Thin-walled hemi-ellipsoidal dome shells. TABLE 3 Buckling Pressures (MPa) for Hemi-ellipsoidal Domes AR

Experiment

4.0 3.5 3.0 2.5 2.0 1-5 1-0 0-7 0.44 0.25

0.097 0.100 0.150 0.201 0.364 0.490 1.280 1-070 0.276 0.060

VMC

0.111 0.135 0.168 0-224 0-323 0.549 1-426 0.944 0.380 0.132

(5) (5) (6) (7) (8) (9) (11) (0) (0) (0)

CCC

0-110 0.134 0.166 0-222 0.321 0.548 1.398 0.865 0.361 0.124

(5) (5) (6) (7) (7) (9) (11) (0) (0) (0)

QQC

0.110 0.134 0.166 0.222 0.321 0.548 1.398 0-866 0.376 0.124

(5) (5) (6) (7) (7) (9) (11) (0) (0) (0)

curvature. A slight increase in meridional curvature will result in a much higher buckling pressure. Most of the domes, especially the oblate domes of aspect ratio 0.44 and 0.25m seemed to have suffered from elastic knockdown, due to slight geometrical imperfections that occurred during manufacture.

Vibration and elastic instability of thin-walled domes

175

TABLE 4 Resonant Frequencies (Hz) in Air for Oblate Domes AR

0-7 0-44 0-25

n=0 Exp

VMC

CCC

QQC

1571 1104 836

1598 1115 711

1900 1271 783

1877 1267 782

Vibrations Comparisons are made in Tables 4 and 5 with the experimentally obtained results o f Ross and Johns, 5 and the three solutions o f Ross, for the 10 hemi-ellipsoidal domes described earlier. F r o m Table 4, it can be seen that the CCC and QQC elements give worse results than the simpler V M C element for the domes o f aspect ratio 0.7 and 0.44, but better results for the dome of aspect ratio 0.25. This may be due to the fact that all three elements adopted four degrees of freedom, namely, u, v, w and 0, but that for these oblate domes, v was zero. F r o m Table 5, it can be seen that the CCC and the QQC elements gave better results than the V M C elements for all modes, except for n = 1. It may be that, for the n -- 1 mode, v was of a much simpler form than assumed. The argument that can be held therefore, is that for n = 0 and n = 1, where v was either zero or of simple form, the assumption that v varied either quadratically or cubically was less justified than assuming that v varied linearly.

CONCLUSIONS The results have shown that by assuming quadratic and cubic variations with respect to the meridian o f hemi-ellipsoidal d o m e shells, for the meridional and circumferential displacements, good results are obtained for the majority of the vessels. The results have also shown that the m o r e sophisticated elements, the CCC element and the QQC element, give worse results than the simpler V M C element, when the circumferential displacement is zero, as in the cases when n = 0 and n = 1. The experimental results have shown that for the buckling o f thinwalled hemi-ellipsoidal domes under uniform external pressure, there was

176

C. T. F. Ross

H

O

e~

H

¢",1

H

Vibration and elastic instability of thin-walled domes

177

elastic k n o c k d o w n due to initial geometrical imperfections in the vessels, especially for the oblate d o m e o f A R = 0.25.

ACKNOWLEDGEMENTS The author would like to thank Prof. Jim Byrne and Dr M o n t y Bedford for encouraging this research. His thanks are extended to Miss Sharon Snook for the care and devotion she showed in typing this manuscript.

REFERENCES 1. Galletly, G. D. & Aylward, R. W., The influence of end closure shape on the buckling of cylinders under external pressure. Trans RINA 117 (1975) 255 268. 2. Blachut, J. & Galletly, G. D., Externally pressurised torispheres - - plastic buckling and collapse. In Buckling of Structures (Edited by I. Elishakoff et al.), pp. 29-45. Elsevier Science Publishers B.V, Amsterdam, 1988. 3. Ross, C. T. F. & Mackney, M. D. A., Deformation and stability studies of thin-walled domes under uniform pressure. Inst Phys Meeting on Stress Analysis of Fabricated Structures and Components, London 21 October, 1981; also J. Strain Analysis, 18 (1983) 167-172. 4. Ross, C. T. F., Pressure Vessels under External Pressure. Chapman & Hall, 1990. 5. Ross, C. T. F. & Johns, T., Vibration of submerged hemi-ellipsoidal domes. J. of Sound and Vib., 91 (1983) 363-373. 6. von Mises, R., Der Kritische aussendruck fur allseits belastate zylindrische rohre. Fest Zum 70, Geburstag von Prof. Dr A Stodola, Zurich, pp. 418-430, 1929. (Also USEMB Translation Report No 366, 1936.) 7. Stircklin, J. A., Haisler, W. E., Macdougall, H. R. & Stebbins, F. J., Nonlinear analysis for shells of revolution by the matrix displacement method. JAIAA, 6 (1968) 2306-2312. 8. Flugge, W., Stresses in Shells. Springer, Berlin, 1973. 9. Cook, R. D., Concepts and Applications of Finite Element Analysis, 2nd Ed. Wiley, New York 1981. 10. Rajagopalan, K., Finite Element Buckling Analysis of Stiffened Cylindrical Shells. Balkema, Rotterdam, 1993. 11. Guyan, R. J., Reduction of stiffness and mass matrices. JAIAA, 3 (1965) 380.