Practice, discrete Ritz-DQ Method for Post-Buckling of Plates¶

This notebook implements the discrete Ritz method using Differential Quadrature (Ritz-DQ) to solve the nonlinear post-buckling behavior of an isotropic plate under uniaxial compression. It compares Yamaki Condition I (straight unloaded edges) and Yamaki Condition III (stress-free warping edges).

1Von Karman Nonlinear Kinematics¶

To capture post-buckling, we must account for the geometric stiffening caused by transverse deflections stretching the midplane. The von Karman nonlinear strain-displacement relations are:

\varepsilon_{xx}=u_{,x}+\frac{1}{2}w_{,x}^2

\varepsilon_{yy}=v_{,y}+\frac{1}{2}w_{,y}^2

\gamma_{xy}=u_{,y}+v_{,x}+w_{,x}w_{,y}

2Total Potential Energy ( $\Pi$ )¶

The total potential energy is the sum of the membrane strain energy ( $U_m$ ) and the bending strain energy ( $U_b$ ). For an isotropic plate with extensional stiffness $C_{ext}=\frac{Eh}{1-\nu^2}$ and bending stiffness $D=\frac{Eh^3}{12(1-\nu^2)}$ :

Membrane Energy:

U_m=\frac{1}{2}\int_A(N_{xx}\varepsilon_{xx}+N_{yy}\varepsilon_{yy}+N_{xy}\gamma_{xy})dA

Substituting the constitutive law ( $N=C_{ext}\varepsilon$ ):

U_m=\frac{C_{ext}}{2}\int_A(\varepsilon_{xx}^2+\varepsilon_{yy}^2+2\nu\varepsilon_{xx}\varepsilon_{yy}+\frac{1-\nu}{2}\gamma_{xy}^2)dA

Bending Energy:

U_b=\frac{1}{2}\int_A(M_{xx}\kappa_{xx}+M_{yy}\kappa_{yy}+M_{xy}\kappa_{xy})dA

Substituting the curvatures ( $\kappa_{xx}=-w_{,xx}$ ):

U_b=\frac{D}{2}\int_A(w_{,xx}^2+w_{,yy}^2+2\nu w_{,xx}w_{,yy}+2(1-\nu)w_{,xy}^2)dA

3The Discrete Ritz-DQ Discretization¶

Instead of solving the strong-form PDEs (which leads to severe boundary condition singularities), we minimize the energy $\Pi$ . In the Discrete Ritz method, the continuous integrals are replaced by numerical quadrature. We use a Gauss-Legendre grid $(\xi_i, \eta_j)$ with weights $W_{ij}$ :

\Pi\approx\sum_{i=1}^{N_c}\sum_{j=1}^{N_c}(U_m(\xi_i,\eta_j)+U_b(\xi_i,\eta_j))W_{ij}

4Penalty Method for Yamaki Condition I¶

Condition III requires stress-free unloaded edges, which is the natural boundary condition of the energy functional (no constraints needed). Condition I requires the unloaded edges ( $y=\pm b/2$ ) to remain straight but free to slide. This is a kinematic constraint: $v(x, \pm b/2)=constant$ .

To enforce this without altering the spectral basis, we add a penalty term to the energy functional that minimizes the statistical variance of the displacement along those edges:

\Pi_{penalty}=\lambda\left(\text{Var}(v|_{y=b/2})+\text{Var}(v|_{y=-b/2})\right)

Note on Penalty Stiffness: As you observed, the penalty factor $\lambda=10^7$ acts as an upper limit. If $\lambda$ is too small, the edges warp. If $\lambda$ is too large (e.g., $\geq 10^{9}$ ), it dominates the true physical energy gradients, creating a highly ill-conditioned Hessian matrix that causes the BFGS optimizer to oscillate or fail to converge.

import numpy as np
from numpy.polynomial.legendre import legval, legder
from scipy.optimize import minimize
from scipy.special import roots_legendre
import plotly.graph_objects as go

