Practice, deflection of a plate using the CLPT¶

!python -m pip install buckling composites structsolve plotly > tmp.txt

import numpy as np
from scipy.special import roots_legendre
from composites import isotropic_plate

from buckling.legendre import vecf, vecfxi, vecfxixi

# approximation order
m1 = m2 = 20
N = 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 (simply supported)
xit1 = 0
xir1 = 1
xit2 = 0
xir2 = 1
etat1 = 0
etar1 = 1
etat2 = 0
etar2 = 1

# Geometric properties
a = 3
b = 7
h = 0.005

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

# Pre-allocate K matrix
K = np.zeros((N, N))

# Constant Jacobian determinant
detJ = a * b / 4

# Stiffness matrix integration using Legendre-Gauss quadrature
for xi, wxi in zip(pts1, weights1):
    P_xi = vecf(m1, xi, xit1, xir1, xit2, xir2)
    Px_xi = vecfxi(m1, xi, xit1, xir1, xit2, xir2)
    Pxx_xi = vecfxixi(m1, xi, xit1, xir1, xit2, xir2)
    
    for eta, weta in zip(pts2, weights2):
        P_eta = vecf(m2, eta, etat1, etar1, etat2, etar2)
        Px_eta = vecfxi(m2, eta, etat1, etar1, etat2, etar2)
        Pxx_eta = vecfxixi(m2, eta, etat1, etar1, etat2, etar2)

        # Vectorized outer products using broadcasting
        Sw   = (P_xi[:, None]   * P_eta[None, :]).ravel()
        Swx  = (Px_xi[:, None]  * P_eta[None, :]   * (2/a)).ravel()
        Swy  = (P_xi[:, None]   * Px_eta[None, :]  * (2/b)).ravel()
        Swxx = (Pxx_xi[:, None] * P_eta[None, :]   * (2/a)**2).ravel()
        Swyy = (P_xi[:, None]   * Pxx_eta[None, :] * (2/b)**2).ravel()
        Swxy = (Px_xi[:, None]  * Px_eta[None, :]  * (2/a) * (2/b)).ravel()
        
        # Kinematics
        e1xx = -Swxx
        e1yy = -Swyy
        e1xy = -2 * Swxy
        
        # Constitutive relations (Bending moments)
        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
        
        weight = wxi * weta

        # Matrix assembly
        K += np.outer(detJ*weight*Mxx, e1xx)
        K += np.outer(detJ*weight*Myy, e1yy)
        K += np.outer(detJ*weight*Mxy, e1xy)

m1 = 20
m2 = 20
N = 400

1External force vector¶

xi = 0 # x = a/2
eta = 0 # y = b/2
P_xi = vecf(m1, xi, xit1, xir1, xit2, xir2)
P_eta = vecf(m2, eta, etat1, etar1, etat2, etar2)
Sw   = (P_xi[:, None]   * P_eta[None, :]).ravel()

Pforce = 1.

Fext = Pforce*Sw

2Solve the linear system Ku = Fext for Ritz coefficients u¶

from scipy.sparse import csc_matrix
from structsolve import solve

u = solve(csc_matrix(K), Fext)

			Removing null columns...
				76 columns removed
			finished!

3Calculating displacement field for plotting¶

nx, ny = 40, 20
X, Y = np.meshgrid(
    np.linspace(0, a, nx),
    np.linspace(0, b, ny),
    indexing='ij'
)

W = np.zeros_like(X)

for ij, x in np.ndenumerate(X):
    xi = 2*x/a - 1
    eta = 2*Y[ij]/b - 1

    P_xi = vecf(m1, xi, xit1, xir1, xit2, xir2)
    P_eta = vecf(m2, eta, etat1, etar1, etat2, etar2)

    Sw = (P_xi[:, None]   * P_eta[None, :]).ravel()

    W[ij] = Sw @ u

4Printing the maximum displacement¶

The FE reference comes from here tests/test_quad4_static_point_load.py

print('W.min(), W.max()', W.min(), W.max())
print('W.max() FE reference, simply supported', 6.594931610258557e-05)

W.min(), W.max() -4.526509226015919e-20 6.411304124327303e-05
W.max() FE reference, simply supported 6.594931610258557e-05

5Plotting in 3D using plotly¶

import plotly.graph_objects as go

fig = go.Figure(data=[go.Surface(x=X, y=Y, z=W)])

fig.update_layout(
    title='Plate Deflection (CLPT)',
    scene=dict(
        xaxis_title='X (m)',
        yaxis_title='Y (m)',
        zaxis_title='W (m)',
        aspectmode='data',
        camera=dict(
            eye=dict(x=0, y=0, z=100*a/2),
            center=dict(x=a/2, y=b/2, z=0),
            up=dict(x=0, y=1, z=0)
        )
    ),
    width=800,
    height=600
)

fig.show()