CHAPTER 5: FINITE ELEMENT VIBRATION ANALYSIS

Introduction

In previous topics we learned how to model the dynamic behavior of multi-DOF systems,
as well as systems possessing infinite numbers of DOF. As the reader may realize, our
discussion was limited to rather simple geometries and boundary conditions, mainly
because more complicated geometries tend to have progressively more complicated exact
solutions, or even no exact solutions at all. The vast majority of geometries which occur in
practical engineering applications are complicated three-dimensional continua which
cannot be represented adequately by the simple closed-form mathematical models. These
theoretical models also cannot portray the appreciable non-linear or anisotropic material
characteristics which are often met with in practice. In such situations therefore, we must
resort to numerical methods in which the continuum of infinitesimal material particles is
represented by an approximately equivalent assembly of inter-connected discrete elements
which are each so simple that they can be treated individually as mathematical continua.
There are a number of methods whereby such models can be analyzed such as numerical
solution of differential equations, finite differences, finite elements, boundary elements,
relaxation techniques, and so on. In this topic, we will demonstrate the Finite Element
Method (FEM) as a typical powerful approach which can handle vibration analysis.

In essence, the FE technique is a numerical technique in which a continuous elastic

structure, or continuum, is divided (discretized) into small but finite substructures, known
as elements. Elements are interconnected at nodes. In this way, a continuum with infinite
number of degrees of freedom can be modeled with a set of elements having a finite number
of degrees of freedom. It is noted that while each finite element represents a continuous
system by itself possessing infinite number of DOF, we can choose the size of the element
to be small enough, so that the deformation within the finite element can be approximated
(interpolated) by relatively low-order polynomials.

1
In this topic, we will discuss the application of FE procedures in vibration analysis, with
emphasis on developing equations of motion for elementary geometries, developing FE
codes on MATLAB and applying commercial FE software to handle more advanced
geometries.

Finite Element Analysis of Rods

Consider an elastic uniform rod of total length LR as shown in Fig.1. Upon applying the FE
technique, the rod is discretzed into a finite number of elements. As the rod under
investigation is uniform, it is assumed that all elements used to mesh the complete rod are
identical. In more advanced problems, of course, this technique may not be possible, as the
entity to be meshed often has a more complicated contour or geometry, which makes it
difficult to employ identical elements. The technique described herein, however, can be
applied to non-identical elements.

LR

1 2 3 4
x
i j
L

Figure 1. Uniform rod.

Let u(x) denote the axial displacement within an element. It is recognized that
displacements are also functions of time, but for the sake of simplicity in writing the
equations, the time t is dropped, but the reader must remember that time dependence is
indeed retained in the upcoming analysis. Note also that the rod element by itself is treated
as a continuous structural member, meaning that axial displacement within the element is
a continuous function. An essence of the finite element technique is that field variables, in
this case axial displacements, within each element is interpolated (or approximated) in
terms of the nodal quantities, that is the displacements measured at the two nodes bounding
the element. In this way, a whole structure is regarded as a multi-DOF system, with field
variables evaluated only the nodes. Once these quantities are known, the displacement field

2
within the entire structure can be evaluated using the interpolation functions initially
assumed for the structure. Several interpolation functions can be chosen to properly satisfy
the physics of the problem.

For the problem at hand, the displacement field within the rod element is approximated by
a linear function:
u ( x)  1   2 x , 0  x  L (1)
This can be written as:
 
u ( x)  1 x  1  (2)
 2 
In this way, the vector of nodal degrees of freedom is expressed as:
 ui 
   u 
e
(3)
 j
But we have:
ui  1
(4)
u j  1   2 L

Therefore:
 ui  1 0   1   1 
   u   1
e
     A  
    (5)
 j  L   2   2

from which:
 1 
    A   
1 e
(6)

 2

Performing the matrix inverse yields:

 1   1 0   ui 
    (7)
 2    1 L 1 L  u j 
Substituting Eq. (7) into (2) gives:
 1 0   ui 
u ( x)  1 x    (8)
  1 L 1 L  u j 

3
or:

 x x   ui 
u ( x)  1    (9)
 L L  u j 

which is written as:

u ( x)   N ( x)  e  (10)

The vector  N ( x) is called the vector of interpolation functions or shape functions, and

is typical in all finite element procedures. Once again, using this vector, we can express the
displacement at any point within the rod element in terms of the nodal degrees of freedom.

