Practice, FSDT semi-analytical plate buckling analysis via the Ritz Method¶

This notebook implements a semi-analytical framework to determine the critical buckling load of plate using the first-order shear deformation theory (FSDT). It utilizes the Ritz method with Legendre polynomials as basis functions.

import numpy as np
from scipy.special import roots_legendre

from composites import isotropic_plate

from buckling.legendre import vecf, vecfxi

# approximation order
m1 = 20
m2 = 10
N = 3*m1*m2
print('m1 = %d' % m1)
print('m2 = %d' % m2)
print('N = %d' % N)
i = np.arange(m1)
j = np.arange(m2)
pts1, weights1 = roots_legendre(2*m1 - 1)
pts2, weights2 = roots_legendre(2*m2 - 1)

# Material properties
E = 200.e9
nu = 0.3
G = E/(2*(1 + nu))

# Boundary conditions w
wxit1 = 0
wxir1 = 1
wxit2 = 0
wxir2 = 1
wetat1 = 0
wetar1 = 1
wetat2 = 0
wetar2 = 1

# Boundary conditions phix phiy
xit1 = 1
xir1 = 1
xit2 = 1
xir2 = 1
etat1 = 1
etar1 = 1
etat2 = 1
etar2 = 1

# Geometric properties
a = 0.3
b = 0.1
h = 0.003

# Calculating ABD matrices
prop = isotropic_plate(thickness=h, E=E, nu=nu)

# Initial stress state
Nxxhat = -100.
Nyyhat = 0
Nxyhat = 0


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

# stiffness matrix
# numerical integration in 2D using Legendre-Gauss quadrature
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
        
        # constitutive stiffness matrix
        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)

        # geometric stiffness matrix
        addouter(KG, detJ*weight*Nxxhat*Swx, Swx)
        addouter(KG, detJ*weight*Nyyhat*Swy, Swy)
        addouter(KG, detJ*weight*Nxyhat*Swx, Swy)
        addouter(KG, detJ*weight*Nxyhat*Swy, Swx)

m1 = 20
m2 = 10
N = 600

1Eigenvalue problem solver¶

We solve the generalized eigenvalue problem $(\mathbf{K} - \lambda \mathbf{K}_G) \boldsymbol{\phi} = \mathbf{0}$ . The output $\lambda_{cr}$ represents the buckling load multiplier. The resulting eigenvalue is normalized against the classical plate buckling parameter $k_x = \frac{N_{xx} b^2}{\pi^2 D}$ for verification.

from scipy.sparse import csc_matrix
from structsolve import lb

eigvals, eigvecs = lb(csc_matrix(K), csc_matrix(KG))

Running linear buckling analysis...
		Eigenvalue solver... 
			eigsh() solver...
				WARNING: Factor is exactly singular
			aborted!
			Removing null columns...
				56 columns removed
			finished!
			eigsh() solver...
			finished!
		finished!
	first 5 eigenvalues:
		19225.238339319305
		20838.118452881405
		22626.041437339274
		24651.100960083982
		29966.53310487754

2Checking results¶

from composites.kassapoglou import calc_Nxx_crit

Nxx_cr = abs(Nxxhat*eigvals[0])

m = None
n = 1 # assuming one half-wave along the y direction (plate width)
Nxx_cr_ref = calc_Nxx_crit(a, b, m, n, prop.D11, prop.D12, prop.D22, prop.D66)

print('Nxx_cr', Nxx_cr)
print('Nxx_cr_ref', Nxx_cr_ref)
print('Nxx_cr/Nxx_cr_ref', Nxx_cr/Nxx_cr_ref)

Nxx_cr 1922523.8339319306
Nxx_cr_ref 1952229.4419737186
Nxx_cr/Nxx_cr_ref 0.9847837516415308

33D FSDT buckling mode visualization algorithm¶

This algorithm renders the 3D deformed shape of a plate based on First-Order Shear Deformation Theory (FSDT) kinematics. It isolates the midplane degrees of freedom and extrapolates them through the thickness for a full volumetric wireframe visualization.

3.11D Domain Mapping¶

The physical dimensions ( $a$ , $b$ ) are discretized into 1D vectors and linearly mapped to the natural parametric coordinates $\xi, \eta \in [-1, 1]$ .
The thickness $h$ is mapped to a discrete $z$ -vector $z \in [-h/2, h/2]$ .

3.2Basis Function Segregation¶

FSDT requires distinct boundary conditions for translation and rotation.
The script evaluates the 1D Legendre polynomial basis functions (vecf) twice: once for the transverse displacement $w$ (using wxit and wxir indices) and once for the rotations $\phi_x, \phi_y$ (using xit and xir indices).

3.3Midplane Tensor Contraction & 3D Extrapolation¶

Modal Extraction: The contiguous global eigenvector is sliced into its respective midplane components: $c_w, c_{\phi x}, c_{\phi y}$ of shape (m1, m2).
2D Assembly: np.einsum computes the continuous 2D midplane spatial fields ( $w(x,y)$ , $\phi_x(x,y)$ , $\phi_y(x,y)$ ) via tensor contraction.
Through-Thickness Mapping: The 3D spatial field is extrapolated over the $z$ -vector using FSDT continuous field assumptions: $U(x,y,z) = z \cdot \phi_x(x,y)$ $V(x,y,z) = z \cdot \phi_y(x,y)$ $W(x,y,z) = w(x,y)$

3.4Wireframe Topology & Rendering¶

The algorithm scales the 3D displacement arrays and perturbs the Cartesian grid coordinates.
The grid is flattened into 1D lists separated by None to prevent Plotly from rendering artificial connecting lines across the mesh boundaries.
To bypass Plotly’s strict type-checking for contour colors, a dummy float 0.0 is injected into the out-of-plane $W$ amplitude array at the None breakpoints.