# ==========================================
# 1. PLATE PROPERTIES
# ==========================================
E, nu, h, a, b = 70e9, 0.3, 0.01, 1.0, 1.0
D = (E * h**3) / (12 * (1 - nu**2))
C_ext = E * h / (1 - nu**2)

# Jacobian mapping from computational [-1, 1] to physical [-a/2, a/2] domain
Jx, Jy = 2.0 / a, 2.0 / b

5Gauss-Legendre Quadrature Grid¶

We use the roots of the Legendre polynomial to define the grid. Gauss-Legendre quadrature provides exact integration for polynomials up to degree $2N_c-1$ , making it highly efficient for energy methods.

# ==========================================
# 2. RITZ-DQ GRID (GAUSS-LEGENDRE)
# ==========================================
N_c = 8  # Polynomial approximation order
xi, w_gl = roots_legendre(N_c)
eta = xi.copy()
XI, ETA = np.meshgrid(xi, eta)
xi_f, eta_f = XI.flatten(), ETA.flatten()

# 2D Quadrature Weights for Energy Integration (scaled by physical domain area)
W_2D = np.outer(w_gl, w_gl).flatten() * (a / 2) * (b / 2)

6Spectral Basis Functions¶

We approximate the displacement fields using a series expansion of continuous polynomials:

w(\xi,\eta)=\sum C^w_{ij}\phi_i(\xi)\phi_j(\eta)

To rapidly satisfy kinematic boundary conditions, we construct a modified_legendre basis:

\phi_m(\xi)=P_{m+2}(\xi)-P_m(\xi)

Because standard Legendre polynomials evaluate to $\pm 1$ at the boundaries, this subtraction ensures $\phi_m(\pm 1)=0$ exactly. We apply this to $w$ on all edges, and $u$ on the loaded edges.

# ==========================================
# 3. SPECTRAL BASIS
# ==========================================
def modified_legendre(pts, m, deriv=0):
    """ Enforces exactly 0 at boundaries: P_{m+2} - P_m """
    c_m2, c_m = np.zeros(m + 3), np.zeros(m + 1)
    c_m2[m + 2], c_m[m] = 1.0, 1.0
    if deriv > 0:
        c_m2, c_m = legder(c_m2, m=deriv), legder(c_m, m=deriv)
    v1 = legval(pts, c_m2) if len(c_m2) > 0 else np.zeros_like(pts)
    v2 = legval(pts, c_m) if len(c_m) > 0 else np.zeros_like(pts)
    return v1 - v2

def standard_legendre(pts, m, deriv=0):
    """ Standard P_m for unconstrained fields """
    c = np.zeros(m + 1); c[m] = 1.0
    if deriv > 0: c = legder(c, m=deriv)
    return legval(pts, c) if len(c) > 0 else np.zeros_like(pts)

# U: Modified in X (enforces u=0 to cleanly apply prescribed delta_u), Standard in Y
P_ux = np.array([[modified_legendre(xi_f, i, d) for i in range(N_c)] for d in range(3)])
P_uy = np.array([[standard_legendre(eta_f, j, d) for j in range(N_c)] for d in range(3)])

# V: Standard in X, Standard in Y (Free Poisson expansion)
P_vx = np.array([[standard_legendre(xi_f, i, d) for i in range(N_c)] for d in range(3)])
P_vy = np.array([[standard_legendre(eta_f, j, d) for j in range(N_c)] for d in range(3)])

# W: Modified in X, Modified in Y (SSSS boundaries w=0)
P_wx = np.array([[modified_legendre(xi_f, i, d) for i in range(N_c)] for d in range(3)])
P_wy = np.array([[modified_legendre(eta_f, j, d) for j in range(N_c)] for d in range(3)])

def get_field(C, bx, by, dx, dy):
    """ Contructs the full field by contracting the coefficient tensor with the 1D bases """
    return np.einsum('ij, ik, jk -> k', C, bx[dx], by[dy]) * (Jx**dx) * (Jy**dy)

