Linear buckling and natural frequencies#

The eigenvalue analyses use the matrices of the free dofs, e.g. KC0uu = KC0[bu, :][:, bu], and are solved with linear_buckling(), natural_frequency() and complex_frequency() of pyfe3d.solver.eigen, which return the eigenvectors with the size of the matrices received, to be expanded by the user, u[bu] = eigvecs[:, i].

Linear buckling#

A linear buckling analysis has three steps:

  1. the linear static solution \(\boldsymbol{u}_0\) of a reference load \(\boldsymbol{F}_{ref}\), \(\boldsymbol{K}_{{C0}_{uu}} \boldsymbol{u}_{0_u} = \boldsymbol{F}_{{ref}_u}\), see Overview and workflow;

  2. the geometric stiffness matrix of the prestress of \(\boldsymbol{u}_0\): for each element, update_probe_xe(x), update_probe_ue(u0) and update_KG(KGr, KGc, KGv, prop). For the shells, a known constant membrane state can be given instead with update_KG_given_stress(Nxx, Nyy, Nxy, KGr, KGc, KGv), with the stress resultants in the element coordinate system of each element, see Conventions;

  3. the eigenvalue problem

    \[(\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G) \boldsymbol{u} = \boldsymbol{0}\]

    solved by eigvals, eigvecs = linear_buckling(KC0uu, KGuu, num_eigvalues=...).

\(\boldsymbol{K}_G\) is linear in \(\boldsymbol{u}_0\), hence in the reference load, and the critical load is \(\lambda_1 \boldsymbol{F}_{ref}\), with \(\lambda_1\) = eigvals[0]. linear_buckling() returns the positive load multipliers first, in increasing order, followed, when fewer than num_eigvalues positive ones exist, by the negative ones closest to zero, which are the critical load multipliers of the reversed reference load. A load that is not scaled by \(\lambda\), e.g. a preload with the geometric stiffness matrix \(\boldsymbol{K}_{G0}\), belongs to the first matrix, linear_buckling(KC0uu + KG0uu, KGuu), which is also how the critical value of one load component is found under the other components kept constant.

The beams have update_KG for the axial force of the prestress, see tests/test_beam_linear_buckling_3D_column.py; Truss and Spring have no geometric stiffness matrix.

Natural frequencies#

The frequency problem

\[(\boldsymbol{K} - \omega_n^2 \boldsymbol{M}) \boldsymbol{u} = \boldsymbol{0}\]

is solved by omegan, eigvecs = natural_frequency(Kuu, Muu, num_eigvalues=...), which returns the circular frequencies \(\omega_n\), in rad/s with consistent units, in increasing order. \(\boldsymbol{K}\) is \(\boldsymbol{K}_{C0}\), plus the geometric stiffness matrix of a prestress when it applies, e.g. KC0uu + KGuu, see tests/test_quad4r_natural_freq_pre_stress.py. The mass matrix is assembled with update_M, whose mtype selects the mass matrix of each element, see Choosing an element.

With an unsymmetric \(\boldsymbol{K}\), e.g. with the aerodynamic stiffness matrix of the piston theory or the load stiffness of a follower load, the frequencies may become complex; omegan is then the principal square root of the complex eigenvalues. Two frequencies that coalesce and become complex indicate the loss of stability by flutter, which linear_buckling() does not detect: this is the dynamic method of Bolotin (1963), see natural_frequency().

natural_frequency() solves the undamped problem.

Damped frequencies#

With a damping matrix \(\boldsymbol{C}\), the free vibrations \(\boldsymbol{u} e^{\lambda t}\) give the quadratic eigenvalue problem

\[(\lambda^2 \boldsymbol{M} + \lambda \boldsymbol{C} + \boldsymbol{K}) \boldsymbol{u} = \boldsymbol{0}\]

