A general-purpose element-based approach to compute dispersion relations in periodic materials with existing finite element codes

In most of standard Finite Element (FE) codes it is not easy to calculate dispersion relations from periodic materials. Here we propose a new strategy to calculate such dispersion relations with available FE codes using user element subroutines. Typically, the Bloch boundary conditions are applied to the global assembled matrices of the structure through a transformation matrix or row-and-column operations. Such a process is difficult to implement in standard FE codes since the user does not have access to the global matrices. In this work, we apply those Bloch boundary conditions directly at the elemental level. The proposed strategy can be easily implemented in any FE code. This strategy can be used either in real or complex algebra solvers. It is general enough to permit any spatial dimension and physical phenomena involving periodic structures. A detailed process of calculation and assembly of the elemental matrices is shown. We verify our method with available analytical solutions and external numerical results, using different material models and unit cell geometries

💡 Research Summary

The paper addresses a long‑standing practical difficulty in commercial finite‑element (FE) software: the inability to directly manipulate global stiffness and mass matrices in order to impose Bloch periodic boundary conditions (BBCs) for the analysis of periodic media. Traditional Bloch analyses require a transformation matrix or row‑and‑column operations on the assembled global matrices to enforce the phase shift e^{i k·a} between opposite sides of the unit cell. Because most commercial packages hide the global matrices from the user, such operations are either impossible or require cumbersome work‑arounds (e.g., multi‑point constraints) that limit the range of physical problems that can be studied.

The authors propose a fundamentally different strategy: impose the Bloch conditions at the element level, thereby eliminating any need to access or modify the global system. The key idea is to split the element connectivity information into two distinct operators. The first, C_c (coordinate connectivity), is used exactly as in a standard FE formulation to compute element stiffness and mass matrices from the nodal coordinates. The second, C_a (assembly connectivity), is employed during the assembly step and is defined so that degrees of freedom belonging to equivalent nodes in adjacent periodic cells are automatically merged with the appropriate complex phase factor. In practice, for each element that contains nodes on a periodic boundary, the authors construct C_c in the usual way, but replace C_a with a list that maps each “slave” node to its “master” counterpart multiplied by the phase shift PS_{i→j}=e^{i k·a}. Consequently, when the FE code assembles the element contributions, the global matrices K_R and M_R that emerge are already reduced and contain the Bloch phase information; no post‑assembly row‑column reduction is required.

This element‑based approach is compatible with both real‑algebra solvers (using the 2‑mesh technique of Aberg & Gudmundson, which splits the complex eigenvalue problem into two real problems) and complex‑algebra solvers that directly handle complex matrices. Moreover, because the method works at the element level, it is agnostic to the underlying physics: scalar wave equations, vector elastodynamics, and even higher‑order theories such as Cosserat or micropolar elasticity can be treated without modification. The authors demonstrate this versatility through a series of verification studies:

1. Classical elastodynamics – 2‑D in‑plane wave propagation in a rectangular unit cell, comparing the computed band diagram with analytical solutions and published numerical results.
2. Cosserat/micropolar media – Both 2‑D and 3‑D unit cells are examined, showing that the additional rotational degrees of freedom and couple‑stress effects are correctly captured in the dispersion curves.
3. Geometric diversity – Unit cells of different shapes (square, triangular, and three‑dimensional brick) and material contrasts are analyzed, confirming that the method works irrespective of mesh topology.
4. Cell‑size sensitivity – By varying the size of the periodic cell while keeping the physical wavelength constant, the authors verify that the approach yields consistent dispersion results, indicating numerical robustness.

All test cases exhibit excellent agreement with reference data from the literature, confirming that the element‑level implementation does not sacrifice accuracy. To facilitate adoption, the authors provide a fully functional version of the FEAP‑pv code that includes a user‑element subroutine implementing the proposed scheme, together with detailed usage instructions and an example of a homogeneous material unit cell.

In summary, the paper delivers a practical, code‑agnostic methodology for Bloch analysis that sidesteps the need for global matrix access. By embedding the Bloch phase directly into the element assembly process, the technique can be deployed in any commercial FE package that supports user‑defined elements. This opens the door for researchers and engineers to perform band‑structure calculations for a wide range of periodic structures—acoustic, elastic, electromagnetic, or multi‑physics—without requiring specialized proprietary software. Future work may extend the framework to nonlinear materials, coupled multiphysics problems, and non‑rectilinear lattices, further broadening its applicability in the emerging fields of metamaterials and phononic crystals.