# 1D Basis arrays used exclusively to extract the V displacement along the unloaded edges for the penalty constraint
P_vx_1d = np.array([standard_legendre(xi, i, 0) for i in range(N_c)])
P_vy_pos = np.array([standard_legendre(np.array([1.0]), j, 0)[0] for j in range(N_c)])
P_vy_neg = np.array([standard_legendre(np.array([-1.0]), j, 0)[0] for j in range(N_c)])

7Total Potential Energy Formulation¶

This is the core objective function. The displacement fields and their derivatives are assembled, the non-linear strains are calculated, and the energy is summed over the quadrature grid.

The energy is normalized by $D$ to scale the objective value down. This drastically improves the condition number of the numerical Hessian matrix computed by the BFGS algorithm.

# ==========================================
# 4. TOTAL POTENTIAL ENERGY FUNCTIONAL
# ==========================================
def energy(C_flat, delta_u, condition):
    Cu, Cv, Cw = C_flat.reshape(3, N_c, N_c)
    
    # 4.1 Kinematics (Note: rigid delta_u is linearly superimposed on the u-field)
    u_x = get_field(Cu, P_ux, P_uy, 1, 0) - (delta_u / a)
    u_y = get_field(Cu, P_ux, P_uy, 0, 1)
    v_x = get_field(Cv, P_vx, P_vy, 1, 0)
    v_y = get_field(Cv, P_vx, P_vy, 0, 1)
    w_x = get_field(Cw, P_wx, P_wy, 1, 0)
    w_y = get_field(Cw, P_wx, P_wy, 0, 1)
    w_xx = get_field(Cw, P_wx, P_wy, 2, 0)
    w_yy = get_field(Cw, P_wx, P_wy, 0, 2)
    w_xy = get_field(Cw, P_wx, P_wy, 1, 1)
    
    # 4.2 Von Karman Strains
    exx = u_x + 0.5 * w_x**2
    eyy = v_y + 0.5 * w_y**2
    gxy = u_y + v_x + w_x * w_y
    
    # 4.3 Strain Energy Densities
    U_m = 0.5 * C_ext * (exx**2 + eyy**2 + 2*nu*exx*eyy + 0.5*(1-nu)*gxy**2)
    U_b = 0.5 * D * (w_xx**2 + w_yy**2 + 2*nu*w_xx*w_yy + 2*(1-nu)*w_xy**2)
    
    # 4.4 DQM Integration & Normalization
    Pi = np.sum((U_m + U_b) * W_2D) / D
    
    # 4.5 Rigid Body Translation Penalty in V
    v_field = get_field(Cv, P_vx, P_vy, 0, 0)
    Pi += 0.5 * 1e4 * np.mean(v_field)**2
    
    # 4.6 Yamaki Condition I: Straight Edge Penalty
    if condition == 'I':
        v_edge_pos = np.einsum('ij, ik, j -> k', Cv, P_vx_1d, P_vy_pos)
        v_edge_neg = np.einsum('ij, ik, j -> k', Cv, P_vx_1d, P_vy_neg)
        # The variance penalty forces the edge to remain straight while allowing the mean to translate
        Pi += 1e7 * (np.var(v_edge_pos) + np.var(v_edge_neg))
        
    return Pi

8Post-Buckling Tracking Loop¶

The solver steps through displacement increments. At each step, a minor initial imperfection is injected if the plate is flat; this symmetry breaker prevents the BFGS algorithm from getting trapped at the $w=0$ saddle point during the bifurcation.

After convergence, the reaction force is evaluated by integrating the boundary stress $N_{xx}$ specifically along $x=-a/2$ .

# ==========================================
# 5. BOUNDARY PRECOMPUTATIONS & SOLVER
# ==========================================
# Precompute basis specifically evaluated at the left edge (xi = -1) to extract reaction forces
xi_e = np.array([-1.0])
XI_e, ETA_e = np.meshgrid(xi_e, eta)
xie_f, etae_f = XI_e.flatten(), ETA_e.flatten()

