- Email: [email protected]
Diagonal mass and associated stiffness matrices for isoparametric 5-node quadrilateral and 9-node hexahedron elements
Finite Elements in Analysis and Design 12 (1993) 37-47
37
Elsevier FINEL 277
Diagonal mass and associated stiffness matrices for isoparametric 5-node quadrilater d and 9-node hexahedron elements Gerhard Sauer Technischer UberwachungsvereinBayern, Miinchen, Germany Received January 1992 Revised December 1992 Abstract. Expanded displacement interpolation functions that generate diagonal mass matrices are presented
for the isoparametric 5-node quadrilateral and 9-node hexabedron. By applying 2X2 or 2 x 2 x 2 Gauss integration, respectively, the modified interpolation functions contract the mass in the central node. The comer nodes are freed from masses. The massless d.o.f, of the comer nodes can be condensed out from the stiffness matrix. The size of the solution matrices is not increased by using elements with central nodes. In contrast, they become smaller. It is shown that the elements defined by the modified interpolation functions and the corresponding integration order are suited to satisfactorily reflect the dynamic structural properties.
Introduction
The competition between consistent and lumped mass matrices is almost as old as the finite element method. Only two years after Clougb had coined the term "finite element" in his famous paper [1] the basic concept of formulating consistent mass matrices was presented by Archer [2]. Although Archer did not explicitly mention the name "finite element", he already introduced the technique for calculating distributed mass influence matrices which is still in use today. In accordance with the Rayleigh-Ritz approach to modal analysis, Archer expressed the structural velocities by the product of displacement interpolation functions and generalized velocities. Putting this expression into the general formula for the kinetic energy, he obtained consistent mass matrices. The Rayleigh-Ritz approach used by Archer for beam elements was later identified to be a form of the displacement-based finite element method. Archer's idea opened the way to unrestrictedly tackle dynamic problems by the finite element method. The traditional scheme of lumping the structural mass in point masses could finally be replaced by a theoretically sounder approach; the calculation of mass and stiffness matrices was based on the same theoretical assumptions. With the formulation of consistent mass matrices it was possible to establish clear convergence criteria for eigenvalues yielded by finite element equations. They were found to be upper bounds to the exact solutions. Applying element displacement interpolation functions to the computation of the mass matrices constituted a significant progress. However, it was soon noted that even by employing the new consistent mass matrix the governing differential equations were approximated only roughly by the discrete equations. Przemieniecki [3] showed that the concept of Correspondence to: Gerhard Sauer, Technischer Uberwachungsverein Bayern, Westendstr. 199, D-8000 Miinchen 21, Germany. 0168-874X/93/$06.00 © 1993 - Elsevier Science Publishers B.V. All rights reserved
38
G. Sauer / Diagonal mass and associated stir]hess matrices
consistent mass matrices had to be extended further to fulfill the underlying differential equations with enhanced accuracy. He suggested to employ frequency-dependent interpolation functions to describe the structural motions. This idea was later advanced by Gupta [4,5] who generalised the idea of Przemieniecki from line elements to two- and three-dimensional elements. This theoretically very appealing idea suffers, however, from a serious disadvantagc: it leads to a quadratic eigenvalue problem that cannot be solved by standard numerical techniques. This disadvantage was probably one of the reasons the dynamic element method of Przemieniecki and Gupta did not find wide acceptance in the finite element community. The old ideas of lumped and frequency-independent consistent mass matrices still enjoy great attraction in finite element analysis. In many cases it is only a matter of taste which type of mass matrix is preferred. Using a finite element mesh with sufficient degrees of freedom, the lowest eigenvalues are found with satisfactory accuracy by both mass matrices. Because of the better theoretical foundations one may be induced to think that the consistent mass matrix is always superior to the lumped mass representation. But this is not the case. In some circumstances the lumped mass matrix delivers more exact eigensolutions than the consistent matrix. In Ref. [6] it is demonstrated, for example, that the lowest eigenfrequencies of an axially vibrating rod are better reproduced by a mass matrix in which the diagonal entries are increased at the expense of the off-diagonal terms than by the traditional consistent mass matrix. A more important field for the application of lumped mass matrices is discovered by examining the methods for the time integration of the equations of motion. One of these methods almost ultimately requires the use of a lumped mass matrix: it is the central difference method [7]. The great efficiency of the central difference method when used together with a diagonal mass matrix has motivated some attempts to find general procedures for direct computation of lumped mass matrices. The methods described in Refs. [8-10] are typical examples of these attempts. First of all they were developed to replace the intuitive procedure of converting a consistent mass matrix into a diagonal matrix. On this intuitive basis, the lumped mass matrix is formed by summing row-wise the mass terms of the consistent mass matrix and storing the sums in the corresponding diagonal matrix positions. The off-diagonal terms are set to zero. This procedure has, however, at least two drawbacks: firstly, it cannot be applied to rotatory mass terms and, secondly, it produces negative diagonal entries in the mass matrix if utilized with elements having mid-side nodes. A way to overcome the inherent limitations of the outlined procedure for diagonalizing a mass matrix is pointed out by Key and Beisinger [8]. They determine scaling factors by comparing exact eigenfrequencies with the frequencies produced by the spatial and time discretizations and calculate the diagonal rotatory mass terms via these scaling factors. An alternative method is proposed by Hinton, Rock and Zienkiewicz [9]. They propose to first compute the diagonal terms of the consistent mass matrix and scale them afterwards so that the total mass is preserved. This procedure always yields a positive definite mass matrix. Schreyer [10] suggests to generate a diagonal mass matrix and the corresponding stiffness matrix by employing orthogonal interpolation functions. This procedure requires the adoption of a mixed variational principle as the orthogonal interpolation functions overestimate the real structural stiffness by several magnitudes if inserted into the customary dynamic variational principle. A method to compute consistent diagonal mass and associated stiffness matrices for isoparametric 4-node quadrilateral and 8-node hexahedron elements is described in Ref. [11] ~. This method is based on an expansion of the usual interpolation functions and the selection of an appropriate integration order. It is demonstrated that the consistent diagonal 1 Equation (2) of [11] should read: a = 1.5(1
~/3).
39
G. Sauer / Diagonal mass and associated stiffness matrices
3
8
7
5
w
5
// 2
~
~
3
2 Fig.1.5-nodequadrilateraland9-nodehexahedron. mass and associated stiffness matrices more accurately approximate the exact eigensolutions than the traditional lumped and consistent matrices. A procedure to obtain diagonal mass and corresponding stiffness matrices for isoparametric 5-node quadrilateral and 9-node hexahedron elements is presented in the following. It is assumed that the last node is always located in the element's center. Again, this procedure exploits the enhanced flexibility given by expanded interpolation functions. This flexibility paves the way for diagonal mass matrices. The goal of this investigation is to find interpolation functions that sweep the mass to the central node. By contracting the total mass in the central node the other nodes become masseless. The d.o.f, of the corner nodes can then be eliminated during assembling of the global stiffness matrix, producing compact solutions matrices. Quadrilateral element with five nodes
If a central node is added to the four nodes of the simplest plane quadrilateral element, the displacement interpolation functions N are no longer bilinear functions of the natural coordinates r and s (Fig. 1). The fifth interpolation function is a quadratic function of r and s. The other four interpolation functions also become quadratically dependent on the natural co-ordinates [7]. To compute a consistent diagonal mass matrix, the functions N are further expanded by terms of a higher order. In this way the interpolation functions acquire a dependence on the third and fourth powers of r and s. The introduction of additional terms into the interpolation functions necessitates the adaptation of the summands which depend linearly on r or s. They are modified in such a way that the function N~ vanishes on all nodes except on node i. On node i it continues to have the value 1. This expansion provides for modified interpolation functions: N 1 =0.2511 - ( 1 - c t ) r - a r 3 ] [ 1 -
(1 - c t ) s
N2 = 0.25[1 + (1 - a ) r + ar3] [1 - (1 - ct)s
-
a s 3] -
0.25N5,
-
a s 3] -
0.25N5,
N3 = 0.25[1 + (1 - a ) r + a r 3 ] [1 + (1 - a ) s + ~ s 3] - 0 . 2 5 N ~ , N4 = 0.25[1 - (1 - a ) r - t ~ r 3 ] [1 + (1 - t ~ ) s
+ a s 3] -
N s = [ 1 - ( 1 - t ~ ) r 2 - ar4] [1 - (1 - a ) s 2 -
~S4].
0.25N5,
(1)
G. Sauer / Diagonal mass and associated stiffness matrices
40
Equations (1) contain the free parameter a whose value remains to be determined. This parameter serves as adjusting quantity. The aim is, of course, to adjust a so that the total mass is concentrated in the central node. The specific value for a depends on the integration order selected to compute the mass matrix. On examining eqns. (1), one observes that an exact integration of the mass matrix would require a 5 x 5 integration order in Gauss quadrature. Calculating the mass matrix in this way would involve a vast amount of computation work, which would largely be wasted as no value for a would be detected that diagonalizes the mass matrix in the desired way. Such a value can be arrived at if a 2 x 2 integration is applied. More specifically, if a is set to a = 1.5,
(2)
the element mass is represented by a diagonal 2 x 2 mass matrix. The total mass is accumulated in the two d.o.f, of the central node. The d.o.f, of the corner nodes are freed from any mass. But the interpolation functions (1) do not only render altered mass matrices; they also render new stiffness matrices. These stiffness matrices are associated with the 2 x 2 diagonal mass matrix; the mass matrix could be termed a consistent matrix. Unfortunately, the new stiffness matrix cannot be used without modification as it produces a singular equation system. It must therefore be modified to eliminate the spurious modes. This modification can easily be accomplished by superposing the new stiffness matrix K~ with the stiffness matrix K 2 resulting from (1) with a = 0. This second matrix K 2 is also computed by 2 x 2 Gauss integration. The calculation of the effective stiffness matrix can formally be expressed by K = f K , + (1 - f ) K 2
where B 1 signifies the strain-displacement matrix when a = 1.5, B z the strain-displacement matrix with a = 0, D the stress-strain matrix and V the volume of the element. The factor f must still be determined. This factor is selected so that the contribution of K 2 to the total stiffness matrix is closely related to a mass matrix that contracts the total mass in the central node. To find a factor that scales K 2 to this matrix, one will remember that a consistent mass matrix is calculated if the usual interpolation f u n c t i o n s - - a = 0 in ( 1 ) f a r e employed. Summing row-wise the elements of this mass matrix according to the above-mentioned lumping scheme, 4 / 9 of the total mass is concentrated in the central node. Now, if the dynamical properties of the element are only to be described by the central mass, omitting the masses of the corner nodes, the mass must either be increased by a factor of 9 / 4 or the stiffness matrix must be scaled by a factor of 4 / 9 . Bearing this in mind, the stiffness matrix scaled by the factor 4 / 9 can be thought of as being associated with a mass matrix defined only at the central node. The factor f in (3) can therefore be set to f = 5/9.
(4)
The quadrilateral element defined by the above interpolation functions and the 2 × 2 Gauss quadrature shall henceforth be identified as OUAOS.
Hexahedron element with nine nodes The development of a diagonal mass and a corresponding stiffness matrix for the isoparametric 9-node hexahedron element is analogous to the development of the mass and stiffness
G. Sauer / Diagonal mass and associated stiffness matrices
41
matrices for the quadrilateral element. It can therefore be described in less detail. The expanded interpolation functions are N l = 0.12511 - ( 1 - a ) r - a r 3 ] [ 1 X [1 - (1 -
3]
a ) t - cet3] -- 0.125N9,
5/2=0.12511 + (1
-a)r
+ ¢xr3] [1 - (1 - 0 0 s - o~s 3]
3] - 0.125N 9,
× [1- (1-a)t-at N3=0.125[1 + (1
- a ) r + a r 3 ] [ 1 + ( 1 - a ) s + a s 3]
×[1-(1-a)t-at
3] -0.125N9,
- a ) r - a r 3 ] [ 1 + (1 - a ) s +as 3]
N 4 =0.12511 - ( 1
3] - 0.125N 9,
× [1- (1-a)t-at N 5 = 0.12511- ( 1 - a ) r - a r × [1 + (1 -
(1 - a ) s - a s
3][1 - ( 1 - a ) s - a s
3]
a)t + at 3] - 0.125N9,
N6=0.12511+ ( 1 - a ) r × [1 + ( 1 - a ) t
+ar3][1-(1-0t)s-as
(5)
3]
+at 3] - 0.125N9,
N7=0.125[1 + ( 1 - a ) r
+ a r 3 ] [ 1 + ( X - a ) s + a s 3]
×[1 + ( 1 - o t ) t + a t 3] -0.125N9, Ns=0.125[1 - (1
- a ) r - a r 3 ] [1 + ( 1 - a ) s +as 3]
×[1 + ( 1 - a ) t
+ a t 3] -0.125N9,
Ng= [ 1 - ( 1 - a ) r 2 - a r 4 ] [ 1
-(1-
0~)$2 -- O~S4]
X[1--(1--Ot)t2--ott4],
where r, s and t are the natural co-ordinates of the hexahedron (Fig. 1). The mass matrix can be obtained from these interpolation functions by a 2 × 2 × 2 Gauss integration. The stiffness matrix is computed according to eqn. (3). The hexahedron element obtained by these assumptions shall be called HEXA9. A mass matrix that contracts the total mass in the central node of the 5-node quadrilateral or 9-node hexahedron can, however, also be generated if a 1-point integration with a = 0 in (1) is employed. The 1-point integration saves considerable computing time and outweighs the time required to determine the stiffness matrix. During the stiffness calculation the Gauss quadrature must be run through twice: In the first run the stiffness contributions due to B 1 are computed. The second run produces the stiffness matrix based on B E. In total, the time needed to generate the stiffness and mass matrices is not higher than the time required for computing the matrices of the traditional 5-node quadrilateral or 9-node hexahedron, respectively, which demand the adoption of 3 × 3 or 3 × 3 × 3 integration. It must now be tested whether the element formulation described above is capable of producing correct results when applied to constant-strain situations and to render improved solutions when used to study the dynamics of structures. A converging element must satisfy the first test in the limit case, i.e. when the element size tends towards zero [12,13]. To carry out this test a patch of elements of arbitrary shape is investigated. The element patch is subjected to boundary conditions that are equivalent to constant-stress states. If the elements
42
G. Sauer / Diagonal mass and associated stiffness matrices
\
\
•12
//, /f
\, !
8"
. . . . . . .
"13 • 9
~7
.1 11[
/
y
J/
.i; Fig. 2. Patch of elements for testing
x
2
QUADS.
in the patch reproduce the constant-stress state, then the element is said to be convergent. The patch test is only carried out for the 5-node quadrilateral.
Patch test for QUAD5
The QUAD5 element is tested on the patch depicted in Fig. 2. The material is assumed to be linear, isotropic elastic. The material's properties are: Young's modulus = 2 × 105 and Poisson's ratio = 0.3. Three different types of constant stress are considered. The boundary conditions are selceted in such a way that the constant stresses o-x, O-y and O~y are alternately produced. H e r e only the results for the stresses ~r and ~ry will be reported. The strains are calculated using either the strain-displacement matrix B 1 or the matrix B 2. In both cases identical results are obtained. In more general situations, however, it seems preferable to e m p l o y B 2 only. The patch test is performed in two steps. In the first step prescribed displacements corresponding to the constant stress states are specified for the patch. The element stresses and the reactions are computed in this step. In the second step only the minimal number of displacements required to suppress rigid-body motions is specified. The nodal forces obtained in the first step are applied as tractions. A constant stress ~x = 100 is induced when the element patch of Fig. 2 is subjected to the displacements u u x = 0.0005x,
Uy = - 0 . 0 0 0 1 5 y .
(6)
Constant shear stress Cry = 100 is caused by the displacements u s = 0.00065y,
Uy = 0.00065x.
(7)
The results for both computations are compiled in Tables 1 and 2. In both steps of the two computations the constant-stress states were reproduced. The nodal forces of both steps were also in agreement. Thus it may be concluded that the element OUAO5 is able to represent constant-stress conditions.
Convergence study
The convergence properties of the elements QUADS and HEXA9 can be evaluated by comparing the eigenfrequencies provided by these elements with exact solutions. To exclude
G. Sauer / Diagonal mass and associated stiffness matrices
43
Table 1 Computed nodal displacements and forces for o"x = 100 Node
1 2 3 4 5 6 7 8 9 10 11 12 13
Coordinates
Displacements
x
y
Ux
0.0000 1.0000 1.0500 0.0000 0.2000 0.9000 0.8000 0.3000 0.5500 0.5250 0.9375 0.5375 0.1250
0.0000 0.0000 0.8500 1.0000 0.2000 0.1000 0.7000 0.7000 0.4250 0.0750 0.4125 0.8125 0.4750
0.000000 0.500000 × 0.525000 X 0.000000 0.100000 x 0.450000 x 0.400000 x 0.150000 x 0.275000 x 0.262500 X 0.468750 x 0.268750 × 0.625000 x
Forces
uy 10- 3 10 T 3
--
10- 3 10- 3 10-3 10 - 3 10-3
-
-
--
--
-
10 -3
--
10- 3 10- 3 10-4
--
0.000000 0.000000 0.127500 x 0.150000 × 0.300000 x 0.150000 x 0.105000 x 0.105000 x 0.637500 x 0.112500 x 0.618750 × 0.121875 x 0.712500 ×
10 - 3 10-3 10-4 10 - 4 10-3 10 - 3 10-4 10 - 4 10-4 10- 3 10- 4
Fx
Fy
- 50.0 42.5 50.0 - 42.5 0,0 0.0 0.0 0.0 0.0 0.0 0,0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
the influence of discretization effects, it is preferable to investigate simple structures with regular meshes. For the elements dealt with here, it is easy to identify such structures: They are a square plate and a cube. Both structures are depicted in Fig. 3. Figure 3 also provides information on the dimensions and Poisson's ratio used in the calculations. It is assumed that both structures can vibrate freely. They are unrestrained . The plate can, of course, only vibrate in its plane. Its stresses obey plane stress conditions. To demonstrate the convergence properties of the new elements, the eigenfrequencies are computed with different meshes. The meshes are refined step-by-step. The solutions obtained with the elements OUAD5 and HEXA9 are compared with the eigenfrequencies rendered by the conventional consistent and diagonal mass matrices of the usual 4-node quadrilateral and 8-node hexahedron, respectively. These two matrices are based on the conventional bilinear interpolation functions. The diagonal mass matrix is constructed by summing up row-wise the terms of the consistent mass matrix as described above. All solutions are finally compared with the eigenfrequencies computed with fine finite element meshes. Exact analytical solu-
Table 2 Computed nodal displacements and forces for trxy = 100 Node
1 2 3 4 5 6 7 8 9 10 11 12 13
Coordinates
Displacements
Forces
X
y
Ux
Uy
F~
F~
0.0000 1.0000 1.0500 0.0000 0.2000 0.9000 0.8000 0.3000 0.5500 0.5250 0.9375 0.5375 0.1250
0.0000 0.0000 0.8500 1.0000 0.2000 0.1000 0.7000 0.7000 0.4250 0.0750 0.4125 0.8125 0.4750
0.000000 0.000000 × 0.552500 × 0.650000 × 0.130000× 0.650000 X 0.455000 × 0.455000 X 0.276250 × 0.487500 X 0.268125 X 0.528125 X 0.308750 X
0.000000 0.650000 × 0.682500 × 0.000000 0.130000× 0.585000 X 0.520000 X 0.195000 X 0.357500 x 0.341250 X 0.609375 X 0.349375 X 0.812500 X
- 50.0 - 52.5 50.0 52.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
- 50.0 42.5 50.0 - 42.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
10 - 3 10- 3 10- 3 10 - 3 10-4 10- 3 10- 3 10- 3 10 - - 4 10-3 10-3 10- 3
10 - 3 10 - 3 10 - 3 10- 3 10- 3 10- 3 10-3 10 - 3 10-3 10-3 10-4
G. Sauer / Diagonal mass and associated stiffness matrices
44
/\
s=de length o f plate and cube plate t h i c k n e s s = 1
l
~ 2
/ \ / \
/ \ / \ /
\
\
\
\
\
\
\
\
\ / \ ,/\ / \ / \ /\ / \/ \ / \ / \ / / \/
\ / / /
/
\ / /
\
/ /
\ / /
\
\/
Fig. 3. Square plate with 4 x 4 and cube with 4 x 4 x 4 mesh.
P o i s s o n ' s ratio = 0 3
tions for the free in-plane vibrations of a plate and a freely vibrating cube do not exist. It is therefore assumed that the lowest frequencies of these structures are reproduced accurately by the fine meshes.
Quadrilateral Solutions obtained for the square plate are shown in Figs. 4-6. Before discussing the test results, an important aspect must be pointed out. In Figs. 4-6, the solutions obtained with models of different d.o.f, are confronted. It must be remembered that, when applying OUADS, the d.o.f, pertaining to the corner nodes are condensed out while assembling the stiffness matrix. Only the d.o.f, of the central nodes are retained. The final solution matrix of a N x N mesh merely comprises 2 N 2 d.o.f. A model made up of 4-node quadrilaterals, however, has
Eigenfrequency
Eigenfrequency
[rad/s]
700
[rad/s]
700 [~
exact
QUAD5
650 ,
\
~
4node
usual consist.
-~
4node
usual diagonal
650 "\-,...
600
600
550
550
500
500
. . . . Lt~
/
450
400
/
__
i 2
/
450
/
/
-----
/
QUADS
/
I i
i
i
{
{
3
4
5
6
7
M e s h size
_i
400 8
m {
2
exact
4node
usual consist
4node
usual diagonal
i
{
t
3
4
5
N
[NxN]
Fig. 4. Convergence of the first non-rigid-body eigenfrequency of the square plate.
M e s h size
k____
6
i
7
N
[NxN]
Fig. 5. Convergence of the second
non-rigid-body eigenfrequency of the square plate.
45
G. Sauer / Diagonal mass and associated stiffness matrices Error in Frequency
[%!
12
8
Error in Frequency
[%1
14 12
m
10
- B - 8nodeusualconsist.
8
-~-
HEXA9
8nodeUsualdiagonal
[ ~
[]
4 4 2
0
0
-
-
--
-4 -4 -6
-8
-8
,\
-10
&
-12 -12 -14 -16 4
5
6
7
8
9
10 11 12 13 14 15 16 17 18
Frequency [Number] Fig. 6. Test of OUADS; errors in the lowest 15 nonrigid-body frequencies of the plate using a 7 x 7 mesh.
7
8
9
10 11 12 13 14 15 16 17 18 19 20 21
Frequency [Number] Fig. 7. Test of HEXA9, errors in the lowest 15 nonrigid-body frequencies of the cube using a 4 × 4 × 4 mesh.
2 ( N + 1) 2 d.o.f. Taking into account this difference, one can properly judge the solutions originating from the use of OUADS. Figures 4 and 5 show the convergence of the two lowest non-rigid-body eigenfrequencies when successively refining the finite element model. Both figures demonstrate that, at the level of coarse meshes, OUAD5 produces solutions closer to the excact values than the other two matrices. The first non-rigid-body eigenfrequency is already approximated with good accuracy at the stage of a 3 × 3 mesh. Advancing to still finer meshes, the first eigenfrequency determined with the usual 4-node diagonal matrix catches up with the solution yielded by OUADS. These two solutions contrast with the eigenfrequency originating from the usual 4-node consistent mass matrix. Even in a 7 × 7 mesh, this mass matrix cannot reproduce the exact value. The eigenfrequency still exhibits an error of 2%. The second non-rigid-body eigenfrequency is best approximated by OUAD5 with all meshes. With a 7 × 7 mesh OUAD5 misses the exact eigenfrequency by 1.4%. The errors in frequency for the fifteen lowest non-rigid-body eigenfrequencies computed with a 7 × 7 mesh are depicted in Fig. 6. Out of the fifteen frequencies a total of eleven frequencies is best approximated by OUADS. For the first and the ninth eigenfrequencies the results of the usual 4-node diagonal mass matrix a r e m o r e accurate. The conventional 4-node consistent mass matrix renders superior solutions for the eigenfrequencies seven and eighteen.
Hexahedron The observations regarding the convergence characteristics of the OUAD5 element can, in principle, be transferred to the HEXA9 element. Typical results for this element are summarized in Table 3 and plotted in Fig. 7. Observing these results, it must again be r e m e m b e r e d that solutions obtained by finite element meshes with different d.o.f, are compared. Employ-
46
G. Sauer / Diagonal mass and associated stiffness matrices
Table 3 Lowest four non-rigid body frequencies of the cube computed with different meshes; multiple frequencies are skipped a Mesh
Mass matrix
oJ 7
co 9
t012
2x 2x 2
HEXA9 8-node usual consistent 8-node usual diagonal
1.0181
1.1363
1.1689
1.1839
1.0742 0.7161
1.3994 0.8523
1.5191 0.8771
1.5191
HEXA9 8-node usual consistent 8-node usual diagonal
O.9606
1. 2225
1. 2405
1.2580
0.9743 0.7978
1.3197 1.0133
1.3517 1.0527
1.4656
HEXA9
O.9308
1.2264
1.2336
1.3075
8-node usual consistent 8-node usual diagonal
0.9375 0.8338
1.2722 1.0816
1.2937 1.1573
1.4353
(/.8945
1.2071
1.2266
1.3890
3x 3x 3
4X4X4
10 x 10 x 10
o315
(I.9162
1.1360
1.2046
a coi = ~oi / EV~-~ ; the closest approximations to the exact values in italics. ing HEXA9, a finite e l e m e n t m e s h o f a N x N x N subdivision c o n d e n s e s to a system o f 3 N 3 e q u a t i o n s which, for l a r g e r N, signifies a c o n s i d e r a b l e r e d u c t i o n in the n u m b e r o f equations. T h e usual 8 - n 0 d e e l e m e n t s b u i l d u p a system o f 3 ( N + 1) 3 equations. HEXA9 c o u l d also be very a d v a n t a g e o u s f r o m t h e p o i n t o f view o f s m a l l e r s o l u t i o n m a t r i c e s , p r o v i d e d it yields accurate eigensolutions. O n i n s p e c t i n g T a b l e 3, o n e n o t e s t h a t this p r e r e q u i s t e is fulfilled. Surprisingly g o o d results a r e o b t a i n e d for t h e lowest frequencies. T h e lowest f r e q u e n c i e s u p to the 14th f r e q u e n c y a r e m o s t closely a p p r o x i m a t e d by HEXA9 with all meshes. F i g u r e 7 d e m o n s t r a t e s that only two e i g e n f r e q u e n c i e s a r e b e t t e r r e p r o d u c e d by o n e of t h e o t h e r two matrices. T h e s e two i d e n t i c a l f r e q u e n c i e s a r e c o m p u t e d m o r e a c c u r a t e l y by t h e u s u a l 8 - n o d e c o n s i s t e n t m a s s matrix. T h e solutions f o u n d by using t h e usual 8 - n o d e d i a g o n a l mass m a t r i x a r e less exact t h a n t h o s e p r o d u c e d by HEXA9 for t h e w h o l e f r e q u e n c y r a n g e c o n s i d e r e d .
Conclusion Q u a d r i l a t e r a l a n d h e x a h e d r o n e l e m e n t s with c e n t r a l n o d e s can be u s e d effectively to c o m p u t e e i g e n f r e q u e n c i e s a n d m o d e s h a p e s o f l i n e a r elastic structures. By c h o o s i n g a p p r o p r i ate i n t e r p o l a t i o n f u n c t i o n s a n d an a s s o c i a t e d i n t e g r a t i o n scheme, d i a g o n a l mass m a t r i c e s a r e g e n e r a t e d . Since t h e y g a t h e r t h e m a s s on t h e c e n t r a l n o d e , t h e d.o.f, o f t h e massless c o r n e r n o d e s can b e c o n d e n s e d o u t f r o m t h e stiffness matrix. T h e size of t h e solutions m a t r i c e s is r e d u c e d a n d now d e p e n d s on t h e n u m b e r o f e l e m e n t s r a t h e r t h a n on t h e n u m b e r o f nodes. T h e c o n v e r g e n c e c h a r a c t e r i s t i c s o f t h e s u g g e s t e d e l e m e n t s a r e satisfactory. T h e s e e l e m e n t s a p p r o x i m a t e t h e exact solutions m o r e closely t h a n e l e m e n t s t h a t e m p l o y t h e usual c o n s i s t e n t o r t h e u s u a l d i a g o n a l m a s s matrices. T h e m a t r i c e s b a s e d on t h e e x p a n d e d i n t e r p o l a t i o n f u n c t i o n s a r e s u i t e d to m o d e l l i n g t h e d y n a m i c s t r u c t u r a l p r o p e r t i e s .
References [1] R.W. C L O U G H , "The finite element in plane stress analysis", Proc. 2nd A S C E Conf. on Electronic Computation, Pittsburgh, PA, September 1960. [2] J.S. ARCHER, "Consistent mass matrix for distributed mass systems", J. Struct. Div. A S C E 89 (ST 4), pp. 161-178, 1963.
G. Sauer / Diagonal mass and associated stiffness matrices
47
[3] J.S. PRZEMIENIECKI, Theory of Matrix Structural Analysis, McGraw-Hill, New York, 1968. [4] K.K. GUPTA, "Numerical formulation for a higher order plane finite dynamic element", Int. J. Numer. Methods Eng. 20, pp. 1407-1414, 1984. [5] K.K. GUPTA, "Development of a solid hexahedron finite dynamic element", Int. J. Numer. Methods Eng. 20, pp. 2143-2150, 1984. [6] G. SAUER and M. WOLF, " A modified beam element mass matrix", Z. Angew. Math. Mech. 68, pp. 483-490, 1988. [7] K.J. BATHE and E.L. WILSON, Numerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1976. [8] S.W. KEY and Z.E. BEISINGER,"Slade D: A computer program for the dynamic analysis of thin shells", Report SLA-73-0079, Sandia Laboratories, Albuqerque, NM, 1973. [9] E. HINTON, T. ROCK and O.C. ZIENKmWlCZ, " A note on mass lumping and related processes in the finite element method", Int. J. Earthquake Eng. Struct. Dyn. 4, pp. 245-249, 1976. [10] H.L. SCHREYER,"Consistent diagonal mass matrices and finite element equations for one-dimensional problems", Int. J. Numer. Methods Eng. 12, pp. 1171-1184, 1978. [11] G. SAUER,"Consistent diagonal mass matrices for the isoparametric 4-node quadrilateral and 8-node hexahedron elements", Commun. Numer. Methods Eng. 9, pp. 35-43, 1993. [12] R.L. TAYLOR,J.C. SIMO,O.C. ZIENKIEWICZ and A.C.H. CrtAN, "The patch t e s t - - A condition for assessing fem convergence", Int. J. Numer. Methods Eng. 22, pp. 39-62, 1986. [13] A. RAZZAOtJE, "The patch test for elements", Int. J. Numer. Methods Eng. 22, pp. 63-71, 1986.
37
Elsevier FINEL 277
Diagonal mass and associated stiffness matrices for isoparametric 5-node quadrilater d and 9-node hexahedron elements Gerhard Sauer Technischer UberwachungsvereinBayern, Miinchen, Germany Received January 1992 Revised December 1992 Abstract. Expanded displacement interpolation functions that generate diagonal mass matrices are presented
for the isoparametric 5-node quadrilateral and 9-node hexabedron. By applying 2X2 or 2 x 2 x 2 Gauss integration, respectively, the modified interpolation functions contract the mass in the central node. The comer nodes are freed from masses. The massless d.o.f, of the comer nodes can be condensed out from the stiffness matrix. The size of the solution matrices is not increased by using elements with central nodes. In contrast, they become smaller. It is shown that the elements defined by the modified interpolation functions and the corresponding integration order are suited to satisfactorily reflect the dynamic structural properties.
Introduction
The competition between consistent and lumped mass matrices is almost as old as the finite element method. Only two years after Clougb had coined the term "finite element" in his famous paper [1] the basic concept of formulating consistent mass matrices was presented by Archer [2]. Although Archer did not explicitly mention the name "finite element", he already introduced the technique for calculating distributed mass influence matrices which is still in use today. In accordance with the Rayleigh-Ritz approach to modal analysis, Archer expressed the structural velocities by the product of displacement interpolation functions and generalized velocities. Putting this expression into the general formula for the kinetic energy, he obtained consistent mass matrices. The Rayleigh-Ritz approach used by Archer for beam elements was later identified to be a form of the displacement-based finite element method. Archer's idea opened the way to unrestrictedly tackle dynamic problems by the finite element method. The traditional scheme of lumping the structural mass in point masses could finally be replaced by a theoretically sounder approach; the calculation of mass and stiffness matrices was based on the same theoretical assumptions. With the formulation of consistent mass matrices it was possible to establish clear convergence criteria for eigenvalues yielded by finite element equations. They were found to be upper bounds to the exact solutions. Applying element displacement interpolation functions to the computation of the mass matrices constituted a significant progress. However, it was soon noted that even by employing the new consistent mass matrix the governing differential equations were approximated only roughly by the discrete equations. Przemieniecki [3] showed that the concept of Correspondence to: Gerhard Sauer, Technischer Uberwachungsverein Bayern, Westendstr. 199, D-8000 Miinchen 21, Germany. 0168-874X/93/$06.00 © 1993 - Elsevier Science Publishers B.V. All rights reserved
38
G. Sauer / Diagonal mass and associated stir]hess matrices
consistent mass matrices had to be extended further to fulfill the underlying differential equations with enhanced accuracy. He suggested to employ frequency-dependent interpolation functions to describe the structural motions. This idea was later advanced by Gupta [4,5] who generalised the idea of Przemieniecki from line elements to two- and three-dimensional elements. This theoretically very appealing idea suffers, however, from a serious disadvantagc: it leads to a quadratic eigenvalue problem that cannot be solved by standard numerical techniques. This disadvantage was probably one of the reasons the dynamic element method of Przemieniecki and Gupta did not find wide acceptance in the finite element community. The old ideas of lumped and frequency-independent consistent mass matrices still enjoy great attraction in finite element analysis. In many cases it is only a matter of taste which type of mass matrix is preferred. Using a finite element mesh with sufficient degrees of freedom, the lowest eigenvalues are found with satisfactory accuracy by both mass matrices. Because of the better theoretical foundations one may be induced to think that the consistent mass matrix is always superior to the lumped mass representation. But this is not the case. In some circumstances the lumped mass matrix delivers more exact eigensolutions than the consistent matrix. In Ref. [6] it is demonstrated, for example, that the lowest eigenfrequencies of an axially vibrating rod are better reproduced by a mass matrix in which the diagonal entries are increased at the expense of the off-diagonal terms than by the traditional consistent mass matrix. A more important field for the application of lumped mass matrices is discovered by examining the methods for the time integration of the equations of motion. One of these methods almost ultimately requires the use of a lumped mass matrix: it is the central difference method [7]. The great efficiency of the central difference method when used together with a diagonal mass matrix has motivated some attempts to find general procedures for direct computation of lumped mass matrices. The methods described in Refs. [8-10] are typical examples of these attempts. First of all they were developed to replace the intuitive procedure of converting a consistent mass matrix into a diagonal matrix. On this intuitive basis, the lumped mass matrix is formed by summing row-wise the mass terms of the consistent mass matrix and storing the sums in the corresponding diagonal matrix positions. The off-diagonal terms are set to zero. This procedure has, however, at least two drawbacks: firstly, it cannot be applied to rotatory mass terms and, secondly, it produces negative diagonal entries in the mass matrix if utilized with elements having mid-side nodes. A way to overcome the inherent limitations of the outlined procedure for diagonalizing a mass matrix is pointed out by Key and Beisinger [8]. They determine scaling factors by comparing exact eigenfrequencies with the frequencies produced by the spatial and time discretizations and calculate the diagonal rotatory mass terms via these scaling factors. An alternative method is proposed by Hinton, Rock and Zienkiewicz [9]. They propose to first compute the diagonal terms of the consistent mass matrix and scale them afterwards so that the total mass is preserved. This procedure always yields a positive definite mass matrix. Schreyer [10] suggests to generate a diagonal mass matrix and the corresponding stiffness matrix by employing orthogonal interpolation functions. This procedure requires the adoption of a mixed variational principle as the orthogonal interpolation functions overestimate the real structural stiffness by several magnitudes if inserted into the customary dynamic variational principle. A method to compute consistent diagonal mass and associated stiffness matrices for isoparametric 4-node quadrilateral and 8-node hexahedron elements is described in Ref. [11] ~. This method is based on an expansion of the usual interpolation functions and the selection of an appropriate integration order. It is demonstrated that the consistent diagonal 1 Equation (2) of [11] should read: a = 1.5(1
~/3).
39
G. Sauer / Diagonal mass and associated stiffness matrices
3
8
7
5
w
5
// 2
~
~
3
2 Fig.1.5-nodequadrilateraland9-nodehexahedron. mass and associated stiffness matrices more accurately approximate the exact eigensolutions than the traditional lumped and consistent matrices. A procedure to obtain diagonal mass and corresponding stiffness matrices for isoparametric 5-node quadrilateral and 9-node hexahedron elements is presented in the following. It is assumed that the last node is always located in the element's center. Again, this procedure exploits the enhanced flexibility given by expanded interpolation functions. This flexibility paves the way for diagonal mass matrices. The goal of this investigation is to find interpolation functions that sweep the mass to the central node. By contracting the total mass in the central node the other nodes become masseless. The d.o.f, of the corner nodes can then be eliminated during assembling of the global stiffness matrix, producing compact solutions matrices. Quadrilateral element with five nodes
If a central node is added to the four nodes of the simplest plane quadrilateral element, the displacement interpolation functions N are no longer bilinear functions of the natural coordinates r and s (Fig. 1). The fifth interpolation function is a quadratic function of r and s. The other four interpolation functions also become quadratically dependent on the natural co-ordinates [7]. To compute a consistent diagonal mass matrix, the functions N are further expanded by terms of a higher order. In this way the interpolation functions acquire a dependence on the third and fourth powers of r and s. The introduction of additional terms into the interpolation functions necessitates the adaptation of the summands which depend linearly on r or s. They are modified in such a way that the function N~ vanishes on all nodes except on node i. On node i it continues to have the value 1. This expansion provides for modified interpolation functions: N 1 =0.2511 - ( 1 - c t ) r - a r 3 ] [ 1 -
(1 - c t ) s
N2 = 0.25[1 + (1 - a ) r + ar3] [1 - (1 - ct)s
-
a s 3] -
0.25N5,
-
a s 3] -
0.25N5,
N3 = 0.25[1 + (1 - a ) r + a r 3 ] [1 + (1 - a ) s + ~ s 3] - 0 . 2 5 N ~ , N4 = 0.25[1 - (1 - a ) r - t ~ r 3 ] [1 + (1 - t ~ ) s
+ a s 3] -
N s = [ 1 - ( 1 - t ~ ) r 2 - ar4] [1 - (1 - a ) s 2 -
~S4].
0.25N5,
(1)
G. Sauer / Diagonal mass and associated stiffness matrices
40
Equations (1) contain the free parameter a whose value remains to be determined. This parameter serves as adjusting quantity. The aim is, of course, to adjust a so that the total mass is concentrated in the central node. The specific value for a depends on the integration order selected to compute the mass matrix. On examining eqns. (1), one observes that an exact integration of the mass matrix would require a 5 x 5 integration order in Gauss quadrature. Calculating the mass matrix in this way would involve a vast amount of computation work, which would largely be wasted as no value for a would be detected that diagonalizes the mass matrix in the desired way. Such a value can be arrived at if a 2 x 2 integration is applied. More specifically, if a is set to a = 1.5,
(2)
the element mass is represented by a diagonal 2 x 2 mass matrix. The total mass is accumulated in the two d.o.f, of the central node. The d.o.f, of the corner nodes are freed from any mass. But the interpolation functions (1) do not only render altered mass matrices; they also render new stiffness matrices. These stiffness matrices are associated with the 2 x 2 diagonal mass matrix; the mass matrix could be termed a consistent matrix. Unfortunately, the new stiffness matrix cannot be used without modification as it produces a singular equation system. It must therefore be modified to eliminate the spurious modes. This modification can easily be accomplished by superposing the new stiffness matrix K~ with the stiffness matrix K 2 resulting from (1) with a = 0. This second matrix K 2 is also computed by 2 x 2 Gauss integration. The calculation of the effective stiffness matrix can formally be expressed by K = f K , + (1 - f ) K 2
where B 1 signifies the strain-displacement matrix when a = 1.5, B z the strain-displacement matrix with a = 0, D the stress-strain matrix and V the volume of the element. The factor f must still be determined. This factor is selected so that the contribution of K 2 to the total stiffness matrix is closely related to a mass matrix that contracts the total mass in the central node. To find a factor that scales K 2 to this matrix, one will remember that a consistent mass matrix is calculated if the usual interpolation f u n c t i o n s - - a = 0 in ( 1 ) f a r e employed. Summing row-wise the elements of this mass matrix according to the above-mentioned lumping scheme, 4 / 9 of the total mass is concentrated in the central node. Now, if the dynamical properties of the element are only to be described by the central mass, omitting the masses of the corner nodes, the mass must either be increased by a factor of 9 / 4 or the stiffness matrix must be scaled by a factor of 4 / 9 . Bearing this in mind, the stiffness matrix scaled by the factor 4 / 9 can be thought of as being associated with a mass matrix defined only at the central node. The factor f in (3) can therefore be set to f = 5/9.
(4)
The quadrilateral element defined by the above interpolation functions and the 2 × 2 Gauss quadrature shall henceforth be identified as OUAOS.
Hexahedron element with nine nodes The development of a diagonal mass and a corresponding stiffness matrix for the isoparametric 9-node hexahedron element is analogous to the development of the mass and stiffness
G. Sauer / Diagonal mass and associated stiffness matrices
41
matrices for the quadrilateral element. It can therefore be described in less detail. The expanded interpolation functions are N l = 0.12511 - ( 1 - a ) r - a r 3 ] [ 1 X [1 - (1 -
3]
a ) t - cet3] -- 0.125N9,
5/2=0.12511 + (1
-a)r
+ ¢xr3] [1 - (1 - 0 0 s - o~s 3]
3] - 0.125N 9,
× [1- (1-a)t-at N3=0.125[1 + (1
- a ) r + a r 3 ] [ 1 + ( 1 - a ) s + a s 3]
×[1-(1-a)t-at
3] -0.125N9,
- a ) r - a r 3 ] [ 1 + (1 - a ) s +as 3]
N 4 =0.12511 - ( 1
3] - 0.125N 9,
× [1- (1-a)t-at N 5 = 0.12511- ( 1 - a ) r - a r × [1 + (1 -
(1 - a ) s - a s
3][1 - ( 1 - a ) s - a s
3]
a)t + at 3] - 0.125N9,
N6=0.12511+ ( 1 - a ) r × [1 + ( 1 - a ) t
+ar3][1-(1-0t)s-as
(5)
3]
+at 3] - 0.125N9,
N7=0.125[1 + ( 1 - a ) r
+ a r 3 ] [ 1 + ( X - a ) s + a s 3]
×[1 + ( 1 - o t ) t + a t 3] -0.125N9, Ns=0.125[1 - (1
- a ) r - a r 3 ] [1 + ( 1 - a ) s +as 3]
×[1 + ( 1 - a ) t
+ a t 3] -0.125N9,
Ng= [ 1 - ( 1 - a ) r 2 - a r 4 ] [ 1
-(1-
0~)$2 -- O~S4]
X[1--(1--Ot)t2--ott4],
where r, s and t are the natural co-ordinates of the hexahedron (Fig. 1). The mass matrix can be obtained from these interpolation functions by a 2 × 2 × 2 Gauss integration. The stiffness matrix is computed according to eqn. (3). The hexahedron element obtained by these assumptions shall be called HEXA9. A mass matrix that contracts the total mass in the central node of the 5-node quadrilateral or 9-node hexahedron can, however, also be generated if a 1-point integration with a = 0 in (1) is employed. The 1-point integration saves considerable computing time and outweighs the time required to determine the stiffness matrix. During the stiffness calculation the Gauss quadrature must be run through twice: In the first run the stiffness contributions due to B 1 are computed. The second run produces the stiffness matrix based on B E. In total, the time needed to generate the stiffness and mass matrices is not higher than the time required for computing the matrices of the traditional 5-node quadrilateral or 9-node hexahedron, respectively, which demand the adoption of 3 × 3 or 3 × 3 × 3 integration. It must now be tested whether the element formulation described above is capable of producing correct results when applied to constant-strain situations and to render improved solutions when used to study the dynamics of structures. A converging element must satisfy the first test in the limit case, i.e. when the element size tends towards zero [12,13]. To carry out this test a patch of elements of arbitrary shape is investigated. The element patch is subjected to boundary conditions that are equivalent to constant-stress states. If the elements
42
G. Sauer / Diagonal mass and associated stiffness matrices
\
\
•12
//, /f
\, !
8"
. . . . . . .
"13 • 9
~7
.1 11[
/
y
J/
.i; Fig. 2. Patch of elements for testing
x
2
QUADS.
in the patch reproduce the constant-stress state, then the element is said to be convergent. The patch test is only carried out for the 5-node quadrilateral.
Patch test for QUAD5
The QUAD5 element is tested on the patch depicted in Fig. 2. The material is assumed to be linear, isotropic elastic. The material's properties are: Young's modulus = 2 × 105 and Poisson's ratio = 0.3. Three different types of constant stress are considered. The boundary conditions are selceted in such a way that the constant stresses o-x, O-y and O~y are alternately produced. H e r e only the results for the stresses ~r and ~ry will be reported. The strains are calculated using either the strain-displacement matrix B 1 or the matrix B 2. In both cases identical results are obtained. In more general situations, however, it seems preferable to e m p l o y B 2 only. The patch test is performed in two steps. In the first step prescribed displacements corresponding to the constant stress states are specified for the patch. The element stresses and the reactions are computed in this step. In the second step only the minimal number of displacements required to suppress rigid-body motions is specified. The nodal forces obtained in the first step are applied as tractions. A constant stress ~x = 100 is induced when the element patch of Fig. 2 is subjected to the displacements u u x = 0.0005x,
Uy = - 0 . 0 0 0 1 5 y .
(6)
Constant shear stress Cry = 100 is caused by the displacements u s = 0.00065y,
Uy = 0.00065x.
(7)
The results for both computations are compiled in Tables 1 and 2. In both steps of the two computations the constant-stress states were reproduced. The nodal forces of both steps were also in agreement. Thus it may be concluded that the element OUAO5 is able to represent constant-stress conditions.
Convergence study
The convergence properties of the elements QUADS and HEXA9 can be evaluated by comparing the eigenfrequencies provided by these elements with exact solutions. To exclude
G. Sauer / Diagonal mass and associated stiffness matrices
43
Table 1 Computed nodal displacements and forces for o"x = 100 Node
1 2 3 4 5 6 7 8 9 10 11 12 13
Coordinates
Displacements
x
y
Ux
0.0000 1.0000 1.0500 0.0000 0.2000 0.9000 0.8000 0.3000 0.5500 0.5250 0.9375 0.5375 0.1250
0.0000 0.0000 0.8500 1.0000 0.2000 0.1000 0.7000 0.7000 0.4250 0.0750 0.4125 0.8125 0.4750
0.000000 0.500000 × 0.525000 X 0.000000 0.100000 x 0.450000 x 0.400000 x 0.150000 x 0.275000 x 0.262500 X 0.468750 x 0.268750 × 0.625000 x
Forces
uy 10- 3 10 T 3
--
10- 3 10- 3 10-3 10 - 3 10-3
-
-
--
--
-
10 -3
--
10- 3 10- 3 10-4
--
0.000000 0.000000 0.127500 x 0.150000 × 0.300000 x 0.150000 x 0.105000 x 0.105000 x 0.637500 x 0.112500 x 0.618750 × 0.121875 x 0.712500 ×
10 - 3 10-3 10-4 10 - 4 10-3 10 - 3 10-4 10 - 4 10-4 10- 3 10- 4
Fx
Fy
- 50.0 42.5 50.0 - 42.5 0,0 0.0 0.0 0.0 0.0 0.0 0,0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
the influence of discretization effects, it is preferable to investigate simple structures with regular meshes. For the elements dealt with here, it is easy to identify such structures: They are a square plate and a cube. Both structures are depicted in Fig. 3. Figure 3 also provides information on the dimensions and Poisson's ratio used in the calculations. It is assumed that both structures can vibrate freely. They are unrestrained . The plate can, of course, only vibrate in its plane. Its stresses obey plane stress conditions. To demonstrate the convergence properties of the new elements, the eigenfrequencies are computed with different meshes. The meshes are refined step-by-step. The solutions obtained with the elements OUAD5 and HEXA9 are compared with the eigenfrequencies rendered by the conventional consistent and diagonal mass matrices of the usual 4-node quadrilateral and 8-node hexahedron, respectively. These two matrices are based on the conventional bilinear interpolation functions. The diagonal mass matrix is constructed by summing up row-wise the terms of the consistent mass matrix as described above. All solutions are finally compared with the eigenfrequencies computed with fine finite element meshes. Exact analytical solu-
Table 2 Computed nodal displacements and forces for trxy = 100 Node
1 2 3 4 5 6 7 8 9 10 11 12 13
Coordinates
Displacements
Forces
X
y
Ux
Uy
F~
F~
0.0000 1.0000 1.0500 0.0000 0.2000 0.9000 0.8000 0.3000 0.5500 0.5250 0.9375 0.5375 0.1250
0.0000 0.0000 0.8500 1.0000 0.2000 0.1000 0.7000 0.7000 0.4250 0.0750 0.4125 0.8125 0.4750
0.000000 0.000000 × 0.552500 × 0.650000 × 0.130000× 0.650000 X 0.455000 × 0.455000 X 0.276250 × 0.487500 X 0.268125 X 0.528125 X 0.308750 X
0.000000 0.650000 × 0.682500 × 0.000000 0.130000× 0.585000 X 0.520000 X 0.195000 X 0.357500 x 0.341250 X 0.609375 X 0.349375 X 0.812500 X
- 50.0 - 52.5 50.0 52.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
- 50.0 42.5 50.0 - 42.5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
10 - 3 10- 3 10- 3 10 - 3 10-4 10- 3 10- 3 10- 3 10 - - 4 10-3 10-3 10- 3
10 - 3 10 - 3 10 - 3 10- 3 10- 3 10- 3 10-3 10 - 3 10-3 10-3 10-4
G. Sauer / Diagonal mass and associated stiffness matrices
44
/\
s=de length o f plate and cube plate t h i c k n e s s = 1
l
~ 2
/ \ / \
/ \ / \ /
\
\
\
\
\
\
\
\
\ / \ ,/\ / \ / \ /\ / \/ \ / \ / \ / / \/
\ / / /
/
\ / /
\
/ /
\ / /
\
\/
Fig. 3. Square plate with 4 x 4 and cube with 4 x 4 x 4 mesh.
P o i s s o n ' s ratio = 0 3
tions for the free in-plane vibrations of a plate and a freely vibrating cube do not exist. It is therefore assumed that the lowest frequencies of these structures are reproduced accurately by the fine meshes.
Quadrilateral Solutions obtained for the square plate are shown in Figs. 4-6. Before discussing the test results, an important aspect must be pointed out. In Figs. 4-6, the solutions obtained with models of different d.o.f, are confronted. It must be remembered that, when applying OUADS, the d.o.f, pertaining to the corner nodes are condensed out while assembling the stiffness matrix. Only the d.o.f, of the central nodes are retained. The final solution matrix of a N x N mesh merely comprises 2 N 2 d.o.f. A model made up of 4-node quadrilaterals, however, has
Eigenfrequency
Eigenfrequency
[rad/s]
700
[rad/s]
700 [~
exact
QUAD5
650 ,
\
~
4node
usual consist.
-~
4node
usual diagonal
650 "\-,...
600
600
550
550
500
500
. . . . Lt~
/
450
400
/
__
i 2
/
450
/
/
-----
/
QUADS
/
I i
i
i
{
{
3
4
5
6
7
M e s h size
_i
400 8
m {
2
exact
4node
usual consist
4node
usual diagonal
i
{
t
3
4
5
N
[NxN]
Fig. 4. Convergence of the first non-rigid-body eigenfrequency of the square plate.
M e s h size
k____
6
i
7
N
[NxN]
Fig. 5. Convergence of the second
non-rigid-body eigenfrequency of the square plate.
45
G. Sauer / Diagonal mass and associated stiffness matrices Error in Frequency
[%!
12
8
Error in Frequency
[%1
14 12
m
10
- B - 8nodeusualconsist.
8
-~-
HEXA9
8nodeUsualdiagonal
[ ~
[]
4 4 2
0
0
-
-
--
-4 -4 -6
-8
-8
,\
-10
&
-12 -12 -14 -16 4
5
6
7
8
9
10 11 12 13 14 15 16 17 18
Frequency [Number] Fig. 6. Test of OUADS; errors in the lowest 15 nonrigid-body frequencies of the plate using a 7 x 7 mesh.
7
8
9
10 11 12 13 14 15 16 17 18 19 20 21
Frequency [Number] Fig. 7. Test of HEXA9, errors in the lowest 15 nonrigid-body frequencies of the cube using a 4 × 4 × 4 mesh.
2 ( N + 1) 2 d.o.f. Taking into account this difference, one can properly judge the solutions originating from the use of OUADS. Figures 4 and 5 show the convergence of the two lowest non-rigid-body eigenfrequencies when successively refining the finite element model. Both figures demonstrate that, at the level of coarse meshes, OUAD5 produces solutions closer to the excact values than the other two matrices. The first non-rigid-body eigenfrequency is already approximated with good accuracy at the stage of a 3 × 3 mesh. Advancing to still finer meshes, the first eigenfrequency determined with the usual 4-node diagonal matrix catches up with the solution yielded by OUADS. These two solutions contrast with the eigenfrequency originating from the usual 4-node consistent mass matrix. Even in a 7 × 7 mesh, this mass matrix cannot reproduce the exact value. The eigenfrequency still exhibits an error of 2%. The second non-rigid-body eigenfrequency is best approximated by OUAD5 with all meshes. With a 7 × 7 mesh OUAD5 misses the exact eigenfrequency by 1.4%. The errors in frequency for the fifteen lowest non-rigid-body eigenfrequencies computed with a 7 × 7 mesh are depicted in Fig. 6. Out of the fifteen frequencies a total of eleven frequencies is best approximated by OUADS. For the first and the ninth eigenfrequencies the results of the usual 4-node diagonal mass matrix a r e m o r e accurate. The conventional 4-node consistent mass matrix renders superior solutions for the eigenfrequencies seven and eighteen.
Hexahedron The observations regarding the convergence characteristics of the OUAD5 element can, in principle, be transferred to the HEXA9 element. Typical results for this element are summarized in Table 3 and plotted in Fig. 7. Observing these results, it must again be r e m e m b e r e d that solutions obtained by finite element meshes with different d.o.f, are compared. Employ-
46
G. Sauer / Diagonal mass and associated stiffness matrices
Table 3 Lowest four non-rigid body frequencies of the cube computed with different meshes; multiple frequencies are skipped a Mesh
Mass matrix
oJ 7
co 9
t012
2x 2x 2
HEXA9 8-node usual consistent 8-node usual diagonal
1.0181
1.1363
1.1689
1.1839
1.0742 0.7161
1.3994 0.8523
1.5191 0.8771
1.5191
HEXA9 8-node usual consistent 8-node usual diagonal
O.9606
1. 2225
1. 2405
1.2580
0.9743 0.7978
1.3197 1.0133
1.3517 1.0527
1.4656
HEXA9
O.9308
1.2264
1.2336
1.3075
8-node usual consistent 8-node usual diagonal
0.9375 0.8338
1.2722 1.0816
1.2937 1.1573
1.4353
(/.8945
1.2071
1.2266
1.3890
3x 3x 3
4X4X4
10 x 10 x 10
o315
(I.9162
1.1360
1.2046
a coi = ~oi / EV~-~ ; the closest approximations to the exact values in italics. ing HEXA9, a finite e l e m e n t m e s h o f a N x N x N subdivision c o n d e n s e s to a system o f 3 N 3 e q u a t i o n s which, for l a r g e r N, signifies a c o n s i d e r a b l e r e d u c t i o n in the n u m b e r o f equations. T h e usual 8 - n 0 d e e l e m e n t s b u i l d u p a system o f 3 ( N + 1) 3 equations. HEXA9 c o u l d also be very a d v a n t a g e o u s f r o m t h e p o i n t o f view o f s m a l l e r s o l u t i o n m a t r i c e s , p r o v i d e d it yields accurate eigensolutions. O n i n s p e c t i n g T a b l e 3, o n e n o t e s t h a t this p r e r e q u i s t e is fulfilled. Surprisingly g o o d results a r e o b t a i n e d for t h e lowest frequencies. T h e lowest f r e q u e n c i e s u p to the 14th f r e q u e n c y a r e m o s t closely a p p r o x i m a t e d by HEXA9 with all meshes. F i g u r e 7 d e m o n s t r a t e s that only two e i g e n f r e q u e n c i e s a r e b e t t e r r e p r o d u c e d by o n e of t h e o t h e r two matrices. T h e s e two i d e n t i c a l f r e q u e n c i e s a r e c o m p u t e d m o r e a c c u r a t e l y by t h e u s u a l 8 - n o d e c o n s i s t e n t m a s s matrix. T h e solutions f o u n d by using t h e usual 8 - n o d e d i a g o n a l mass m a t r i x a r e less exact t h a n t h o s e p r o d u c e d by HEXA9 for t h e w h o l e f r e q u e n c y r a n g e c o n s i d e r e d .
Conclusion Q u a d r i l a t e r a l a n d h e x a h e d r o n e l e m e n t s with c e n t r a l n o d e s can be u s e d effectively to c o m p u t e e i g e n f r e q u e n c i e s a n d m o d e s h a p e s o f l i n e a r elastic structures. By c h o o s i n g a p p r o p r i ate i n t e r p o l a t i o n f u n c t i o n s a n d an a s s o c i a t e d i n t e g r a t i o n scheme, d i a g o n a l mass m a t r i c e s a r e g e n e r a t e d . Since t h e y g a t h e r t h e m a s s on t h e c e n t r a l n o d e , t h e d.o.f, o f t h e massless c o r n e r n o d e s can b e c o n d e n s e d o u t f r o m t h e stiffness matrix. T h e size of t h e solutions m a t r i c e s is r e d u c e d a n d now d e p e n d s on t h e n u m b e r o f e l e m e n t s r a t h e r t h a n on t h e n u m b e r o f nodes. T h e c o n v e r g e n c e c h a r a c t e r i s t i c s o f t h e s u g g e s t e d e l e m e n t s a r e satisfactory. T h e s e e l e m e n t s a p p r o x i m a t e t h e exact solutions m o r e closely t h a n e l e m e n t s t h a t e m p l o y t h e usual c o n s i s t e n t o r t h e u s u a l d i a g o n a l m a s s matrices. T h e m a t r i c e s b a s e d on t h e e x p a n d e d i n t e r p o l a t i o n f u n c t i o n s a r e s u i t e d to m o d e l l i n g t h e d y n a m i c s t r u c t u r a l p r o p e r t i e s .
References [1] R.W. C L O U G H , "The finite element in plane stress analysis", Proc. 2nd A S C E Conf. on Electronic Computation, Pittsburgh, PA, September 1960. [2] J.S. ARCHER, "Consistent mass matrix for distributed mass systems", J. Struct. Div. A S C E 89 (ST 4), pp. 161-178, 1963.
G. Sauer / Diagonal mass and associated stiffness matrices
47
[3] J.S. PRZEMIENIECKI, Theory of Matrix Structural Analysis, McGraw-Hill, New York, 1968. [4] K.K. GUPTA, "Numerical formulation for a higher order plane finite dynamic element", Int. J. Numer. Methods Eng. 20, pp. 1407-1414, 1984. [5] K.K. GUPTA, "Development of a solid hexahedron finite dynamic element", Int. J. Numer. Methods Eng. 20, pp. 2143-2150, 1984. [6] G. SAUER and M. WOLF, " A modified beam element mass matrix", Z. Angew. Math. Mech. 68, pp. 483-490, 1988. [7] K.J. BATHE and E.L. WILSON, Numerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1976. [8] S.W. KEY and Z.E. BEISINGER,"Slade D: A computer program for the dynamic analysis of thin shells", Report SLA-73-0079, Sandia Laboratories, Albuqerque, NM, 1973. [9] E. HINTON, T. ROCK and O.C. ZIENKmWlCZ, " A note on mass lumping and related processes in the finite element method", Int. J. Earthquake Eng. Struct. Dyn. 4, pp. 245-249, 1976. [10] H.L. SCHREYER,"Consistent diagonal mass matrices and finite element equations for one-dimensional problems", Int. J. Numer. Methods Eng. 12, pp. 1171-1184, 1978. [11] G. SAUER,"Consistent diagonal mass matrices for the isoparametric 4-node quadrilateral and 8-node hexahedron elements", Commun. Numer. Methods Eng. 9, pp. 35-43, 1993. [12] R.L. TAYLOR,J.C. SIMO,O.C. ZIENKIEWICZ and A.C.H. CrtAN, "The patch t e s t - - A condition for assessing fem convergence", Int. J. Numer. Methods Eng. 22, pp. 39-62, 1986. [13] A. RAZZAOtJE, "The patch test for elements", Int. J. Numer. Methods Eng. 22, pp. 63-71, 1986.