Linear buckling analysis#
The function lb() solves the linear buckling eigenvalue problem
\(([K] + \lambda [K_G])\{u\} = \{0\}\), returning the load multipliers
\(\lambda\).
The example below calculates the critical load of a simply supported
isotropic plate under uniaxial compression, using a Ritz model based on the
first-order shear deformation theory (FSDT). The approximation functions are
the hierarchical Legendre functions of Bardell, the same as those of
buckling, and the plate
stiffnesses are those of
composites; both are
implemented with NumPy in the module tests/ritz.py:
"""Approximation functions and plate properties of the Ritz models of the tests
They replace ``buckling.legendre`` and ``composites.isotropic_plate``, such
that the tests only depend on NumPy and SciPy.
"""
from functools import lru_cache
from types import SimpleNamespace
import numpy as np
from numpy.polynomial import legendre
#: Power series coefficients of the cubic Hermite functions, which control
#: the displacement and the rotation at each extremity
_HERMITE = np.array([[0.5, -0.75, 0., 0.25],
[0.125, -0.125, -0.125, 0.125],
[0.5, 0.75, 0., -0.25],
[-0.125, -0.125, 0.125, 0.125]])
@lru_cache(maxsize=None)
def _coefficients(n, d1t, d1r, d2t, d2r, der):
"""Legendre series coefficients, one column per function"""
C = np.zeros((max(n, 4), max(n, 4)))
for i, flag in enumerate((d1t, d1r, d2t, d2r)):
C[:4, i] = flag*legendre.poly2leg(_HERMITE[i])
for i in range(4, n):
Pi_2 = np.zeros(i - 1)
Pi_2[-1] = 1.
C[:i+1, i] = legendre.legint(Pi_2, m=2, lbnd=-1)
C = legendre.legder(C[:, :n], m=der, axis=0) if der else C[:, :n]
C.flags.writeable = False
return C
def vecf(n, xi, d1t, d1r, d2t, d2r):
r"""Hierarchical functions of Bardell at the natural coordinate ``xi``
The first four are cubic Hermite functions, whose displacement and
rotation at `\xi=-1` and `\xi=+1` are controlled by the flags ``d1t``,
``d1r``, ``d2t`` and ``d2r``. The following ones are the Legendre
polynomials `P_{i-2}` integrated twice from `-1`, which vanish with
their first derivative at both extremities. Same as
``buckling.legendre.vecf``.
"""
return legendre.legval(xi, _coefficients(n, d1t, d1r, d2t, d2r, 0))
def vecfxi(n, xi, d1t, d1r, d2t, d2r):
"""First derivative of :func:`vecf` with respect to ``xi``"""
return legendre.legval(xi, _coefficients(n, d1t, d1r, d2t, d2r, 1))
def vecfxixi(n, xi, d1t, d1r, d2t, d2r):
"""Second derivative of :func:`vecf` with respect to ``xi``"""
return legendre.legval(xi, _coefficients(n, d1t, d1r, d2t, d2r, 2))
def isotropic_plate(thickness, E, nu):
"""Bending and transverse shear stiffnesses of an isotropic plate
The transverse shear stiffnesses use the shear correction factor 5/6.
Same as the ``D`` and the transverse shear terms of the ``Laminate``
returned by ``composites.isotropic_plate``.
"""
D = E*thickness**3/(12*(1 - nu**2))
G = E/(2*(1 + nu))
return SimpleNamespace(
D11=D, D12=nu*D, D16=0., D22=D, D26=0., D66=(1 - nu)*D/2,
A44=5/6*G*thickness, A45=0., A55=5/6*G*thickness)
The stiffness matrices are integrated numerically and the critical load is
compared with the analytical solution. The code is extracted from one of the
structsolve unit tests:
def test_lb_plate_buckling_fsdt():
"""Test plate buckling (FSDT) based on semi-analytical Ritz method.
Reference: notebook BucklingPlates-FSDT.ipynb from the buckling repository.
Uses Legendre polynomials as basis functions with Gauss-Legendre quadrature.
"""
from scipy.special import roots_legendre
m1 = 20
m2 = 10
N = 3 * m1 * m2
pts1, weights1 = roots_legendre(2 * m1 - 1)
pts2, weights2 = roots_legendre(2 * m2 - 1)
E = 200.e9
nu = 0.3
# BCs: w simply supported, phi free
wxit1, wxir1, wxit2, wxir2 = 0, 1, 0, 1
wetat1, wetar1, wetat2, wetar2 = 0, 1, 0, 1
xit1, xir1, xit2, xir2 = 1, 1, 1, 1
etat1, etar1, etat2, etar2 = 1, 1, 1, 1
a, b, h = 0.3, 0.1, 0.003
prop = isotropic_plate(thickness=h, E=E, nu=nu)
Nxxhat = -100.
Swx = np.zeros(N)
Swy = np.zeros(N)
Sphix = np.zeros(N)
Sphiy = np.zeros(N)
Sphixx = np.zeros(N)
Sphixy = np.zeros(N)
Sphiyx = np.zeros(N)
Sphiyy = np.zeros(N)
buff = np.zeros((N, N))
K = np.zeros((N, N))
KG = np.zeros((N, N))
def addouter(matrix, vec1, vec2):
np.outer(vec1, vec2, out=buff)
matrix += buff
for xi, wxi in zip(pts1, weights1):
wP_xi = vecf(m1, xi, wxit1, wxir1, wxit2, wxir2)
wPx_xi = vecfxi(m1, xi, wxit1, wxir1, wxit2, wxir2)
P_xi = vecf(m1, xi, xit1, xir1, xit2, xir2)
Px_xi = vecfxi(m1, xi, xit1, xir1, xit2, xir2)
for eta, weta in zip(pts2, weights2):
wP_eta = vecf(m2, eta, wetat1, wetar1, wetat2, wetar2)
wPx_eta = vecfxi(m2, eta, wetat1, wetar1, wetat2, wetar2)
P_eta = vecf(m2, eta, etat1, etar1, etat2, etar2)
Px_eta = vecfxi(m2, eta, etat1, etar1, etat2, etar2)
weight = wxi * weta
Pi, Pj = np.meshgrid(P_xi, P_eta, indexing='ij')
Sphix[m1*m2:2*m1*m2] = (Pi * Pj).flatten()
Sphiy[2*m1*m2:] = (Pi * Pj).flatten()
Pxi, Pj = np.meshgrid(wPx_xi, wP_eta, indexing='ij')
Swx[:m1*m2] = (Pxi * Pj * (2 / a)).flatten()
Pxi, Pj = np.meshgrid(Px_xi, P_eta, indexing='ij')
Sphixx[m1*m2:2*m1*m2] = (Pxi * Pj * (2 / a)).flatten()
Sphiyx[2*m1*m2:] = (Pxi * Pj * (2 / a)).flatten()
Pi, Pxj = np.meshgrid(wP_xi, wPx_eta, indexing='ij')
Swy[:m1*m2] = (Pi * Pxj * (2 / b)).flatten()
Pi, Pxj = np.meshgrid(P_xi, Px_eta, indexing='ij')
Sphixy[m1*m2:2*m1*m2] = (Pi * Pxj * (2 / b)).flatten()
Sphiyy[2*m1*m2:] = (Pi * Pxj * (2 / b)).flatten()
e1xx = Sphixx
e1yy = Sphiyy
e1xy = Sphixy + Sphiyx
g0yz = Sphiy + Swy
g0xz = Sphix + Swx
Mxx = prop.D11 * e1xx + prop.D12 * e1yy + prop.D16 * e1xy
Myy = prop.D12 * e1xx + prop.D22 * e1yy + prop.D26 * e1xy
Mxy = prop.D16 * e1xx + prop.D26 * e1yy + prop.D66 * e1xy
Qy = prop.A44 * g0yz + prop.A45 * g0xz
Qx = prop.A45 * g0yz + prop.A55 * g0xz
detJ = a * b / 4
addouter(K, detJ * weight * Mxx, e1xx)
addouter(K, detJ * weight * Myy, e1yy)
addouter(K, detJ * weight * Mxy, e1xy)
addouter(K, detJ * weight * Qy, g0yz)
addouter(K, detJ * weight * Qx, g0xz)
addouter(KG, detJ * weight * Nxxhat * Swx, Swx)
eigvals, eigvecs = lb(csc_matrix(K), csc_matrix(KG), silent=True)
D = E * h**3 / (12 * (1 - nu**2))
Ncr_analytical = 4 * np.pi**2 * D / b**2
Ncr_computed = eigvals[0] * abs(Nxxhat)
# FSDT gives slightly lower Ncr than CLPT due to shear deformation
np.testing.assert_allclose(Ncr_computed, Ncr_analytical, rtol=0.02)