Follower loads and unsymmetric matrices#
Loads that depend on the configuration, such as follower (hydrostatic)
pressures, are supported through optional keyword arguments of the callables
of Analysis, see callbacks: calc_fint, calc_kC and
calc_kG receive the load factor inc when they declare it, and
calc_fext receives the current configuration c, returning the load
vector \(\partial R/\partial \lambda\) used by the arc-length methods. Their
load stiffness makes the tangent stiffness matrix unsymmetric in general,
which the solvers support: solve() and the non-linear solvers use LU
factorizations, lb() selects the solvers of unsymmetric matrices when
needed, symmetric=None, and freq() uses general eigenvalue solvers
by default.
Non-linear analyses#
Hardening springs under a load \(\lambda(\{q_0\} + [Q]\{c\})\) with an unsymmetric \([Q]\): the callables receive the load factor and the configuration, and Newton-Raphson and the arc-length methods converge to the equilibrium path. The Riks method converges quadratically only with the load vector of the current configuration:
class FollowerProblem(object):
def __init__(self, K, a, q0, Q, with_c=True):
self.K = K
self.a = a
self.q0 = q0
self.Q = Q
self.with_c = with_c
self.an = None
# (inc, ||R||/||lbd*q||, number of converged increments) of each
# calc_fint call
self.calls = []
def fint_el(self, c):
return self.K @ c + self.a*c**3
def residual(self, c, lbd):
return lbd*(self.q0 + self.Q @ c) - self.fint_el(c)
def calc_fext(self, inc=1., silent=True, c=None):
if c is None:
return inc*self.q0
return inc*(self.q0 + self.Q @ c)
def calc_fint(self, c, inc=1., silent=True):
R = self.residual(c, inc)
ref = max(np.linalg.norm(inc*(self.q0 + self.Q @ c)), 1e-300)
step = len(self.an.increments) if self.an is not None else 0
self.calls.append((inc, np.linalg.norm(R)/ref, step))
return self.fint_el(c) - inc*(self.Q @ c)
def calc_kC(self, c=None, NLgeom=False, inc=1., silent=True):
if not NLgeom:
return csr_matrix(self.K)
return csr_matrix(self.K + np.diag(3*self.a*c**2) - inc*self.Q)
def calc_kG(self, c=None, NLgeom=False, silent=True):
return csr_matrix(self.K.shape)
def analysis(self, **kwargs):
if self.with_c:
calc_fext = self.calc_fext
else:
calc_fext = lambda inc=1., silent=True: self.calc_fext(inc)
an = Analysis(calc_fext, self.calc_fint, self.calc_kC, self.calc_kG)
for k, v in kwargs.items():
setattr(an, k, v)
self.an = an
return an
def test_riks_quadratic_only_with_configuration_load_vector():
r"""The bordered system of Riks is the Newton-Raphson method of the
augmented system only when dR/dlbd = q0 + Q c, i.e. with ``calc_fext``
accepting ``c``; with the constant reference vector the corrections
converge linearly"""
worst = {}
for with_c in (True, False):
pb = make_problem(with_c=with_c)
an = pb.analysis(NL_method='arc_length_riks', initialInc=0.5,
relTOL=1e-13, maxNumIter=60)
an.static(NLgeom=True, silent=True)
assert np.isclose(an.increments[-1], 1.)
# residuals of the iterations of the first step, whose first value
# is the one of the tangent predictor
first = [e for inc, e, step in pb.calls if step == 0]
worst[with_c] = (min(orders(first)) if orders(first) else None,
first)
p_c, errors_c = worst[True]
assert p_c is not None and p_c > 1.7, errors_c
p_const, errors_const = worst[False]
assert len(errors_const) > len(errors_c), (errors_c, errors_const)
Linear static analysis#
With a configuration-dependent load, Analysis.static(NLgeom=False) solves
the geometrically linear problem \([K_0]\{c\} = \{F(c)\}\):
def test_analysis_linear_static_with_configuration_dependent_load():
pb = linear_follower_problem()
an = pb.analysis()
an.relTOL = 1e-12
increments, cs = an.static(NLgeom=False, silent=True)
ref = np.linalg.solve(pb.K0 - pb.Q, pb.F0)
assert np.allclose(cs[0], ref, rtol=1e-9)
# without c the load is the reference one, as before
pb_old = linear_follower_problem(with_c=False)
increments, cs = pb_old.analysis().static(NLgeom=False, silent=True)
assert np.allclose(cs[0], np.linalg.solve(pb.K0, pb.F0))
Linear buckling with unsymmetric matrices#
@pytest.mark.parametrize('sparse_solver,max_dense_size',
[(True, 2000), (True, 0), (False, 2000)])
def test_lb_unsymmetric_real(sparse_solver, max_dense_size):
r"""Unsymmetric KG with real eigenvalues, against scipy.linalg.eig"""
rng = np.random.default_rng(3)
n = 40
K = spd(n, rng)
# KG = S + small skew part keeps the eigenvalues real
S = -spd(n, rng)/n
W = rng.standard_normal((n, n))
KG = S + 1e-3*(W - W.T)
assert not is_symmetric(KG)
eigvals, eigvecs = lb(csr_matrix(K), csr_matrix(KG), silent=True,
sparse_solver=sparse_solver, num_eigvalues=6,
max_dense_size=max_dense_size)
assert np.isrealobj(eigvals) and np.isrealobj(eigvecs)
mu = eig(KG, K, right=False)
lam = np.sort((-1/mu).real)
lam = lam[lam > 0]
assert np.allclose(eigvals[:5], lam[:5], rtol=1e-8)
for i in range(5):
u = eigvecs[:, i]
r = K @ u + eigvals[i]*(KG @ u)
assert np.linalg.norm(r) <= 1e-8*np.linalg.norm(K @ u)
Under follower loads, many negative or complex load multipliers may have a
smaller modulus than the critical one, e.g. when the pressure is reversed.
The sparse solver of lb() searches the lowest positive real load
multiplier along the positive real axis, see Eigenvalue solvers. The
cantilevers of Schweizerhof and Ramm (1984), with the matrices of
panels, original and reversed:
@pytest.mark.parametrize('name, sign, approx', FOLLOWER)
def test_lb_follower_shift_search(name, sign, approx):
"""Regression test: the sparse fallback of version 0.5.3, eigs() on
K^-1 KG, silently returned only negative or complex load multipliers for
the reversed problems. The shift search returns the lowest positive real
load multiplier of the full dense eig(), first, without the dense
fallback"""
import warnings
from structsolve.linear_buckling import DenseFallbackWarning
K, KG = load_follower(name, sign)
ref = lowest_positive_real(lb(K, KG, silent=True, sparse_solver=False)[0])
np.testing.assert_allclose(ref, approx, rtol=1e-5)
with warnings.catch_warnings():
warnings.simplefilter('error', DenseFallbackWarning)
eigvals, eigvecs = lb(K, KG, silent=True, max_dense_size=0)
assert eigvals.shape == (25,) and eigvecs.shape == (K.shape[0], 25)
assert np.isreal(eigvals[0]) and eigvals[0].real > 0
np.testing.assert_allclose(eigvals[0].real, ref, rtol=1e-8)
assert lowest_positive_real(eigvals) == eigvals[0].real
res = relative_residuals(K, KG, eigvals, eigvecs)
assert res[0] < 1e-8 and np.all(res < 1e-5)
Flutter: the kinetic criterion#
Ziegler’s double pendulum under a follower force never diverges, its
stiffness matrix is never singular, and lb() finds no critical load.
It loses stability by flutter at \(P_f = (7/2 - \sqrt{2}) k/l\), where two
natural frequencies of freq() coalesce and become complex:
def ziegler(P, extra=0):
r"""``K(P)`` and ``M``, with ``extra`` decoupled unit oscillators added
such that the sparse solvers, limited to ``k < size - 1``, can be used"""
k, l, m = K_SPRING, L_BAR, MASS
K = np.array([[2*k - P*l, -k + P*l], [-k, k]])
M = m*l**2*np.array([[3., 1.], [1., 1.]])
if extra:
rng = np.random.default_rng(0)
K = block_diag(K, np.diag(1e3 + rng.random(extra)))
M = block_diag(M, np.eye(extra))
return K, M
@pytest.mark.parametrize('sparse_solver,extra', [(False, 0), (True, 10)])
def test_freq_ziegler_flutter(sparse_solver, extra):
r"""The general solver of :func:`freq` gives real frequencies below the
flutter load and complex ones above it; the flutter load found by
bisection is the analytical one"""
kw = dict(silent=True, sparse_solver=sparse_solver, num_eigvalues=4,
sort=False)
below = freq(*map(csr_matrix, ziegler(0.9*P_FLUTTER, extra)), **kw)[0]
above = freq(*map(csr_matrix, ziegler(1.1*P_FLUTTER, extra)), **kw)[0]
assert max_imag(below) <= 1e-8*np.abs(below).max()
assert max_imag(above) > 1e-3*np.abs(above).max()
a, b = 0.9*P_FLUTTER, 1.1*P_FLUTTER
for _ in range(40):
c = 0.5*(a + b)
lam2 = freq(*map(csr_matrix, ziegler(c, extra)), **kw)[0]
if max_imag(lam2) > 1e-6*np.abs(lam2).max():
b = c
else:
a = c
assert np.isclose(0.5*(a + b), P_FLUTTER, rtol=1e-4)
The follower pressure loads of panels
use all these features, see Shell.add_pressure_load(..., follower=True).