Several methods can be employed to obtain the equation of motion of rod elements. In this
section we will use Lagrange’s equation. In this regard, we need to express the strain energy
associated with the rod element. The reader is referred to any elementary textbook in
strength of materials to express the strain energy as:
2
 u 
L
1
U  EA   dx (11)
2 0
x 

Inserting the shape functions, as determined in (10) and noting that (see example in box
below):

 u 
2

      N x   N x  
e T T e
(12)
 x 
we obtain:
L
U
2
  EA  N x   N x dx  e 
1 e T T
(13)
0

where
N ( x) 
 N x     (14)
 x 
The strain energy can then be written as:

U
2
   K e   e 
1 e T
(15)

4
where  K e  is the element stiffness matrix, given by:
2 2

L
 K e   EA  N x   N x dx
T
(16)
0

Performing the above calculation yields:

EA  1 1
 K e  
L  1 1 
(17)

which is the element stiffness matrix for a uniform rod element.

Example:
5
Let a  2 3 , b   
8
5
Then u  ab  2 3    34
8
Squaring the above equation is expressed as:
2 5
u 2  b a ab  5 8  2 3     34  34 
T T

3 8

The element kinetic energy is then evaluated for the rod element and can be expressed as:
L
1
T   A u 2 dx (18)
2 0

Inserting the shape functions, as determined in Eq. (10) yields:

L
1
   A  N   N dx  e 
T
T  e
T
(19)
2 0

which can also be written as:

1 e
   
T
T   M e   e (20)
2

where  M e  is the element mass matrix, given by:

2 2

L
 M e    A  N   N dx
T
(21)
0

5
Performing the above calculation yields:
1 3 1 6 
 M e    AL   (22)
1 6 1 3 
which is the element mass matrix for a uniform rod element. Now we can apply
Lagrange’s equation, which is in the form:
d  T  T V
   0 (23)
dt  q  q q
or more precisely:

d  T  T V
   0 (24)
  
dt    e       e 

e

The resulting equations of motion are therefore:

 
 M e   e   K e   e   0 (25)

which is the element equation of motion for free vibration. Upon assembly of the element
equations of motion (see next section), we can determine the equations of motion for the
entire structure in the form:

 M     K    0 (26)

We can easily extend the derivation to include a damping term and a forcing term in the
equation of motion:

 M    C     K    f (t ) (27)

Assembly of Element Matrices

Assembly of the element matrices can be accomplished by the routine assembly process
used in finite element procedures. As an example, consider a rod consisting of 2 elements
as shown in Fig.2.

1 2
1 2 3

Figure 2. Rod with 2 elements.

6
The strain energy of the rod is equal to the sum of strain energy in each element, i.e.:

   K e1   e1   e 2   K e 2   e 2 
1 e1 T 1 T
U (28)
2 2
This equation can be written as:
1 EA  1 1  u1  1 EA  1 1 u2 
U  u1 u2     u2 u3   
L  1 1  u2  2 L  1 1  u3 
(29)
2

For simplicity we assumed both elements have the same properties, but this can be easily
changed if otherwise. Now the strain energy for the whole structure can be written as:
1
U     K  
T

2
 1 1 0   u1  (30)
1    
 u1 u2 u3   1 1  1 1 u2 
2
 0 1 1  u3 

So the assembled stiffness of the whole structure is:

 1 1 0 
 K    1 2 1 (31)
 0 1 1 

The assembly process works as follows: start with a zero matrix having the size of the
total number of degrees of freedom (3 in the example above), then loop over each element
populating the overall stiffness matrix by inserting the entries in their appropriate
locations. In the above example, the assembled stiffness matrix will be a banded matrix,
with entries close to the diagonal and zero elsewhere. The same procedure is used to
assemble the mass matrix (using expressions for the kinetic energy), and for any number
of elements. Elements having different material properties or geometrical parameters (cross
sectional areas) can easily be assembled. For example, if a rod is discretized into 6 elements
and 7 nodes, we have:
1 2 3 4 5 6 7

1 2 3 4 5 6

The assembly process works by inserting the stiffness matrix of each element (bounded by
the proper nodes) until the overall stiffness matrix is built. The same applies for the mass
matrix.

7
 1 2 3 4 5 6 7

1 1 -1 0 0 0 0 0