solved by eigvals, eigvecs = complex_frequency(Kuu, Cuu, Muu, num_eigvalues=...), which returns the complex eigenvalues \(\lambda\), one per complex conjugate pair, with \(Im(\lambda) \ge 0\), and the complex eigenvectors. The damped circular frequency is \(\omega_d = |Im(\lambda)|\), abs(eigvals.imag), the damping ratio \(\zeta = -Re(\lambda)/|\lambda|\), -eigvals.real/abs(eigvals), and a positive \(Re(\lambda)\) is a growing motion. With \(\boldsymbol{C} = \boldsymbol{0}\), \(\lambda = i \omega_n\), with the \(\omega_n\) of natural_frequency().

The damping matrix is assembled by the user, e.g. Rayleigh damping \(\boldsymbol{C} = \alpha \boldsymbol{M} + \beta \boldsymbol{K}\), or the aerodynamic damping of the piston theory, Cuu = -mu*CAuu with update_CA of Quad4 and Quad4R, with the aerodynamic stiffness in Kuu = KC0uu - beta*KA_betauu - gamma*KA_gammauu, see update_KA_beta. A flutter analysis then follows the largest \(Re(\lambda)\) over the flow speed: the aerodynamic damping moves the flutter onset beyond the coalescence of the frequencies of natural_frequency(), see tests/test_quad4_piston_theory.py.

Verification and dense fallback#

The three functions check that \(\boldsymbol{K}\) has no null row or column, a dof without stiffness, which raises a numpy.linalg.LinAlgError naming the dofs, see check_null_rows_cols(), whereas \(\boldsymbol{K}_G\), \(\boldsymbol{C}\) and \(\boldsymbol{M}\) may have null rows, e.g. a lumped mass matrix without rotary inertia. They equilibrate the diagonal of the matrices with diagonal_preconditioner(), which does not change the eigenvalues, and linear_buckling() and natural_frequency() select the solvers of symmetric or unsymmetric matrices with is_symmetric(). Then:

  • the sparse solvers of SciPy run first, in shift-invert or Cayley mode, with a shift estimated from the matrices for the linear buckling, estimate_cayley_sigma();

  • every result is verified: the relative residual of each eigenpair and, for symmetric matrices, the absence of an eigenvalue missed below the highest one returned, from the count of the negative pivots of a shifted matrix, see check_eigenpairs() and is_positive_definite(); an eigenpair of natural_frequency() or complex_frequency() whose relative residual is too large still passes when its backward error is at the round-off, e.g. an eigenvalue close to zero of a model pre-stressed close to buckling, see pyfe3d.solver.eigen._backward_errors(); for the complex eigenvalues of complex_frequency() no such count is implemented, and only the residuals are verified;

  • when the sparse solvers fail, a dense solver of SciPy is used as a last resort, when the dense problem has at most max_dense_size dofs, with a DenseFallbackWarning giving the reasons; otherwise a RuntimeError lists every failure, such that no unverified result is returned silently;

  • sparse_solver=False solves the full problem with a dense solver, returning all the eigenvalues, as a reference for validations.

The fallback can be turned into an error, e.g. in a test suite:

import warnings
from pyfe3d.solver import DenseFallbackWarning
warnings.simplefilter('error', DenseFallbackWarning)

With Intel MKL 2024.2.0 to 2025.0.0, the symmetric solver of ARPACK returns wrong eigenpairs when its number of Lanczos vectors exceeds 32; these versions are detected and the number is capped, see pyfe3d.solver.arpackutils, also for the override with the environment variable PYFE3D_ARPACK_MAX_NCV.

Examples#

The reference load, a preload and the reversed reference load of a Quad4 plate:

