Modal analysis and periodic boundaries

This chapter covers two workflows that sit on top of the standard assembly and solver path:

• extracting modal data from the generalized eigenvalue problem,

• enforcing periodic constraints for representative volume element (RVE) calculations and homogenization.

The code paths are separate, but they share the same global matrix conventions and one-based legacy tables. solve_modal() works with reduced free-DOF blocks; the periodic module builds a constraint system and solves the resulting saddle-point problem with Lagrange multipliers.

1. Modal Analysis

1.1 Generalized eigenvalue problem

The free vibration equation without damping is:

\[ \mathbf{M}\ddot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{0} \]

Assuming harmonic motion \(\mathbf{u}(t) = \boldsymbol{\phi} e^{i\omega t}\):

\[ -\omega^2 \mathbf{M}\boldsymbol{\phi} e^{i\omega t} + \mathbf{K}\boldsymbol{\phi} e^{i\omega t} = \mathbf{0} \]

This yields the generalized eigenvalue problem:

\[ \mathbf{K}\boldsymbol{\phi} = \omega^2 \mathbf{M}\boldsymbol{\phi} \]

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

1.2 Orthogonality of modes

For distinct eigenvalues \(\omega_m^2 \neq \omega_n^2\), the mode shapes satisfy:

\[ \boldsymbol{\phi}_m^T \mathbf{M} \boldsymbol{\phi}_n = 0, \quad \boldsymbol{\phi}_m^T \mathbf{K} \boldsymbol{\phi}_n = 0 \quad (m \neq n) \]

Proof: From the eigenvalue equations: $\( \mathbf{K}\boldsymbol{\phi}_m = \omega_m^2 \mathbf{M}\boldsymbol{\phi}_m, \quad \mathbf{K}\boldsymbol{\phi}_n = \omega_n^2 \mathbf{M}\boldsymbol{\phi}_n \)$

Pre-multiply the first by \(\boldsymbol{\phi}_n^T\) and the second by \(\boldsymbol{\phi}_m^T\): $\( \boldsymbol{\phi}_n^T \mathbf{K}\boldsymbol{\phi}_m = \omega_m^2 \boldsymbol{\phi}_n^T \mathbf{M}\boldsymbol{\phi}_m \)\( \)\( \boldsymbol{\phi}_m^T \mathbf{K}\boldsymbol{\phi}_n = \omega_n^2 \boldsymbol{\phi}_m^T \mathbf{M}\boldsymbol{\phi}_n \)$

Since \(\mathbf{K}\) and \(\mathbf{M}\) are symmetric, transposing gives equal LHS. Subtracting: $\( (\omega_m^2 - \omega_n^2) \boldsymbol{\phi}_n^T \mathbf{M}\boldsymbol{\phi}_m = 0 \)$

For \(\omega_m \neq \omega_n\): \(\boldsymbol{\phi}_n^T \mathbf{M}\boldsymbol{\phi}_m = 0\) ∎

1.3 Mass normalization

Modes are typically normalized such that:

\[ \boldsymbol{\phi}_n^T \mathbf{M} \boldsymbol{\phi}_n = 1 \]

This implies: $\( \boldsymbol{\phi}_n^T \mathbf{K} \boldsymbol{\phi}_n = \omega_n^2 \)$

1.4 Modal participation factors

The participation factor for mode \(n\) in direction \(j\) quantifies how much the mode contributes to the response under uniform excitation:

\[ \Gamma_{n,j} = \frac{\boldsymbol{\phi}_n^T \mathbf{M} \mathbf{r}_j}{\boldsymbol{\phi}_n^T \mathbf{M}\boldsymbol{\phi}_n} \]

where \(\mathbf{r}_j\) is a unit influence vector (ones for DOFs in direction \(j\), zeros elsewhere).

For mass-normalized modes (\(\boldsymbol{\phi}_n^T \mathbf{M}\boldsymbol{\phi}_n = 1\)):

\[ \Gamma_{n,j} = \boldsymbol{\phi}_n^T \mathbf{M} \mathbf{r}_j \]

1.5 Effective modal mass

The effective modal mass in direction \(j\) is:

\[ M_{eff,n,j} = \Gamma_{n,j}^2 \cdot (\boldsymbol{\phi}_n^T \mathbf{M}\boldsymbol{\phi}_n) = \Gamma_{n,j}^2 \]

The sum of effective masses equals the total mass:

\[ \sum_{n=1}^{N} M_{eff,n,j} = \mathbf{r}_j^T \mathbf{M} \mathbf{r}_j = M_{total,j} \]