2 -1 1+1 -1 0 0 0 0

3 0 -1 1+1 -1 0 0 0

EA
K   4
0 0 -1 1+1 -1 0 0

L 5
0 0 0 -1 1+1 -1 0

6
0 0 0 0 -1 1+1 -1

7
0 0 0 0 0 -1 1

Boundary Conditions

Inspection of the overall stiffness matrix indicated in the previous section reveals that it is
essentially singular, meaning that upon application of an external force one cannot solve
for the unknown displacements (in static analysis, for instance) using  K     F  by

multiplying by the inverse of the stiffness matrix, since it is non-existent. That is so because
the structure in this case is not restrained, and hence application of an external force causes
instability in the analysis. Such an error is often encountered in commercial FE software,
where the user is alerted of a “singular stiffness matrix” error, which is an indication that
the structure has insufficient constraints.

Application of boundary conditions can be accomplished using the same procedures

typically adopted in static FE analysis. As an illustrative example, consider the rod shown
in Fig.3, where one end is rigidly fixed to a wall, while the other is subjected to a force.

1 2 F
1 2 3

Figure 3. Rod with one end fixed, other end loaded.

8
In this case, we have:
 F1   1 1 0   u1 
  EA   
 F2    1 2 1 u2  (32)
 F  L  0 1 1  u 
 3   3
However, F1 is an unknown reaction force that is yet to be determined. Although we can
directly obtain the value of this reaction from physical intuition (equal and opposite to the
applied tip load), we cannot make this intervention before using the FE analysis to solve
the problem! Moreover, in more complicated cases, such as statically indeterminate
problems, we cannot obtain the reaction forces using only the equilibrium equations. As
far as displacements are concerned, we will also regard u2 and u3 as unknown
displacements. Hence equation (32) represents three equations in three unknowns. In order
to proceed with the solution, we will first eliminate the equations pertaining to unknown
reaction forces, in this case the first equation. This is accomplished by removing the first
element in both the load vector and displacement vector, as well as removing the first row
and column of the stiffness matrix. To this end, we will get:
 F2  EA  2 1 u2 
     (33)
 F3  L  1 1  u3 

u 2 
which can be solved for the unknown displacements   as the vector of externally
 u3 

F 
applied loads  2  is known. Once the displacements are known, we can then substitute
 F3 
into equation (32) with:
 F1   1 1 0   0 
  EA   u 
 F2    1 2 1 2 (34)
 F  L  0 1 1  u 
 3   3
in order to obtain the unknown reaction force. For a cantilever rod with 6 elements:
1 2 3 4 5 6 7

1 2 3 4 5 6

9
The reduced stiffness matrix becomes 6  6 and is given by:
1 2 3 4 5 6 7

1 1 -1 0 0 0 0 0

2 -1 2 -1 0 0 0 0

3 0 -1 2 -1 0 0 0

EA
K   4
0 0 -1 2 -1 0 0

L 5
0 0 0 -1 2 -1 0

6
0 0 0 0 -1 2 -1

7
0 0 0 0 0 -1 1

The same applies for the mass matrix, and all calculations are done on the reduced matrices.

Vibration Analysis of Rods

Let us now study the free and forced vibration analysis of rods using the finite element
procedure outline in the previous sections. For free vibration, we have:

 M     K    0 (35)

Imposing harmonic motion in the form:

    sin t (36)

we get:
 K    M  2  (37)

hence:

 M   K    2 
1
(38)

which is clearly an eigenvalue problem that can be solved to obtain the natural frequencies
and mode shapes of the rod under investigation. It is perhaps worthy to mention at this
point that such an analysis can be readily made for any boundary conditions imposed on

10
the structure, which is in clear contrast with exact analytical solutions that largely depend
on the boundary conditions.

For steady-state harmonic analysis, on the other hand, we have:

 M     K     f (t ) (39)

with:
 f (t )  F  sin t (40)

The resulting displacement, described by equation (28) is inserted into (32) yielding:
 K    2  M    F  (41)

which can be solved for the steady-state displacement amplitude at each excitation
frequency:
1
   K    2  M  F  (42)

where the matrix on the right hand side is commonly known as the dynamic stiffness
matrix.

As an illustrative example, consider a fixed-free rod with the following properties:

Modulus of elasticity 80 GPa
Cross-sectional area 0.01 m2
Length 8m
Density 7800 kg/m3