3.5Interactive UI¶

A toggle widget (updatemenus) is generated to pass boolean visibility masks to the figure traces, allowing instantaneous switching between precomputed buckling modes without re-executing the Python kernel.

# Create a new cell at the bottom of your notebook and paste this code

import numpy as np
import plotly.graph_objects as go

# --- Extrapolating FSDT 2D Midplane to 3D Space ---
num_modes = 3  
scale = 5  # Adjusted scale; modify based on the a, b, h aspect ratio

nx, ny, nz = 40, 20, 3

# 1. Define 1D Coordinate Arrays
x_1d = np.linspace(0, a, nx)
y_1d = np.linspace(0, b, ny)
z_1d = np.linspace(-h/2, h/2, nz)

# Map to natural coordinates [-1, 1]
xi_1d = 2 * x_1d / a - 1
eta_1d = 2 * y_1d / b - 1

# 2. Evaluate 1D Basis Functions
# For w DOFs (transverse displacement), we use the w-specific boundary condition indices
wP_xi = np.array([vecf(m1, xi, wxit1, wxir1, wxit2, wxir2) for xi in xi_1d])
wP_eta = np.array([vecf(m2, eta, wetat1, wetar1, wetat2, wetar2) for eta in eta_1d])

# For phi_x and phi_y DOFs (rotations), we use the phi-specific boundary condition indices
P_xi = np.array([vecf(m1, xi, xit1, xir1, xit2, xir2) for xi in xi_1d])
P_eta = np.array([vecf(m2, eta, etat1, etar1, etat2, etar2) for eta in eta_1d])

U = np.zeros((num_modes, nx, ny, nz))
V = np.zeros((num_modes, nx, ny, nz))
W = np.zeros((num_modes, nx, ny, nz))

n_dof = m1 * m2

# 3. Assemble Fields via FSDT Kinematics
for mode in range(num_modes):
    
    # Extract DOF blocks. The script assumes DOFs are packed as: [w, phi_x, phi_y]
    c_w = eigvecs[0 : n_dof, mode].reshape((m1, m2))
    c_phix = eigvecs[n_dof : 2*n_dof, mode].reshape((m1, m2))
    c_phiy = eigvecs[2*n_dof : 3*n_dof, mode].reshape((m1, m2))
    
    # Reconstruct the 2D midplane continuous fields using tensor contraction
    w_2d = np.einsum('ij, xi, yj -> xy', c_w, wP_xi, wP_eta)
    phix_2d = np.einsum('ij, xi, yj -> xy', c_phix, P_xi, P_eta)
    phiy_2d = np.einsum('ij, xi, yj -> xy', c_phiy, P_xi, P_eta)
    
    # Extrapolate to 3D via FSDT assumptions: U = z*phi_x, V = z*phi_y, W = w
    for k, z in enumerate(z_1d):
        U[mode, :, :, k] = z * phix_2d
        V[mode, :, :, k] = z * phiy_2d
        W[mode, :, :, k] = w_2d

X, Y, Z = np.meshgrid(x_1d, y_1d, z_1d, indexing='ij')

fig = go.Figure()

# 4. Build wireframe traces
for mode in range(num_modes):
    X_def = X + scale * U[mode]
    Y_def = Y + scale * V[mode]
    Z_def = Z + scale * W[mode]
    W_amp = W[mode]
    
    lines_x, lines_y, lines_z, lines_w = [], [], [], []

    # X-direction lines
    for j in range(ny):
        for k in range(nz):
            lines_x.extend(X_def[:, j, k].tolist() + [None])
            lines_y.extend(Y_def[:, j, k].tolist() + [None])
            lines_z.extend(Z_def[:, j, k].tolist() + [None])
            lines_w.extend(W_amp[:, j, k].tolist() + [0.0])

    # Y-direction lines
    for i in range(nx):
        for k in range(nz):
            lines_x.extend(X_def[i, :, k].tolist() + [None])
            lines_y.extend(Y_def[i, :, k].tolist() + [None])
            lines_z.extend(Z_def[i, :, k].tolist() + [None])
            lines_w.extend(W_amp[i, :, k].tolist() + [0.0])

    # Z-direction lines
    for i in range(nx):
        for j in range(ny):
            lines_x.extend(X_def[i, j, :].tolist() + [None])
            lines_y.extend(Y_def[i, j, :].tolist() + [None])
            lines_z.extend(Z_def[i, j, :].tolist() + [None])
            lines_w.extend(W_amp[i, j, :].tolist() + [0.0])

    fig.add_trace(go.Scatter3d(
        x=lines_x, y=lines_y, z=lines_z, mode='lines',
        line=dict(color=lines_w, colorscale='Jet', width=3, showscale=True, colorbar=dict(title=f"W (Mode {mode+1})", x=0.85)),
        name=f'Mode {mode+1}', opacity=0.9, visible=(mode == 0)
    ))

# 5. Toggle Widget
buttons = []
for mode in range(num_modes):
    visibility = [False] * num_modes
    visibility[mode] = True
    buttons.append(dict(
        label=f'Mode {mode+1}', method='update',
        args=[{'visible': visibility}, {'title': f'FSDT Buckling Mode Shape: Mode {mode+1}'}]
    ))

fig.update_layout(
    updatemenus=[dict(active=0, buttons=buttons, x=0.05, xanchor='left', y=1.05, yanchor='top', direction='down', showactive=True)],
    scene=dict(xaxis_title='X', yaxis_title='Y', zaxis_title='Z', aspectmode='data'),
    title='FSDT Buckling Mode Shape: Mode 1', showlegend=False, margin=dict(l=0, r=0, b=0, t=40)
)

fig.show()