Harmonic (frequency-response) analysis#

harmonic_solve() and the underlying femorph_solver.solvers.harmonic.solve_harmonic_modal() evaluate the steady-state frequency-response function

\[\bigl(K \;-\; \omega^{2}\,M \;+\; i\,\omega\,C\bigr)\, u(\omega) \;=\; F\]

at every excitation frequency in a user-supplied list. The solve is modal-superposition: the load and the response are projected onto the lowest n_modes mass-orthonormal eigenmodes of the undamped system, and a per-mode complex receptance is back-substituted.

Algorithm#

The undamped eigenproblem \(K \varphi_{n} = \omega_{n}^{2} M \varphi_{n}\) is solved once (via Modal analysis) and the eigenvectors are normalised so \(\Phi^{T} M \Phi = I\). The modal response coordinate at each excitation frequency \(\omega\) is

\[\eta_{n}(\omega) \;=\; \frac{\varphi_{n}^{T}\,F} {\omega_{n}^{2} - \omega^{2} + 2\,j\,\zeta_{n}\,\omega_{n}\,\omega}\]

and the physical response is the modal sum \(u(\omega) = \Phi\,\eta(\omega)\). Every excitation frequency re-uses the same eigenpairs — the cost of a 200-point sweep is dominated by the one-time modal solve, not the per-frequency back-substitution.

API#

The shortest path lives on Model. Build the model, fix the BCs, stage the load, then call harmonic_solve():

import numpy as np
import femorph_solver as fs
from femorph_solver import ELEMENTS

m = fs.Model.from_grid(grid)
m.assign(ELEMENTS.HEX8(integration="enhanced_strain"),
         {"EX": 2.0e11, "PRXY": 0.30, "DENS": 7850.0})
m.fix(where=clamped_face, dof=["UX", "UY", "UZ"])
m.apply_force(tip_node, fz=-1.0)

res = m.harmonic_solve(
    frequencies=np.linspace(0.0, 1000.0, 200),  # Hz
    n_modes=20,
    modal_damping_ratio=0.02,
)

res.frequency        # (n_freq,) Hz
res.displacement     # (n_freq, N) complex
res.frf_at(dof_idx)  # (n_freq,) complex receptance at one DOF
res.modal            # the underlying ModalResult that spanned the basis

For low-level work directly on the K / M matrices, the same solver is exposed as femorph_solver.solvers.harmonic.solve_harmonic_modal() and takes the prescribed-DOF dictionary explicitly.

Damping#

Two damping models are accepted; both can be combined and their contributions add.

Rayleigh damping#

Pass an (α, β) pair so the per-mode ratio becomes \(\zeta_{n} = \alpha/(2\,\omega_{n}) + \beta\,\omega_{n}/2\):

res = m.harmonic_solve(frequencies=freq, n_modes=20,
                       rayleigh=(0.05, 1e-4))

Pure rigid-body modes (\(\omega_{n} = 0\)) drop the \(\alpha/(2\omega_{n})\) term — singular by construction — and keep whatever modal_damping_ratio was supplied at that index.

Convergence and basis size#

The truncation error at frequency \(\omega\) is governed by every mode whose \(\omega_{n}\) exceeds \(\omega\). In practice:

Retain at least the first mode whose \(f_{n}\) exceeds the largest excitation frequency in your sweep.
Add another factor-of-two cushion when the response near a higher pole matters (truncation manifests as a soft underestimate of the off-resonance amplitude, not a hard miss).
Static-correction terms are not yet implemented; below the first retained pole the truncation residual is bounded but not zero.

Persistence#

HarmonicResult.save() writes the complex per-frequency response to a single .pv file alongside the underlying frequency grid and applied damping ratios. The on-disk format splits each per-frequency complex array into separate displacement_re_NNNN and displacement_im_NNNN point-data entries since pyvista’s field arrays don’t carry complex dtype natively; reading the file back via femorph_solver.result.read() returns a HarmonicResultDisk with a displacement() accessor that recombines the real and imaginary halves.