Table 1 lists the first three natural frequencies as we increase the number of elements used,
together with the analytical solution. Convergence is seen to take place as we increase the
number of elements used for meshing the rod.

11
 Table 1. Natural frequencies and comparison with analytical solution.
Number of elements [Hz] [Hz] [Hz]
1 110.4 --- ---
2 102.7 358.7 ---
3 101.2 331.1 600.6
4 100.7 317.7 577.2
5 100.5 311.4 551.8
6 100.4 308 536.4
10 100.2 303 513.3
20 100.1 300.9 503.6
40 100.1 300.4 501.2
Analytical Solution 100.0801 300.2 500.4

Figure 4 shows the first three mode shapes, as predicted by the present FE analysis.

Figure 4. First three mode shapes of a fixed-free rod.

Consider now the case when the rod is subjected to a harmonic force acting at its tip and
having a magnitude of 100N. We can determine the steady-state displacement amplitude
at any point on the rod in accordance with the technique outlined previously. Figure 5
shows such frequency response, where the resonance is shown to occur at the natural
frequencies indicated.

12
 Figure 5. Frequency response of harmonically excited rod.

Damping can easily be implemented in the present model. Figure 6 shows the frequency
response with damping included in the form of Rayleigh (proportional) damping.

Figure 6. Frequency response of harmonically excited rod with damping.

13
Finite Element Analysis of Beams
The finite element analysis of beams follows the same procedure as that described for rods,
except for the energy expressions, shape functions and element degrees of freedom. Figure
7 shows a beam in flexural vibrations, where the transverse displacement is expressed as
w( x) . Once again, displacements are in fact functions of time, but the time t is dropped for
brevity in the following analysis.

w(x)
x
L
Figure 7. Beam in flexural vibrations.

For beams, the displacement field within the rod element is approximated by a cubic
polynomial:
w( x)  1   2 x   3 x 2   4 x3 , 0  x  L (43)
The choice of this function is based upon relations from strength of materials governing
the Euler-Bernoulli beam theory, in which an element loaded by concentrated forces at its
nodes features a linear variation in bending moment. Since bending moment is a function
of the second derivative of transverse displacement, the above polynomial is justified.
Equation (43) can be written as:
1 
 
3  2
w( x)  1 x x2 x   (44)
 3 
 4 

The vector of nodal degrees of freedom is then expressed as:

 wi 
 
   wi 
e
(45)
 j
 j 

14
where
w
 (46)
x
is the slope. Substituting to get the values of the element degrees of freedom at each node,
we have:
 wi  1 
   
 i  2 
    A    (47)
w j   3 
 j   4 

from which:
1   wi 
   
 2 1   i 
    A    (48)
 3  w j 
 4   j 

Substituting (48) into (44) gives:

 wi 
 
1   i 
w( x)  1 x x2 x 3   A    (49)
w j 
 j 

which is written as:

w( x)   N ( x) e  (50)

Once again, the vector  N ( x) is the vector of interpolation functions or shape functions

for beam elements, which is used to express the displacement at any point within the rod
element in terms of the nodal degrees of freedom.

The strain energy of the beam element is expressed as:

2
1
L
 2w 
U  EI   2  dx (51)
2 0  x 

Inserting the shape functions, as determined in (50) yields:

15
 L
U   e  EI   N xx   N xx dx  e 
1 T T
(52)
2 0

where
  2 N ( x) 
 N xx     (53)
 x 
2

The strain energy of the beam element can be written as:

U
2
   K e   e 
1 e T
(54)

where  K e  is the element stiffness matrix, given by:

4 4

L
 K e   EI   N xx   N xx dx
T
(55)
0

Performing the above calculation yields:

12 L 6 12 L 6 
 4L 6 2 L 
EI 
 K   2
e
(56)
L  12 L 6 
 
 4L
which is the element stiffness matrix for a uniform beam element. The element kinetic
energy of a beam element is expressed as:
L
1
T  A w2 dx (57)
2 0

Inserting the shape functions, as determined in (50) yields:

L
1
   A  N   N dx  e 
T
T  e
T
(58)
2 0

which can also be written as:

1 e
   
T
T   M e   e (59)
2

where  M e  is the element mass matrix, given by:

4 4

L
 M    A  N   N dx
e T
(60)
0

