Practice, deflection of a plate using the FSDT¶

!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

# approximation order
m1 = m2 = 20
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 (simply supported)
wxit1 = 0
wxir1 = 1
wxit2 = 0
wxir2 = 1
wetat1 = 0
wetar1 = 1
wetat2 = 0
wetar2 = 1

# Boundary conditions phix (simply supported: constrained at y=0, b)
phix_xit1 = 1
phix_xir1 = 1
phix_xit2 = 1
phix_xir2 = 1
phix_etat1 = 0
phix_etar1 = 1  # Derivative producing twist, must be free
phix_etat2 = 0
phix_etar2 = 1  # Derivative producing twist, must be free

# Boundary conditions phiy (simply supported: constrained at x=0, a)
phiy_xit1 = 0
phiy_xir1 = 1  # Derivative producing twist, must be free
phiy_xit2 = 0
phiy_xir2 = 1  # Derivative producing twist, must be free
phiy_etat1 = 1
phiy_etar1 = 1
phiy_etat2 = 1
phiy_etar2 = 1

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

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

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)


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

detJ = a * b / 4

# 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)
    
    phix_P_xi = vecf(m1, xi, phix_xit1, phix_xir1, phix_xit2, phix_xir2)
    phix_Px_xi = vecfxi(m1, xi, phix_xit1, phix_xir1, phix_xit2, phix_xir2)
    
    phiy_P_xi = vecf(m1, xi, phiy_xit1, phiy_xir1, phiy_xit2, phiy_xir2)
    phiy_Px_xi = vecfxi(m1, xi, phiy_xit1, phiy_xir1, phiy_xit2, phiy_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)

        phix_P_eta = vecf(m2, eta, phix_etat1, phix_etar1, phix_etat2, phix_etar2)
        phix_Px_eta = vecfxi(m2, eta, phix_etat1, phix_etar1, phix_etat2, phix_etar2)

        phiy_P_eta = vecf(m2, eta, phiy_etat1, phiy_etar1, phiy_etat2, phiy_etar2)
        phiy_Px_eta = vecfxi(m2, eta, phiy_etat1, phiy_etar1, phiy_etat2, phiy_etar2)

        Sphix[m1*m2:2*m1*m2]  = (phix_P_xi[:, None] * phix_P_eta[None, :]).ravel()
        Sphiy[2*m1*m2:]       = (phiy_P_xi[:, None] * phiy_P_eta[None, :]).ravel()
        
        Swx[:m1*m2]           = (wPx_xi[:, None] * wP_eta[None, :] * (2/a)).ravel()
        Sphixx[m1*m2:2*m1*m2] = (phix_Px_xi[:, None] * phix_P_eta[None, :] * (2/a)).ravel()
        Sphiyx[2*m1*m2:]      = (phiy_Px_xi[:, None] * phiy_P_eta[None, :] * (2/a)).ravel()
        
        Swy[:m1*m2]           = (wP_xi[:, None] * wPx_eta[None, :] * (2/b)).ravel()
        Sphixy[m1*m2:2*m1*m2] = (phix_P_xi[:, None] * phix_Px_eta[None, :] * (2/b)).ravel()
        Sphiyy[2*m1*m2:]      = (phiy_P_xi[:, None] * phiy_Px_eta[None, :] * (2/b)).ravel()
        
        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 = 5/6 * (prop.A44*g0yz + prop.A45*g0xz)
        Qx = 5/6 * (prop.A45*g0yz + prop.A55*g0xz)
        
        weight = wxi * weta

        # stiffness matrix
        K += np.outer(detJ*weight*Mxx, e1xx)
        K += np.outer(detJ*weight*Myy, e1yy)
        K += np.outer(detJ*weight*Mxy, e1xy)
        K += np.outer(detJ*weight*Qy, g0yz)
        K += np.outer(detJ*weight*Qx, g0xz)

m1 = 20
m2 = 20
N = 1200

1External force vector¶

Sw = np.zeros(N)

xi = 0 # x = a/2
eta = 0 # y = b/2

wP_xi = vecf(m1, xi, wxit1, wxir1, wxit2, wxir2)
wP_eta = vecf(m2, eta,  wetat1, wetar1, wetat2, wetar2)
Sw[:m1*m2] = (wP_xi[:, None] * wP_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...
				156 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)

Sw = np.zeros(N)

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

    wP_xi = vecf(m1, xi, wxit1, wxir1, wxit2, wxir2)
    wP_eta = vecf(m2, eta, wetat1, wetar1, wetat2, wetar2)

    Sw[:m1*m2] = (wP_xi[:, None] * wP_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.488089304725298e-20 6.411427399380127e-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 (FSDT)',
    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()