Physical Intuition and Continuum Dynamics Context

In the previous chapters, we focused entirely on static finite element analysis, where we assumed that external loads were applied so slowly that the structural acceleration and velocity terms were negligible. In static systems, the structure is in perfect equilibrium at every instant: the internal elastic restoring forces exactly balance the external loads. However, real-world structural systems are frequently subjected to transient, cyclic, or sudden dynamic excitations. Examples include:

When loads vary rapidly with time, the structure does not remain in static equilibrium. The mass of the structure resists acceleration, giving rise to inertial forces. Additionally, internal friction, material hysteresis, and joints dissipate energy, giving rise to damping forces. To model these physical phenomena, we must transition from static equilibrium equations to the dynamic equations of motion.

For a continuous, three-dimensional solid body, the local dynamic equilibrium is governed by the Cauchy equations of motion:

where:

To solve these continuous partial differential equations (PDEs) for complex geometries, we apply the finite element discretization process. This converts the continuous field variables into discrete nodal displacements. The resulting system-level system of ordinary differential equations (ODEs) is written as:

where:

Core Concepts

1. The Physical Meaning of Mass Matrices

The global mass matrix $\mathbf{M}$ is a mathematical representation of how the structure's mass is distributed among its nodal degrees of freedom. Just as the stiffness matrix $\mathbf{K}$ maps nodal displacements to elastic nodal forces, the mass matrix $\mathbf{M}$ maps nodal accelerations to inertial nodal forces:

In structural dynamics, the mass matrix must satisfy several fundamental physical properties:

Thus, $\mathbf{M}$ is a symmetric positive-definite (SPD) matrix.

2. Lumped vs. Consistent Mass Formulations

There are two primary methods used in finite element analysis to formulate the mass matrix: the Lumped Mass Matrix and the Consistent Mass Matrix.

A. Lumped Mass Matrix

The lumped mass matrix method is a simplified approach where the continuous mass of the element is divided and mathematically "lumped" directly onto its nodes. All off-diagonal terms are assumed to be zero, resulting in a diagonal mass matrix:

For example, for a 1D uniform bar element of length $L$, cross-sectional area $A$, and density $\rho$, the total mass is $m = \rho A L$. In a lumped mass formulation, half of this mass is placed at local Node 1, and the other half is placed at local Node 2:

*Advantages of Lumped Mass matrices*:

If the mass matrix is diagonal, no system of equations needs to be solved at each time step, making explicit solvers incredibly fast.

B. Consistent Mass Matrix

The consistent mass matrix is derived using the same variational principles and interpolation shape functions that were used to formulate the element stiffness matrix. The term "consistent" means that the mass distribution is mathematically consistent with the assumed displacement field within the element.

For an element, the kinetic energy $T^e$ is given by:

Within the element, the continuous velocity field $\dot{\mathbf{u}}(x,y,z,t)$ is interpolated from the nodal velocities $\dot{\mathbf{d}}^e(t)$ using the shape functions $\mathbf{N}$:

Substituting this interpolation into the kinetic energy equation yields:

Since the kinetic energy of a discrete system is expressed as $T^e = \frac{1}{2} \dot{\mathbf{d}}^{e\,T} \mathbf{m}^e \dot{\mathbf{d}}^e$, we identify the term in parentheses as the element consistent mass matrix:

*Advantages of Consistent Mass matrices*:

Step-by-Step Derivations

1. Consistent Mass Matrix for a 1D Linear Bar Element

Let us derive the consistent mass matrix for a uniform, linear 1D bar element of length $L$, cross-sectional area $A$, and mass density $\rho$.

The element has two nodes. The coordinate $x$ ranges from $0$ to $L$. The linear shape functions are:

The vector of shape functions is:

Using the consistent mass matrix formula:

Assuming $\rho$ and $A$ are constant over the element:

Let us evaluate each integral component step-by-step:

Let $u = 1 - \frac{x}{L}$, then $du = -\frac{1}{L} dx \implies dx = -L \, du$. Changing limits: when $x=0 \implies u=1$; when $x=L \implies u=0$.

Assembling these components back into the matrix:

Compare this with the lumped mass matrix:

Notice that the sum of the rows in both formulations equals the total mass of the element, $\rho A L$, confirming conservation of mass.

2. Consistent Mass Matrix for an Euler-Bernoulli Beam Element

For a 2-node Euler-Bernoulli beam element with four degrees of freedom (deflection $v$ and slope $\theta$ at each node), we write the shape functions using Hermite cubic polynomials:

The vector of shape functions is:

The consistent mass matrix is:

For a uniform beam, this requires integrating 16 terms of products of cubic polynomials. Let us evaluate one off-diagonal term, $m_{12}$:

Expanding the product inside the integral:

Now, integrating from $0$ to $L$:

Finding a common denominator ($3 \times 5 \times 6 \times 7 \to 210$):

Multiplying by $\rho A$:

By evaluating all 16 integrations, we obtain the full $4 \times 4$ element consistent mass matrix for a flexural Euler-Bernoulli beam:

Solved Example: Natural Frequencies and Modal Vectors of a Cantilever Bar

Problem Statement: Consider a uniform elastic bar of length $L$, cross-sectional area $A$, Young's Modulus $E$, and mass density $\rho$. The bar is fixed at its left end ($x=0$) and free at its right end ($x=L$).

Using a two-element model (which contains 3 nodes), calculate the natural frequencies and eigenvectors (mode shapes) using the Consistent Mass Matrix. Compare the results with the exact analytical frequencies.

Step 1: Discretize the domain

We divide the bar of length $L$ into two equal elements, each of length $L_e = L/2$.

Thus, the system has 2 active degrees of freedom: $\mathbf{u} = \begin{bmatrix} u_2 & u_3 \end{bmatrix}^T$.

Step 2: Formulate element matrices

The stiffness matrix for a 1D bar element of length $L_e$ is:

Since $L_e = L/2$, we have:

The consistent mass matrix for a 1D bar element of length $L_e$ is:

Step 3: Assemble the global stiffness and mass matrices

Let us assemble the matrices by mapping local nodes to global nodes:

Global Stiffness Matrix $\mathbf{K}_{global}$ (size $3 \times 3$ before boundary conditions):

Global Mass Matrix $\mathbf{M}_{global}$ (size $3 \times 3$ before boundary conditions):

Step 4: Apply boundary conditions

At Node 1, the bar is fixed, so $u_1 = 0$. We eliminate the first row and column of our global matrices to solve for the active degrees of freedom $u_2$ and $u_3$:

Step 5: Set up the eigenvalue problem

For free undamped vibration, the dynamic equations of motion simplify to:

We assume a harmonic solution of the form:

where $\boldsymbol{\phi}$ is the eigenvector (mode shape) and $\omega$ is the natural frequency.

Differentiating twice with respect to time:

Substituting these expressions back into the equation of motion yields:

For a non-trivial solution (where $\boldsymbol{\phi} \neq \mathbf{0}$), the determinant of the characteristic matrix must be zero:

Let us define a dimensionless eigenvalue parameter $\lambda$:

This definition allows us to rewrite $\omega^2 \mathbf{M}$ as:

Now, substitute this into the characteristic equation:

Dividing out the constant scalar factor $\frac{2AE}{L}$:

Step 6: Solve the characteristic equation for eigenvalues

Expanding the determinant:

We solve this quadratic equation using the quadratic formula:

Thus, the two roots are:

Step 7: Calculate natural frequencies $\omega_1, \omega_2$

Using our definition $\omega^2 = \frac{24 E \lambda}{\rho L^2} \implies \omega = \frac{\sqrt{24\lambda}}{L} \sqrt{\frac{E}{\rho}}$:

First Natural Frequency $\omega_1$:

Second Natural Frequency $\omega_2$:

Step 8: Calculate eigenvectors (mode shapes)

We substitute each eigenvalue back into the equation $(\mathbf{K} - \omega^2 \mathbf{M}) \boldsymbol{\phi} = \mathbf{0}$ to find the corresponding eigenvectors:

1. First Mode Shape $\boldsymbol{\phi}_1$ (using $\lambda_1 = 0.07418$):

From the first row equation:

If we set the deflection at Node 3 (free end) as $\phi_{31} = 1.0$, then:

Thus, the first mode shape is:

This mode shape shows that both Node 2 and Node 3 displace in the same direction, which represents the first fundamental breathing mode of the bar.

2. Second Mode Shape $\boldsymbol{\phi}_2$ (using $\lambda_2 = 1.92582$):

From the first row equation:

If we set $\phi_{32} = 1.0$, then:

Thus, the second mode shape is:

In this mode, Node 2 and Node 3 displace in opposite directions. The sign change indicates the presence of a node point (a point along the length of the bar with zero displacement) within the second element, representing a higher-frequency vibration mode.

Step 9: Compare with the exact analytical solution

The exact analytical frequencies of a continuous cantilever bar are given by:

Our 2-element consistent mass calculation: $\omega_1 \approx \frac{1.334}{L} \sqrt{\frac{E}{\rho}}$ (Error = $-15.1\%$)

Our 2-element consistent mass calculation: $\omega_2 \approx \frac{6.798}{L} \sqrt{\frac{E}{\rho}}$ (Error = $+44.3\%$)

*Note*: Why is the error relatively high? Because we only used 2 elements to model a continuous bar! To capture higher modes (like $n=2$), we need more elements to model the high curvature in the displacement field. If we increase the number of elements, our FEA natural frequencies will rapidly converge to the exact analytical solutions.