P_ux_e = np.array([[modified_legendre(xie_f, i, d) for i in range(N_c)] for d in range(2)])
P_uy_e = np.array([[standard_legendre(etae_f, j, d) for j in range(N_c)] for d in range(2)])
P_vx_e = np.array([[standard_legendre(xie_f, i, d) for i in range(N_c)] for d in range(2)])
P_vy_e = np.array([[standard_legendre(etae_f, j, d) for j in range(N_c)] for d in range(2)])
P_wx_e = np.array([[modified_legendre(xie_f, i, d) for i in range(N_c)] for d in range(2)])
P_wy_e = np.array([[modified_legendre(etae_f, j, d) for j in range(N_c)] for d in range(2)])

P_cr_theoretical = 4 * np.pi**2 * D / (b**2)

def run_postbuckling(condition, steps=20, max_du=1.5e-3):
    print(f"\nRunning Yamaki Condition {condition}...")
    du_array = np.linspace(0, max_du, steps)
    P_results, w_results = [], []
    C_opt = np.zeros(3 * N_c**2)

    for i, du in enumerate(du_array):
        
        # SYMMETRY BREAKER: Re-seed fundamental mode if unbuckled to escape the saddle point
        if np.max(np.abs(C_opt[2*N_c**2:])) < 1e-4:
            C_opt[2*N_c**2] = 1e-3 
            
        res = minimize(energy, C_opt, args=(du, condition), method='BFGS', options={'disp': False, 'gtol': 1e-6})
        
        if res.success or res.status == 2:
            C_opt = res.x
            Cu, Cv, Cw = C_opt.reshape(3, N_c, N_c)
            
            w_field = get_field(Cw, P_wx, P_wy, 0, 0)
            w_results.append(np.max(np.abs(w_field)) / h)
            
            # Extract strains precisely at the boundary to integrate stress
            u_x_e = get_field(Cu, P_ux_e, P_uy_e, 1, 0) - (du / a)
            v_y_e = get_field(Cv, P_vx_e, P_vy_e, 0, 1)
            w_x_e = get_field(Cw, P_wx_e, P_wy_e, 1, 0)
            w_y_e = get_field(Cw, P_wx_e, P_wy_e, 0, 1)
            
            Nxx_e = C_ext * ((u_x_e + 0.5*w_x_e**2) + nu*(v_y_e + 0.5*w_y_e**2))
            P_avg = np.sum(Nxx_e * w_gl) * 0.5
            P_results.append(np.abs(P_avg) / P_cr_theoretical)
            
            print(f"Step {i+1:02d} | du: {du:.2e} | P/P_cr: {P_results[-1]:.3f} | W/h: {w_results[-1]:.3f}")
        else:
            print(f"Failed at du = {du}")
            break
            
    return P_results, w_results, C_opt

# Execute both boundary configurations
P_I, w_I, C_opt_I = run_postbuckling('I')
P_III, w_III, C_opt_III = run_postbuckling('III')