16
Performing the above calculation yields:
13 35 11L 210 9 70 13L 420 
 L 105 13L 420  L2 140 
2
 M e    AL  (61)
 13 35 11L 210 
 
 L2 105 
which is the element mass matrix for a uniform beam element. Applying Lagrange’s
equation yields the element equation of motion for free vibration as:

 
 M e   e   K e   e   0 (62)

Upon assembly of the element equations of, we can determine the equations of motion for
the entire structure in the form:

 M     K    0 (63)

Or more generally:

 M    C     K    f (t ) (64)

5.7 Example

Consider a cantilever beam with the following properties:

Modulus of elasticity 80 GPa
Cross-section (width x height) 20 mm × 0.4 mm
Length 0.2 m
Density 2700 kg/m3

Table 2 lists the first three natural frequencies as we increase the number of elements used,
together with the analytical solution. Once again, convergence is seen to take place as we
increase the number of elements used for meshing the rod.

17
 Table 2. Natural frequencies and comparison with analytical solution.
Number of elements [Hz] [Hz] [Hz]
2 8.7974 55.5732 187.9587
3 8.7940 55.2866 156.2198
4 8.7934 55.1698 155.4919
5 8.7932 55.1331 154.8515
6 8.7932 55.1192 154.5796
10 8.7931 55.1074 154.3365
Analytical Solution 8.7922 55.1034 154.3461

Figure 8 shows the first three mode shapes, as predicted by the present FE analysis.

Figure 8. Cantilever beam mode shapes.

18
Finite Element Analysis of Strings
This section presents the finite element model of a stretched string undergoing transverse
vibration. Figure 9 shows a schematic diagram of a string stretched between two rigid ends.

Figure 9. Taught string.

Consider an element of the stretched string, as indicated in Fig. 10. The string is displaced
by a transverse displacement dy The tension T in the foil is assumed to be constant
assuming the string displacement is comparatively smaller than its length.

dy
dx
y
Figure 10. String element – deformed and undeformed.

The kinetic energy of the string element is given by:

L
T   y dx
1 2
2 (65)
0

The elongation of the string element due to the transverse deformation is given by:
1

   dx 2  dy 2   dx
2
(66)

For small displacements, the elongation reduces to:

 1  dy 2 
1

 dy 2  2  1 dy 2 
  dx 1  2   dx  dx 1  2
 dx  dx     (67)
 dx   2 dx   2  dx  

Hence the strain is given by:

2
1  dy 
   (68)
2  dx 

In this way, the strain energy of the string is given by:

19
 1
L
 dy 2 
2 0  dx  
V T   dx (69)

The finite element model developed in this work relies upon discretizing the string domain
into two-noded elements with one degree of freedom per node denoting transverse
displacement. The transverse displacement at any point within an element is interpolated
by:
y  1   2 x (70)

Accordingly, the transverse displacement at any point in the string element can be
expressed in terms of the nodal degrees of freedom as:
y   N y  e  (71)

where  N y  is the vectors of shape functions and  e  is the nodal degrees of freedom
vector. Imposing the shape functions yields:
L
1 e T
    N y   N y dx  e  
T
T  
e
(72)
2 0

Accordingly, the mass matrix of the foil element is given by:

L
 M e      N y  N dx
T
y (73)
0

Performing this integral yields:

1 3 1 6 
 M e    L   (74)
1 6 1 3 
Similarly, the strain energy of the string element can be expressed as:
L   dN y T  dN y  
 K    T  
e
   dx (75)
0   dx   dx  

Performing this integral yields:

T  1 1
 K e   
L  1 1 
(76)

20
Finite element analysis of trusses

A truss is made up of a number of links. Each link is connected by two hinges at its ends
and carries either tension or compression, as illustrated in Figure 1. Each link can therefore
be treated as a rod under uniaxial loading. The one-dimensional rod elements developed
earlier, however, are not suitable to model a truss like the one shown below, since different
links are inclined at different angles, and every link undergoes two-dimensional motion.
Hence it is necessary to develop a two-dimensional truss element. We can use finite
element analysis to model such trusses, not only to calculate how much they deform under
static load, but also to study their dynamics, such as natural frequencies and dynamic
response.

Figure 1. Truss under load.

