Dynamic analysis of bar and frame structures constitutes one of the fundamental areas of structural mechanics and the finite element method (FEM). In contrast to static problems, where loads are constant in time, dynamic analysis considers the influence of time-dependent forces, inertia effects, and vibration damping. Such problems occur in civil engineering, mechanical engineering, aerospace engineering, and structural dynamics. The fundamental equation of motion of a discretized structural system obtained using the finite element method is
$$\mathbf{B}\ddot{\mathbf{q}} + \mathbf{C}\dot{\mathbf{q}} + \mathbf{K}\mathbf{q} = \mathbf{F}(t)$$where:
• $\mathbf{B}$ -- mass matrix,
• $\mathbf{C}$ -- damping matrix,
• $\mathbf{K}$ -- stiffness matrix,
• $\mathbf{q}$ -- vector of nodal displacements,
• $\dot{\mathbf{q}}$ -- vector of nodal velocities,
• $\ddot{\mathbf{q}}$ -- vector of nodal accelerations,
• $\mathbf{F}(t)$ -- vector of external time-dependent loads.
This equation describes the motion of a multi-degree-of-freedom system and forms the basis for numerical vibration analysis of bar and frame structures.
For a single element, the displacement field can be written as
$$u(x,t)=\mathbf{N}(x)\mathbf{q}(t)$$where:
• $\mathbf{N}(x)$ -- matrix of shape functions,
• $\mathbf{q}_e(t)$ -- vector of element nodal displacements.
The stiffness matrix relates nodal forces and bending moments to nodal transverse displacements and rotations.
$$\mathbf{K}= \frac{EI}{L^3} \begin{bmatrix} 12 & 6L & -12 & 6L \\ 6L & 4L^2 & -6L & 2L^2 \\ -12 & -6L & 12 & -6L \\ 6L & 2L^2 & -6L & 4L^2 \end{bmatrix}$$In dynamic analysis, proper modeling of inertia effects is essential. Two principal forms of the mass matrix are commonly used in FEM:
• lumped mass matrix,
• consistent mass matrix.
The choice of the mass matrix formulation significantly influences the accuracy and computational efficiency of dynamic simulations.
The consistent mass matrix is derived using the same interpolation functions that are employed in the displacement approximation. Therefore, the inertia distribution is fully consistent with the finite element interpolation.
The derivation begins from the kinetic energy of the element
$$T=\frac{1}{2} \int_V \rho \dot{u}^T \dot{u} \, dV$$where:
• $\rho$ -- material density,
• $V$ -- element volume,
• $\dot{u}$ -- velocity field.
Using the finite element approximation
\begin{equation}
u(x,t)=\mathbf{N}(x)\mathbf{q}(t)
\end{equation}
the velocity field becomes
\begin{equation}
\dot{u}(x,t)=\mathbf{N}(x)\dot{\mathbf{q}}(t)
\end{equation}
Substituting this relation into the kinetic energy expression yields
$$T= \frac{1}{2} \int_0^L \rho A \left( \mathbf{N}\dot{\mathbf{q}} \right)^T \left( \mathbf{N}\dot{\mathbf{q}} \right) dx$$which becomes
$$T= \frac{1}{2} \dot{\mathbf{q}}^T \left( \int_0^L \rho A \mathbf{N}^T \mathbf{N} \, dx \right) \dot{\mathbf{q}}$$Comparing this expression with the standard quadratic form of kinetic energy
\begin{equation}
T=
\frac{1}{2}
\dot{\mathbf{q}}^T
\mathbf{B}_e
\dot{\mathbf{q}}
\end{equation}
the consistent mass matrix of the beam element is obtained as
\begin{equation}
\mathbf{B}=
\int_0^L
\rho A
\mathbf{N}^T
\mathbf{N}
\, dx
\end{equation}
After carrying out the integration, the consistent mass matrix for the Euler--Bernoulli beam element becomes
\begin{equation}
\mathbf{B}=
\frac{\rho AL}{420}
\begin{bmatrix}
156 & 22L & 54 & -13L \\
22L & 4L^2 & 13L & -3L^2 \\
54 & 13L & 156 & -22L \\
-13L & -3L^2 & -22L & 4L^2
\end{bmatrix}
\end{equation}
where:
• $\rho$ -- material density,
• $A$ -- cross-sectional area,
• $L$ -- beam element length.
The consistent mass matrix accurately represents distributed inertia effects and preserves dynamic coupling between translational and rotational degrees of freedom.
For a one-dimensional bar element with linear shape functions
\begin{equation}
N_1(x)=1-\frac{x}{L},
\qquad
N_2(x)=\frac{x}{L}
\end{equation}
the shape function matrix is
$$\mathbf{N}(x)= \begin{bmatrix} N_1(x) & N_2(x) \end{bmatrix}$$which leads to the consistent mass matrix of the bar element:
$$\mathbf{B}= \frac{\rho AL}{6} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$$The lumped mass matrix is obtained by concentrating the total mass of the element at its nodes. In this approach, inertia coupling between nodes is neglected, which leads to a diagonal mass matrix.
The total mass of the beam element is
$$m=\rho AL$$A simplified lumped mass matrix for the beam element can be written as
$$\mathbf{B}= \frac{\rho AL}{2} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}$$The reason for neglecting inertia coupling in the lumped mass formulation can be understood by starting from the consistent mass expression
$$\mathbf{B}=\int_0^L \rho A \mathbf{N}^T \mathbf{N}\, dx$$In this formulation, the off-diagonal terms arise from products of different shape functions, i.e.
\begin{equation}
B_{ij} \sim \int_0^L N_i(x) N_j(x)\, dx \quad (i \neq j)
\end{equation}
These terms represent inertia coupling between nodal degrees of freedom, meaning that acceleration at one node induces inertia forces at the other node.
In the lumped mass approach, the mass distribution is replaced by equivalent nodal point masses located at the element nodes. Mathematically, this corresponds to the approximation
$$\int_0^L \rho A N_i(x)N_j(x)\,dx \;\approx\; 0 \quad \text{for } i \neq j$$and
$$\int_0^L \rho A N_i(x)^2\,dx \;\rightarrow\; \frac{m}{2}$$This approximation is equivalent to assuming that the inertia force field is concentrated at discrete points rather than distributed along the element length. As a consequence, the kinetic energy becomes
$$T \approx \frac{1}{2}\sum_{k=1}^{n} m_k \dot{q}_k^2$$which contains no cross-terms of the form $\dot{q}_i \dot{q}_j$. Therefore, no dynamic coupling between nodes is present in the mass representation, and the resulting mass matrix becomes diagonal.
For a bar element:
$$\mathbf{B}= \frac{\rho AL}{2} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$$The lumped formulation can also be obtained from the consistent matrix using mass diagonalization techniques such as:
• row summation,
• diagonal scaling,
• nodal quadrature methods.
The principal advantage of the lumped mass matrix is its computational efficiency. Since the matrix is diagonal:
• matrix inversion becomes trivial,
• explicit time integration schemes become highly efficient,
• storage requirements are reduced.
For this reason, lumped mass matrices are widely used in large-scale transient dynamic simulations and explicit finite element methods.
However, the simplification reduces accuracy in representing distributed inertia and may lead to less accurate predictions of higher natural frequencies and mode shapes. In practical finite element analysis, the choice between the two formulations depends on the type of dynamic problem, desired accuracy, and available computational resources.
In practical engineering applications, the mass distribution along a beam or bar is often non-uniform due to variations in cross-section, material properties, or attached concentrated masses. Let the distributed mass per unit length be described by
$$\mu(x)=\rho(x)A(x)$$so that the total mass of the element is
$$m=\int_0^L \mu(x)\,dx$$In the lumped mass approach, this distributed mass is replaced by equivalent nodal masses located at the element nodes. For a two-node beam or bar element, the nodal masses $m_1$ and $m_2$ are determined such that the total mass and the first moment of mass are preserved approximately.
A common approach is based on consistent weighting of the shape functions:
\begin{equation}
m_1 = \int_0^L N_1(x)\,\mu(x)\,dx,
\qquad
m_2 = \int_0^L N_2(x)\,\mu(x)\,dx
\end{equation}
where $N_1(x)$ and $N_2(x)$ are the linear shape functions of the element.
These definitions ensure that the equivalent nodal masses reproduce the same inertia work as the distributed system under rigid body motion.
In practice, when $\mu(x)$ is not known analytically, it is often approximated numerically using Gauss integration or midpoint evaluation. A simple two-point quadrature gives
$$m_1 \approx \frac{L}{2}\,\mu(x_1), \qquad m_2 \approx \frac{L}{2}\,\mu(x_2)$$where $x_1$ and $x_2$ are representative sampling points within the element (often the nodes or Gauss points).
A particularly common engineering approximation assumes that each node carries half of the total element mass:
$$m_1 = m_2 = \frac{1}{2}\int_0^L \mu(x)\,dx$$This approach preserves the total mass exactly but neglects the first moment of inertia distribution, which may introduce errors in dynamic response when the mass is strongly non-uniform.
Damping describes the dissipation of vibration energy. In engineering practice, Rayleigh damping is most commonly applied:
$$\mathbf{C}=\alpha \mathbf{B}+\beta \mathbf{K}$$where:
• $\alpha$ -- mass proportional damping coefficient,
• $\beta$ -- stiffness proportional damping coefficient.
This model provides a convenient approximation of the damping properties of real structures.
In the absence of external excitation and damping, the equation of motion reduces to
\begin{equation}
\mathbf{B}\ddot{\mathbf{q}}+\mathbf{K}\mathbf{q}=0
\end{equation}
Assuming a harmonic solution of the form
\begin{equation}
\mathbf{q}(t)=\boldsymbol{\phi}e^{i\omega t}
\end{equation}
leads to the eigenvalue problem
$$(\mathbf{K}-\omega^2\mathbf{B})\boldsymbol{\phi}=0$$where:
• $\omega$ -- natural frequencies,
• $\boldsymbol{\phi}$ -- vibration mode shapes.
Solving this problem allows determination of the dynamic characteristics of the structure.
When the structure is subjected to time-dependent external loads, the excitation term must be included:
\begin{equation}
\mathbf{F}(t)\neq 0
\end{equation}
Dynamic excitation may include:
• harmonic excitation,
• impulsive loads,
• seismic excitation,
• random loads.
For harmonic excitation
$$\mathbf{F}(t)=\mathbf{F}_0\sin(\Omega t)$$the resonance phenomenon may occur when the excitation frequency approaches one of the natural frequencies of the structure.
The numerical solution of the dynamic equation of motion requires time integration procedures. In FEM practice, both explicit and implicit methods are used.
One of the most widely applied implicit schemes is the Newmark method. It assumes the following relations:
$$\mathbf{q}_{n+1}= \mathbf{q}_n+ \Delta t \dot{\mathbf{q}}_n+ \Delta t^2 \left[ \left( \frac{1}{2}-\beta \right) \ddot{\mathbf{q}}_n+ \beta \ddot{\mathbf{q}}_{n+1} \right]$$and
$$\dot{\mathbf{q}}_{n+1}= \dot{\mathbf{q}}_n+ \Delta t \left[ (1-\gamma)\ddot{\mathbf{q}}_n+ \gamma \ddot{\mathbf{q}}_{n+1} \right]$$where:
• $\Delta t$ -- time step,
• $\beta$, $\gamma$ -- algorithm parameters.
The Newmark method is characterized by good numerical stability and is widely used in structural dynamics.
Modal analysis consists in decomposing the dynamic response into a superposition of vibration mode shapes. The displacement vector can be expressed as
$$\mathbf{q}(t)=\mathbf{\Phi}\boldsymbol{\eta}(t)$$where:
• $\mathbf{\Phi}$ -- matrix of mode shapes,
• $\boldsymbol{\eta}(t)$ -- modal coordinates.
After transformation, the equation of motion becomes decoupled into independent modal equations, significantly reducing computational effort.
The modal approach is especially effective for large structural systems.
Dynamic analysis of bar and frame structures is widely used in engineering applications, including:
• seismic analysis of buildings,
• design of bridges and towers,
• vibration analysis of machines,
• aerospace structural dynamics,
• impact and impulsive load analysis,
• investigation of dynamic stability.