Running Yamaki Condition I...
Step 01 | du: 0.00e+00 | P/P_cr: 0.000 | W/h: 0.000
Step 02 | du: 7.89e-05 | P/P_cr: 0.218 | W/h: 0.000
Step 03 | du: 1.58e-04 | P/P_cr: 0.436 | W/h: 0.000
Step 04 | du: 2.37e-04 | P/P_cr: 0.655 | W/h: 0.000
Step 05 | du: 3.16e-04 | P/P_cr: 0.873 | W/h: 0.000
Step 06 | du: 3.95e-04 | P/P_cr: 1.045 | W/h: 0.338
Step 07 | du: 4.74e-04 | P/P_cr: 1.152 | W/h: 0.619
Step 08 | du: 5.53e-04 | P/P_cr: 1.257 | W/h: 0.807
Step 09 | du: 6.32e-04 | P/P_cr: 1.359 | W/h: 0.957
Step 10 | du: 7.11e-04 | P/P_cr: 1.460 | W/h: 1.086
Step 11 | du: 7.89e-04 | P/P_cr: 1.558 | W/h: 1.200
Step 12 | du: 8.68e-04 | P/P_cr: 1.654 | W/h: 1.304
Step 13 | du: 9.47e-04 | P/P_cr: 1.749 | W/h: 1.399
Step 14 | du: 1.03e-03 | P/P_cr: 1.841 | W/h: 1.487
Step 15 | du: 1.11e-03 | P/P_cr: 1.932 | W/h: 1.569
Step 16 | du: 1.18e-03 | P/P_cr: 2.020 | W/h: 1.646
Step 17 | du: 1.26e-03 | P/P_cr: 2.107 | W/h: 1.719
Step 18 | du: 1.34e-03 | P/P_cr: 2.192 | W/h: 1.788
Step 19 | du: 1.42e-03 | P/P_cr: 2.276 | W/h: 1.854
Step 20 | du: 1.50e-03 | P/P_cr: 2.358 | W/h: 1.917

Running Yamaki Condition III...
Step 01 | du: 0.00e+00 | P/P_cr: 0.000 | W/h: 0.000
Step 02 | du: 7.89e-05 | P/P_cr: 0.218 | W/h: 0.000
Step 03 | du: 1.58e-04 | P/P_cr: 0.436 | W/h: 0.000
Step 04 | du: 2.37e-04 | P/P_cr: 0.655 | W/h: 0.000
Step 05 | du: 3.16e-04 | P/P_cr: 0.873 | W/h: 0.000
Step 06 | du: 3.95e-04 | P/P_cr: 1.037 | W/h: 0.367
Step 07 | du: 4.74e-04 | P/P_cr: 1.122 | W/h: 0.674
Step 08 | du: 5.53e-04 | P/P_cr: 1.204 | W/h: 0.878
Step 09 | du: 6.32e-04 | P/P_cr: 1.282 | W/h: 1.043
Step 10 | du: 7.11e-04 | P/P_cr: 1.357 | W/h: 1.185
Step 11 | du: 7.89e-04 | P/P_cr: 1.429 | W/h: 1.310
Step 12 | du: 8.68e-04 | P/P_cr: 1.498 | W/h: 1.424
Step 13 | du: 9.47e-04 | P/P_cr: 1.565 | W/h: 1.529
Step 14 | du: 1.03e-03 | P/P_cr: 1.630 | W/h: 1.627
Step 15 | du: 1.11e-03 | P/P_cr: 1.694 | W/h: 1.719
Step 16 | du: 1.18e-03 | P/P_cr: 1.755 | W/h: 1.806
Step 17 | du: 1.26e-03 | P/P_cr: 1.815 | W/h: 1.888
Step 18 | du: 1.34e-03 | P/P_cr: 1.873 | W/h: 1.967
Step 19 | du: 1.42e-03 | P/P_cr: 1.930 | W/h: 2.042
Step 20 | du: 1.50e-03 | P/P_cr: 1.986 | W/h: 2.114

9Verification and Plotting¶

The results are plotted against the exact analytical benchmark for Yamaki Condition I (Levy, 1942).

# ==========================================
# 6. PLOT
# ==========================================
fig = go.Figure()
fig.add_trace(go.Scatter(x=P_I, y=w_I, mode='lines+markers', name='Ritz Model (Yamaki I: Straight)'))
fig.add_trace(go.Scatter(x=P_III, y=w_III, mode='lines+markers', name='Ritz Model (Yamaki III: Warping)'))

levy_P = [1.0, 1.15, 1.55, 2.15, 3.00]
levy_W = [0.0, 0.50, 1.00, 1.50, 2.00]
fig.add_trace(go.Scatter(x=levy_P, y=levy_W, mode='markers', marker=dict(symbol='x', size=10, color='red'), name='Levy/Yamaki I Exact Benchmark'))