Truss links can be modeled as rods since each member carries either tension or
compression. Each member has two connection points (hinges) which are also called
nodes. Let us consider a truss link inclined at an angle , as illustrated in Figure 2. The
local degrees of freedom are u1 and u2, which designate the axial displacement of each
node. The global degrees of freedom are designated U1, V1, U2 and V2, and represent the
displacements of the nodes in the global Cartesian coordinates. In order to develop the
mass and stiffness matrices of such an element, we shall assume that displacements are
small with respect to the element length, hence the inclination angles do not change
significantly due to the applied loads.

21
 V Deformed
v
v1 u

V1 Undeformed
u1


Node 1
U1
U
Figure 2. Truss element

From Figure 2, we can express the local displacements in terms of the global ones as:
u1  U1 cos   V1 sin 
(1)
u2  U 2 cos   V2 sin 
and
v1  V1 cos   U1 sin 
(2)
v2  V2 cos   U 2 sin 
Equation (1) can be written in a matrix form as:
U1 
 u1  cos  sin  0 0   V1 
   
sin   U 2 
(3)
u2   0 0 cos 
 V2 

Similarly, Equation (2) can be written as:

U1 
 v1    sin  cos  0 0   V1 
   
cos   U 2 
(4)
v2   0 0  sin 
 V2 

Recall that:
 1 0   u1 
u ( x)  1 x    (5)
  1 L 1 L  u2 
which means a linear interpolation function can be used to obtain the displacement u at any
given x as a function of the nodal displacements u1 and u2. The same linear interpolation
function can be assumed for the transverse displacement v. Hence we can write:

22
  1 0   v1 
v( x)  1 x    (6)
 1 L 1 L  v2 
Combining Equations (3) and (5) yields:
U1 
 1 0  cos  sin  0 0   V1 
u ( x)  1 x    
 1 L 1 L   0 0 cos  sin   U 2 
 V2 
  x  x x cos  x sin   e
 cos  1   sin  1      (7)
  L  L L L 
  N ( x)  e 

Note that the vector  N ( x) is the vector of interpolation functions for axial displacements

u(x) a truss element of length L, inclined at an angle , and the vector of nodal degrees of
freedom is now given by:
U1 
V 
   U1 
e
(8)
 2
 V2 

Similarly, combining Equations (4) and (6) yields:

 1 0   v1 
v( x)  1 x   
  1 L 1 L  v2 
U1 
 1 0    sin  cos  0 0   V1 
 1 x    
 1 L 1 L   0 0  sin  cos   U 2  (9)
 V2 
 x  x   x sin  x cos   e
 sin    1  cos    1   
 L  L  L L 
  N v ( x)  e 

The vector  N v ( x) is the vector of interpolation functions for transverse displacements

v(x) for a truss element. The strain energy of the truss element is identical to the rod, and
is given by:

23
 2
 u 
L
1
U  EA   dx (10)
2 0
x 

Which reduces to:

L
U   e  EA  N x   N x dx  e 
1 T T
(11)
2 0

Where:
N ( x) 
 N x     (12)
 x 
The strain energy can then be written as:

U
2
   K e   e 
1 e T
(13)

where  K e  is the stiffness matrix of the truss element, which is given by:
4 4

L
 K e   EA  N x   N x dx
T
(14)
0

Performing the above integration yields:

 cos 2  Sym 
 
EA  cos  sin  sin 2

 K e   (15)
L   cos 2   cos  sin  cos 2  
 
  cos  sin   sin 2  cos  sin  sin 2  

which is the element stiffness matrix for a uniform truss element inclined at an angle .
In order to calculate the kinetic energy of the element, we have to observe that the
transverse motion contributes to the kinetic energy, unlike the case of a one-dimensional
rod, so the kinetic energy can be expressed as:
L
 A  u 2  v 2 dx
1
T (16)
2 0

Inserting the interpolation functions given by Equations (7) and (9) yields:
L L

 
1 e
 A  N   N dx  e    e   A  N v   N v dx  e 
1
T T
T 
T T
(17)
2 0
2 0

From which we obtain the mass matrix as:

L L
 M e  T  A  N   N dx   A  N v   N v dx
T T
(18)
0 0

24
Performing these integrals yields:
1 3 0 1 6 0 
 0 1 3 0 1 6
 M    AL 
e  (19)
1 6 0 1 3 0 
 
 0 1 6 0 1 3
Which is the mass matrix of the truss element. Note that it does not depend on the
element orientation.

25