Linear buckling analysis#
The function lb() solves the linear buckling eigenvalue problem
\((\boldsymbol{K} + \lambda \boldsymbol{K}_G) \boldsymbol{u} = \boldsymbol{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, 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
The Legendre-based approximation functions of Castro and Donadon (2017) and
a NumPy version of ``composites.isotropic_plate``, such that the tests only
depend on NumPy and SciPy.
Reference: Saullo G. P. Castro, Maurício V. Donadon. "Assembly of
Semi-Analytical models to Address Linear Buckling and Vibration of Stiffened
Composite Panels with Debonding Defect". Composite Structures, 2017.
https://doi.org/10.1016/j.compstruct.2016.10.026
"""
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.
"""
return legendre.legval(xi, _coefficients(n, d1t, d1r, d2t, d2r, 0))
def vecf_xi(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 vecf_xixi(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.
Uses Legendre polynomials as basis functions with Gauss-Legendre quadrature.
Compares the critical load with the classical solution ``4*pi**2*D/b**2``,
with a relative tolerance of 2 %.
"""
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.
Sw_x = np.zeros(N)
Sw_y = np.zeros(N)
Sphix = np.zeros(N)
Sphiy = np.zeros(N)
Sphix_x = np.zeros(N)
Sphix_y = np.zeros(N)
Sphiy_x = np.zeros(N)
Sphiy_y = 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):
wPxi = vecf(m1, xi, wxit1, wxir1, wxit2, wxir2)
wPxi_xi = vecf_xi(m1, xi, wxit1, wxir1, wxit2, wxir2)
Pxi = vecf(m1, xi, xit1, xir1, xit2, xir2)
Pxi_xi = vecf_xi(m1, xi, xit1, xir1, xit2, xir2)
for eta, weta in zip(pts2, weights2):
wPeta = vecf(m2, eta, wetat1, wetar1, wetat2, wetar2)
wPeta_eta = vecf_xi(m2, eta, wetat1, wetar1, wetat2, wetar2)
Peta = vecf(m2, eta, etat1, etar1, etat2, etar2)
Peta_eta = vecf_xi(m2, eta, etat1, etar1, etat2, etar2)
weight = wxi * weta
Pi, Pj = np.meshgrid(Pxi, Peta, indexing='ij')
Sphix[m1*m2:2*m1*m2] = (Pi * Pj).flatten()
Sphiy[2*m1*m2:] = (Pi * Pj).flatten()
Pi_xi, Pj = np.meshgrid(wPxi_xi, wPeta, indexing='ij')
Sw_x[:m1*m2] = (Pi_xi * Pj * (2 / a)).flatten()
Pi_xi, Pj = np.meshgrid(Pxi_xi, Peta, indexing='ij')
Sphix_x[m1*m2:2*m1*m2] = (Pi_xi * Pj * (2 / a)).flatten()
Sphiy_x[2*m1*m2:] = (Pi_xi * Pj * (2 / a)).flatten()
Pi, Pj_eta = np.meshgrid(wPxi, wPeta_eta, indexing='ij')
Sw_y[:m1*m2] = (Pi * Pj_eta * (2 / b)).flatten()
Pi, Pj_eta = np.meshgrid(Pxi, Peta_eta, indexing='ij')
Sphix_y[m1*m2:2*m1*m2] = (Pi * Pj_eta * (2 / b)).flatten()
Sphiy_y[2*m1*m2:] = (Pi * Pj_eta * (2 / b)).flatten()
e1xx = Sphix_x
e1yy = Sphiy_y
e1xy = Sphix_y + Sphiy_x
g0yz = Sphiy + Sw_y
g0xz = Sphix + Sw_x
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 * Sw_x, Sw_x)
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)