1.6 Implementation details

solve_modal(K, M, n_modes=10, C_bc=None, dof=2, sigma=0.0) returns a ModalResult with:

• eigenvalues: \(\omega^2\)

• omega: circular frequencies (rad/s)

• freq_hz: \(f = \omega/(2\pi)\) (Hz)

• period: \(T = 1/f\) (seconds)

• mode_shapes: mass-normalized eigenvectors

• participation: \(\Gamma_{n,j}\) per direction

• effective_mass: \(M_{eff,n,j}\) per direction

The solver reduces the system by removing constrained DOFs before solving, then expands results back to full size. it.

2. Periodic boundaries

The periodic module is organized around three steps:

1. identify matching nodes on opposite faces,

2. build the linear constraint system,

3. solve the augmented problem and average the result.

Pairing nodes

find_periodic_pairs(X, axis, tol=1e-6) compares the minimum and maximum faces of a domain along a selected axis and returns 1-based node pairs. The implementation uses cKDTree when available and the node count is large enough; otherwise it falls back to a direct distance search. The tolerance is scaled by the domain size, which makes the matching rule stable across different units.

find_all_periodic_pairs(X, periodic_axes, tol=1e-6) repeats that process for multiple axes and returns a dictionary keyed by axis.

check_periodic_mesh(X, axis, tol=1e-6) returns a plain report dictionary with valid, n_left, n_right, max_mismatch, and message. This is the first function to call if a mesh should be periodic but pairing fails.

Fixing the rigid body mode

Pure periodic constraints do not remove rigid translations by themselves. fix_corner(X, C_existing, dof) adds zero-displacement rows at the corner node closest to the minimum coordinate corner so a separate direct solve does not remain singular. solve_periodic() does not take a boundary-condition table; it handles the periodic augmentation through solve_lag_general().

Constraint matrices and macro strain

periodic_constraints(X, pairs, dof, eps_macro=None) converts the node pairs into the algebraic system

\[ \mathbf{G}\mathbf{u} = \mathbf{Q} \]

where each row enforces one DOF pair. The right node gets +1, the left node gets -1, and the macro strain contribution is written into Q. The Voigt strain vector is converted to a tensor by _voigt_to_tensor(), so the shear components are split consistently between the off-diagonal tensor entries.

apply_macro_strain() is just a convenience wrapper that returns the Q vector from the same logic.

solve_periodic(K, p, X, pairs, dof, eps_macro=None, return_lagrange=False) calls periodic_constraints() and then passes K, p, G, and Q to solve_lag_general(). That means the periodic solve is not a special-case solver. It is a standard saddle-point solve with a different constraint matrix.

3. Volume averages and homogenization

volume_average_stress() and volume_average_strain() are the reduction layer that converts element-level fields into domain averages. The current implementation supports t3 and q4 elements.

• For t3, the code reuses the triangular batch geometry helper and weights element values by area.

• For q4, the code averages the four Gauss-point values and estimates the area with a shoelace formula on the nodal polygon.

homogenize(K, T, X, G_mat, pairs, dof, element_type="q4") then drives the RVE through unit macro-strain cases and assembles the effective constitutive matrix column by column. In 2D the result is 3 x 3; in 3D it is 6 x 6.

The important part of the implementation is the flow:

1. build a unit macro strain vector,

2. call solve_periodic() with zero external load,

3. average the resulting stress with volume_average_stress(),

4. insert the result into one column of C_eff.

That structure is why the function is easy to extend and easy to audit.

4. Practical workflow

The normal sequence for an RVE solve is:

import numpy as np
from femlabpy.periodic import (
    check_periodic_mesh,
    find_periodic_pairs,
    solve_periodic,
    homogenize,
)

report = check_periodic_mesh(X, axis=0)
if not report["valid"]:
    raise ValueError(report["message"])

pairs_x = find_periodic_pairs(X, axis=0)
pairs_y = find_periodic_pairs(X, axis=1)
pairs = np.vstack([pairs_x, pairs_y])
u = solve_periodic(K, p_zero, X, pairs, dof=2, eps_macro=np.array([1.0, 0.0, 0.0]))
C_eff = homogenize(K, T, X, G_mat, pairs, dof=2, element_type="q4")

The same pattern is used in the repository examples: verify the mesh, pair the faces, solve the constrained system, and average the stress response into an effective stiffness tensor. If you need a separate reference constraint for a non-periodic solve, add fix_corner() to your own C table before calling the linear solver.