Frequency analysis#
The function freq() solves the eigenvalue problem
\((\boldsymbol{K} + \lambda^2 \boldsymbol{M}) \boldsymbol{u} = \boldsymbol{0}\), where \(\lambda^2 = -\omega_n^2\)
and \(\omega_n\) is the natural frequency in rad/s. The returned lambda2
are therefore negative for a stable structure, and
omega = np.sqrt(-lambda2) # rad/s
f = omega/(2*np.pi) # Hz
with np.sqrt(-lambda2 + 0j) when lambda2 may be complex or positive,
e.g. beyond a flutter or a divergence load, see Follower loads and unsymmetric matrices.
The code below is extracted from one of the structsolve unit tests:
def test_freq_simple():
"""Test frequency analysis with known natural frequencies
For a diagonal system: [K]{u} + omega^2 [M]{u} = 0
K = diag(k1, k2, ...), M = diag(m1, m2, ...)
omega_i = sqrt(ki/mi)
"""
n = 10
k_vals = np.arange(1, n + 1, dtype=float) * 100
m_vals = np.ones(n)
K = csc_matrix(np.diag(k_vals))
M = csc_matrix(np.diag(m_vals))
eigvals, eigvecs = freq(K, M, silent=True, num_eigvalues=5,
sparse_solver=True)
# freq returns lambda2 where [K] + lambda2 [M] = 0
# so lambda2 = -omega^2 = -k/m
# The sorted natural frequencies should match sqrt(k_vals)
assert eigvecs.shape[0] == n
Options#
symmetric:False, the default, uses the general solvers, valid for unsymmetric matrices, and returns complex arrays.Trueuses the solvers of symmetric matrices, which require a positive definiteM, and returns real arrays;Nonechecks the symmetry of the matrices.Normalization: with
symmetric=True, the eigenvectors are normalized with respect toM, \(\boldsymbol{u}_i^\top \boldsymbol{M} \boldsymbol{u}_j = \delta_{ij}\). Otherwise the sparse solver returns eigenvectors normalized with respect toMand the dense solver eigenvectors with unit Euclidean norm, seefreq(); an application that needs a given normalization applies it to the returned eigenvectors.sort:sort=True, the default, sortslambda2and the eigenvectors together by increasing natural frequency and removes the eigenpairs with \(\omega_n \le 10^{-6}\) rad/s, e.g. rigid-body modes, such that fewer thannum_eigvaluesmay be returned.sort=Falsekeeps the order of the solver.num_eigvalues: the number of eigenvalues computed by the sparse solver, those closest to a negative shift, i.e. the lowest frequencies. See Accuracy of the eigenvalue solutions for the verification of the results.
The sorted eigenpairs, with and without a rigid-body mode:
@pytest.mark.parametrize('sparse_solver', [True, False])
@pytest.mark.parametrize('symmetric', [True, False])
def test_freq_sort(sparse_solver, symmetric):
"""sort=True sorts and filters the eigenvalues together with the
eigenvectors: increasing natural frequencies, without the rigid-body
mode"""
rng = np.random.RandomState(3)
n = 12
k = rng.permutation(np.arange(1, n + 1, dtype=float))*100
K = np.diag(k)
# a free-free pair of springs: one rigid-body mode
K[:2, :2] = 50.*np.array([[1., -1.], [-1., 1.]])
M = np.eye(n)
K, M = csc_matrix(K), csc_matrix(M)
lambda2, eigvecs = freq(K, M, silent=True, num_eigvalues=6, sort=True,
sparse_solver=sparse_solver, symmetric=symmetric)
omega2 = (-lambda2).real
assert lambda2.shape[0] == eigvecs.shape[1]
assert np.all(omega2 > 1e-6)
# exactly non-decreasing, which implies the non-decreasing rounded
# omega_n of the primary sort key: the double eigenvalue omega2 = 100,
# of the springs and of k = 100, which the dense eig() returns with a
# round-off difference in an order depending on the LAPACK build, is
# ordered by the exact tie-break of the sort key
assert np.all(np.diff(omega2) >= 0)
lambda2_u, _ = freq(K, M, silent=True, num_eigvalues=6, sort=False,
sparse_solver=sparse_solver, symmetric=symmetric)
w2_u = np.sort((-lambda2_u).real)
np.testing.assert_allclose(omega2, w2_u[w2_u > 1e-6][:omega2.size],
rtol=1e-8)
# each eigenvalue with its eigenvector
res = np.linalg.norm(K @ eigvecs + (M @ eigvecs)*lambda2, axis=0)
assert np.all(res <= 1e-6*np.linalg.norm(K @ eigvecs, axis=0))
Prestressed structures#
The frequencies of a prestressed structure come from the stiffness matrix
of the prestressed state, e.g. K + lam*KG with the geometric stiffness
matrix KG of a reference load, see Linear buckling with a prestress, or the tangent
stiffness matrix of a non-linear analysis. A compressive prestress reduces
the frequencies, and the lowest one vanishes at the buckling load.
For the plate of panels saved in
tests/data, with a lumped mass matrix, the frequencies decrease under a
compressive prestress and the lowest one vanishes at the critical load
multiplier of lb():
def test_freq_prestressed_plate():
"""Natural frequencies of the plate of panels saved in tests/data, with
a lumped mass matrix, under a compressive prestress lam*KG: they
decrease with lam, and the lowest one vanishes linearly at the critical
load multiplier of lb"""
import os
from scipy.sparse import diags, load_npz
from structsolve import lb
data = os.path.join(os.path.dirname(__file__), 'data')
K = load_npz(os.path.join(data, 'plate_ssss_Nxx_K.npz')).tocsr()
KG = load_npz(os.path.join(data, 'plate_ssss_Nxx_KG.npz')).tocsr()
rng = np.random.RandomState(3)
M = diags(1. + rng.rand(K.shape[0])).tocsr()
# uniaxial compression: KG is negative semi-definite
assert np.linalg.eigvalsh(KG.toarray()).max() <= 1e-10*abs(KG).max()
lambda_1 = lb(K, KG, silent=True, sparse_solver=False)[0][0]
omega2 = []
for frac in [0., 0.5, 0.9, 1 - 1e-6, 1 - 1e-7]:
lambda2, eigvecs = freq(K + frac*lambda_1*KG, M, silent=True,
num_eigvalues=5, symmetric=True)
omega = np.sqrt(-lambda2) # rad/s
f_Hz = omega/(2*np.pi) # Hz
assert np.all(np.diff(f_Hz) >= 0)
# symmetric=True: eigenvectors normalized with respect to M
np.testing.assert_allclose(eigvecs.T @ M @ eigvecs, np.eye(5),
atol=1e-8)
omega2.append(omega**2)
omega2 = np.array(omega2)
# the frequencies decrease with the compressive prestress
assert np.all(np.diff(omega2, axis=0) <= 1e-12*omega2[0, 0])
# the lowest one vanishes linearly at the critical load multiplier
np.testing.assert_allclose(omega2[-2, 0]/omega2[-1, 0], 10., rtol=1e-2)