fig.update_layout(
    title='Post-Buckling Kinematics: Straight vs. Warping Edges',
    xaxis_title='Load Ratio P / P_cr',
    yaxis_title='Normalized Deflection W_max / h',
    width=800
)
fig.show()

10Comparing solutions for straight (Yamaki I) with warped (Yamaki III) edges¶

from plotly.subplots import make_subplots

# ==========================================
# 6. IN-PLANE WARPING VISUALIZATION
# ==========================================
def plot_warping_comparison(C_flat_I, C_flat_III, scale_inplane=50):
    """
    Renders 3D surface plots with in-plane deformations applied to the physical grid.
    Yamaki I (Straight) vs Yamaki III (Warped).
    """
    res = 60
    xi_d = np.linspace(-1, 1, res)
    eta_d = np.linspace(-1, 1, res)
    XI_d, ETA_d = np.meshgrid(xi_d, eta_d, indexing='ij')
    
    X_phys = XI_d * (a / 2)
    Y_phys = ETA_d * (b / 2)
    
    # Dense basis evaluation (0th derivative)
    P_ux_d = np.array([modified_legendre(xi_d, i, 0) for i in range(N_c)])
    P_uy_d = np.array([standard_legendre(eta_d, j, 0) for j in range(N_c)])
    P_vx_d = np.array([standard_legendre(xi_d, i, 0) for i in range(N_c)])
    P_vy_d = np.array([standard_legendre(eta_d, j, 0) for j in range(N_c)])
    P_wx_d = np.array([modified_legendre(xi_d, i, 0) for i in range(N_c)])
    P_wy_d = np.array([modified_legendre(eta_d, j, 0) for j in range(N_c)])
    
    def get_deformed_grid(C_flat):
        Cu, Cv, Cw = C_flat.reshape(3, N_c, N_c)
        U_d = np.einsum('ij, ix, jy -> xy', Cu, P_ux_d, P_uy_d)
        V_d = np.einsum('ij, ix, jy -> xy', Cv, P_vx_d, P_vy_d)
        W_d = np.einsum('ij, ix, jy -> xy', Cw, P_wx_d, P_wy_d)
        
        # Apply amplified in-plane displacements to grid
        X_def = X_phys + scale_inplane * U_d
        Y_def = Y_phys + scale_inplane * V_d
        return X_def, Y_def, W_d

    X_I, Y_I, W_I = get_deformed_grid(C_flat_I)
    X_III, Y_III, W_III = get_deformed_grid(C_flat_III)
    
    fig = make_subplots(
        rows=1, cols=2, 
        specs=[[{'type': 'surface'}, {'type': 'surface'}]],
        subplot_titles=('Yamaki I: Straight Edges', 'Yamaki III: Warped Edges')
    )
    
    fig.add_trace(go.Surface(x=X_I, y=Y_I, z=W_I, colorscale='Viridis', showscale=False), row=1, col=1)
    fig.add_trace(go.Surface(x=X_III, y=Y_III, z=W_III, colorscale='Plasma', colorbar=dict(title='W (m)', x=1.05)), row=1, col=2)
    
    camera_top_down = dict(eye=dict(x=0.0, y=0.0, z=2.5))
    
    fig.update_layout(
        title=f"In-Plane Edge Warping Comparison (In-plane scaled {scale_inplane}x)",
        scene=dict(xaxis_title='X + U', yaxis_title='Y + V', zaxis_title='W', aspectratio=dict(x=1, y=b/a, z=0.3), camera=camera_top_down),
        scene2=dict(xaxis_title='X + U', yaxis_title='Y + V', zaxis_title='W', aspectratio=dict(x=1, y=b/a, z=0.3), camera=camera_top_down),
        width=1200, height=600
    )
    fig.show()

plot_warping_comparison(C_opt_I, C_opt_III, scale_inplane=50)

Practice, discrete Ritz-DQ Method for Post-Buckling of Plates