def test_quad4_plate_reference_load():
    r"""The Quad4 plate of ``test_solver.py``: a reference load twice as
    large halves the load multipliers, a constant preload belongs to ``K``,
    and the tensile reference load gives the negative load multipliers
    closest to zero, without the dense fallback"""
    KC0uu, KGuu = plate_matrices(Nxx=-1.)
    with warnings.catch_warnings():
        warnings.simplefilter('error', DenseFallbackWarning)
        eigvals, _ = linear_buckling(KC0uu, KGuu, num_eigvalues=4)
        lambda_1 = eigvals[0]
        eigvals2, _ = linear_buckling(KC0uu, 2*KGuu, num_eigvalues=4)
        eigvals_pre, _ = linear_buckling(KC0uu + 0.4*lambda_1*KGuu, KGuu,
                                         num_eigvalues=4)
        eigvals_rev, eigvecs_rev = linear_buckling(KC0uu, -KGuu,
                                                   num_eigvalues=4)
    assert np.all(eigvals > 0)
    np.testing.assert_allclose(eigvals2, eigvals/2, rtol=1e-8)
    np.testing.assert_allclose(eigvals_pre[0], 0.6*lambda_1, rtol=1e-8)
    np.testing.assert_allclose(eigvals_rev, -eigvals, rtol=1e-8)
    assert np.all(relative_residuals(KC0uu, -KGuu, eigvals_rev, eigvecs_rev)
                  < 1e-6)

The frequencies of the same plate under a compressive prestress, whose lowest one vanishes at the critical load multiplier:

def test_quad4_plate_prestress():
    r"""The Quad4 plate of ``test_solver.py`` under a compressive prestress
    `\lambda \boldsymbol{K}_G`: the frequencies decrease with `\lambda`, and the
    squared lowest one vanishes linearly at the critical load multiplier of
    :func:`pyfe3d.solver.linear_buckling`; with ``symmetric=True`` the
    eigenvectors are normalized with respect to the mass matrix"""
    from pyfe3d.solver import linear_buckling
    KC0uu, KGuu, Muu = plate_matrices(Nxx=-1., with_mass=True)
    lambda_1 = linear_buckling(KC0uu, KGuu, num_eigvalues=2)[0][0]
    omega2 = []
    with warnings.catch_warnings():
        warnings.simplefilter('error', DenseFallbackWarning)
        for frac in [0., 0.5, 0.9, 0.99, 0.999]:
            omegan, eigvecs = natural_frequency(KC0uu + frac*lambda_1*KGuu,
                                                Muu, num_eigvalues=4,
                                                symmetric=True)
            assert np.all(np.diff(omegan) >= 0)
            np.testing.assert_allclose(eigvecs.T @ Muu @ eigvecs, np.eye(4),
                                       atol=1e-8)
            omega2.append(omegan**2)
    omega2 = np.array(omega2)
    assert np.all(np.diff(omega2, axis=0) <= 1e-12*omega2[0, 0])
    np.testing.assert_allclose(omega2[-2, 0]/omega2[-1, 0], 10., rtol=1e-2)

The damped eigenvalues of a chain of masses with Rayleigh damping, against the analytical modal values:

def test_rayleigh_damping_chain():
    r"""Rayleigh damping of a chain of masses: the analytical modal
    eigenvalues and modes"""
    n, k, m = 40, 100., 2.
    K, M = chain(n, k, m)
    alpha, beta = 0.2, 2e-3
    C = alpha*M + beta*K
    j = np.arange(1, 7)
    omega = 2*np.sqrt(k/m)*np.sin(j*np.pi/(2*(n + 1)))
    zeta = alpha/(2*omega) + beta*omega/2
    lam_ref = -zeta*omega + 1j*omega*np.sqrt(1 - zeta**2)
    with warnings.catch_warnings():
        warnings.simplefilter('error', DenseFallbackWarning)
        lam, U = complex_frequency(K, C, M, num_eigvalues=6)
    np.testing.assert_allclose(lam, lam_ref, rtol=1e-10)
    np.testing.assert_allclose(-lam.real/np.abs(lam), zeta, rtol=1e-10)
    np.testing.assert_allclose(np.abs(lam), omega, rtol=1e-10)
    i = np.arange(1, n + 1)
    for col, jj in enumerate(j):
        mode = np.sin(i*jj*np.pi/(n + 1))
        mode /= np.linalg.norm(mode)
        np.testing.assert_allclose(np.abs(mode @ U[:, col]), 1., rtol=1e-10)
    assert np.all(residuals(K, C, M, lam, U) < 1e-10)

Complete models are in Linear buckling analysis of a plate, Linear buckling analysis of a cylinder under torsion, Natural frequency analysis of a plate and Natural frequency analysis of a cylinder.