Post-buckling of plates

While historically buckling was seen as failure, modern engineering recognizes the postbuckling reserve of stiffened panels. These panels can withstand loads significantly exceeding their initial buckling threshold by allowing local buckling of the skin while stiffeners maintain global integrity.

1Differential quadrature¶

Methodological Comparison: DQM, iDQM, and MiDQM

When solving high-order boundary value problems such as the coupled Föppl–von Kármán (FvK) equations, the choice of spatial discretization dictates the numerical stability of the solver. We can categorize the Differential Quadrature (DQ) family into three distinct approaches.

1.1Pure direct Differential Quadrature Method (DQM)¶

TODO highlight the strong form

Pure DQM directly approximates the derivatives of a function at a set of grid points using a weighted linear sum of the function values at all other points. Higher-order derivative matrices are computed by simple matrix multiplication of the first-order weighting matrix ($D^{(m)} = D^{(1)} \dots D^{(1)}$).

Characteristics: Extremely simple to implement.

Drawback: As noted by Raju et al. (2013) Raju et al., 2013, Pure DQM is notorious for numerical instability in highly nonlinear post-buckling problems. The condition number of $D^{(m)}$ scales exponentially with $m$ and the number of grid points $N$.

1.2Pure inverse Differential Quadrature Method (iDQM)¶

Proposed by Ojo et al. (2021) Ojo et al., 2021, the Pure iDQM inverts the mathematical logic of DQM. Instead of approximating the function and differentiating, it directly approximates the highest-order derivative (e.g., $W_{,xxxx}$) using standard DQ. The lower-order derivatives and the function itself are recovered via integration weighting matrices ($H^{(m)}$), utilizing boundary conditions as constants of integration.

Characteristics: Requires the formulation of inverse boundary value problems to define the integration constants analytically.

Advantage: Integration is inherently a smoothing operation, meaning the condition number of the resulting algebraic system is significantly reduced.

Let’s critically evaluate the motivation of Ojo et al. (2021) Ojo et al., 2021. Ojo et al. proposed the iDQM fundamentally to bypass the extreme ill-conditioning and round-off error amplification caused by high-order differentiation matrices ($D^{(4)}$) in Pure DQM. When the condition number hits $\mathcal{O}(10^8)$, standard numerical solvers fail to converge. While Ojo’s mathematical diagnosis of the matrix condition number is correct, their conclusion that one must switch to integral matrices to achieve stability in high-order PDEs is a misdiagnosis of where the failure actually occurs in nonlinear structural solvers.

Take an analytical Airy stress-based DQM, such as the one implemented in this practice, where Ojo’s motivation is fundamentally bypassed. The instability in standard Pure DQM for FvK equations is rarely caused by a static matrix-vector multiplication $D^{(4)}W$. The failure actually occurs inside the Newton-Raphson solver during the finite-difference approximation of the Jacobian. When a standard optimizer perturbs the state vector by $\epsilon = 10^{-7}$, the $\mathcal{O}(10^8)$ condition number of $D^{(4)}$ amplifies this perturbation into massive $\mathcal{O}(0.1)$ numeric noise, destroying the descent direction. By deriving the exact analytical Jacobian using Kronecker tensor products, rhe finite-difference noise is entirely eliminated. Consequently, the Pure DQM achieves machine-precision convergence instantly, rendering the heavy machinery of integral $H$-matrices (iDQM) unnecessary.

When using the Ritz-DQM, Ojo’s motivation also becomes invalid. Ojo’s premise relies on the instability of $4^{\text{th}}$-order derivatives in the strong form. The Ritz method, however, operates on the Total Potential Energy ($\Pi$), which is the weak form. Furthermore, a lower derivative order is achieved with integration by parts when compared to the strong from, reducing the highest spatial derivative to $m=2$ (the bending curvatures $\kappa_{x,y}$), fundamentally reducing the numerical instability. Finally, the Ritz-DQM does not use dense, global differentiation matrices $D^{(m)}$. Instead, in the practice herein proposed, an analytical, hierarchical Legendre basis functions are proposed, which allow exact analytical derivatives to be calculated at Gauss-Legendre integration points. Because Ritz-DQM relies on exact polynomial derivatives and exact numerical quadrature, it entirely circumvents the round-off amplification issues that Ojo et al. aimed to solve with the inverse quadrature method.

1.3Mixed iDQM (MiDQM)¶

TODO find more about this

MiDQM is a hybrid approach discussed in contemporary literature (including extensions of Ojo et al.). It uses direct differentiation matrices ($D^{(m)}$) for lower-order derivative terms and integration matrices ($H^{(m)}$) for the highest-order terms.

Characteristics: Balances the smoothing benefits of integral operators with the straightforward boundary condition enforcement of standard differentiation.

2Differential quadrature (OLD)¶

The Differential Quadrature Method (DQM) is a global numerical technique used to solve partial differential equations (PDEs). Unlike the Finite Element Method (FEM) which relies on localized piecewise interpolation, DQM approximates the derivative of a function at a specific grid point as a weighted linear sum of the function values at all discrete points in the domain.

For a 1D function $f(x)$ discretized over $N$ points, the $n$-th order derivative at point $x_i$ is given by:

where $W_{ij}^{(n)}$ are the weighting coefficients. To suppress Runge’s phenomenon at higher polynomial orders, the grid points are not distributed uniformly, but rather follow the roots of orthogonal polynomials, such as Chebyshev-Gauss-Lobatto or Gauss-Legendre grids. While traditional DQM solves the strong form of the governing PDEs directly at the collocation points, this approach is mathematically brittle for 2D plate boundaries (specifically corners) due to linearly dependent constraints.

2.1The Ritz-DQ Method¶

To bypass corner singularities and maintain the exponential convergence of spectral methods, it’s possible to transition to the weak form via the Ritz-DQ method, which hybridises the classical Ritz variational method with DQM’s numerical integration rules. Instead of directly evaluating the PDEs, the Ritz-DQ method minimizes the Total Potential Energy ($\Pi$) of the system. The continuous displacement fields $(u, v, w)$ are expanded using a spectral basis (Legendre polynomials). The DQM machinery is then leveraged to evaluate the continuous energy integrals using exact Gauss-Legendre quadrature.

2.1.1Approximation functions¶

The transverse deflection $w(x,y)$ and in-plane displacements $u(x,y)$, $v(x,y)$ are approximated by expanding unknown Ritz coefficients $\mathbf{C}$ over 1D Legendre polynomials. To explicitly satisfy the simply supported (SSSS) kinematic boundary conditions ($w=0$ at edges), the following modified Legendre basis can be constructued:

Because standard Legendre polynomials $P_m(\pm 1) = (\pm 1)^m$, the difference exactly vanishes at the domain edges. The transverse field is thus:

2.1.2Strain energy¶

To capture post-buckling behavior, the strain-displacement relations must account for geometric nonlinearity. The von Karman mid-plane strains isolate the dominant transverse stiffening effects:

The Total Potential Energy ($\Pi$) is the sum of the membrane ($U_m$) and bending ($U_b$) strain energies. Defining extensional stiffness $C_{ext} = \frac{Eh}{1-\nu^2}$ and bending stiffness $D = \frac{Eh^3}{12(1-\nu^2)}$:

In the Ritz-DQM, the continuous integral $\int_A (\dots) dA$ is replaced by the 2D Gauss-Legendre quadrature summation $\sum \sum (\dots) W_{ij}$.

An example of the Ritz-DQ method can be seen in this notebook.

2.2Practice, Ritz-DQ Methodomputational Workflow¶

In the Ritz-DQ practice, the following algorithm is proposed:

1. Initialization: Generate the Gauss-Legendre quadrature roots $(\xi_i, \eta_j)$ and weights $W_{ij}$ for polynomial degree $N_c$.

2. Precomputation: Evaluate the 1D basis polynomials and their spatial derivatives at the grid points. Form the tensor products.

3. Displacement Mapping: Define a loop over incremental edge shortenings $\Delta u$.

4. Energy Minimization: At each load step, project the current coefficient guess $\mathbf{C}$ through the precomputed tensors to obtain physical strains. Calculate $\Pi$ using DQM quadrature. Use an optimizer (e.g., BFGS or SLSQP) to find the coefficient array that minimizes $\Pi$.

5. Symmetry Breaking: If the plate is unbuckled ($w \approx 0$), inject an artificial perturbation (e.g., $C^w_{00} = 10^{-3}$) prior to minimization to push the gradient off the unstable saddle point.

6. Force Recovery: Post-multiply the converged displacement gradients to extract the boundary membrane stress $N_{xx}$, integrating it via 1D quadrature to output the applied macroscopic load $P$.

2.2.1Boundary Condition Variations: Yamaki I vs. Yamaki III¶

The post-buckling stiffness of a plate is highly sensitive to the in-plane boundary conditions along the unloaded edges ($y = \pm b/2$). We analyze two foundational cases defined by Yamaki Yamaki, 1959.

2.2.1.1Yamaki Condition III: Stress-Free (Warping) Edges¶

• Physics: The unloaded edges are completely unrestrained in the $y$-direction. As the plate buckles and deflects in $z$, the edges are physically pulled inward. Because they are unrestrained, they warp and wave, locally relieving membrane tension.

• Implementation: This is a natural boundary condition. The optimizer minimizes the energy unconstrained, inherently finding the stress-free warped state. This results in lower geometric stiffness.

2.2.1.2Yamaki Condition I: Straight Edges¶

• Physics: The unloaded edges are supported by stiffeners that allow them to slide inward macroscopically (zero average stress, $\int N_{yy} dx = 0$), but force them to remain perfectly straight. This forces the entire edge to displace by the same amount, generating massive transverse membrane tension in the center of the plate.

• Implementation: This requires an explicit kinematic equality constraint. The variance of the $v$-displacement along the edges must be zero:

  $$v(x, \pm b/2) - \bar{v}_{\pm b/2} = 0$$

In this example this is enforced either via a severe energy penalty factor.

3Effective width¶

The postbuckling behavior and stress/strain distribution of stiffened panels is complex and non-linear. Complicated non-linear numerical calculation methods that employ sig- nificant computational resources are laborious and are required to confidently predict the panels ultimate load capacity Pevzner et al., 2008. To alleviate the calculations, a relatively simplified model, the so called “effective width” approach, has been proposed by von K'arman et al. Kármán et al., 1932 and subsequently modified by Cox Cox, 1933 and Sechler Sechler, 1937. This approach has provided a good average approximation for calculation of the effective width, $w_e$, i.e. the portion of the between adjacent stringers buckled skin, that together with the stringer constitute the integral skin-stringer combi- nation that participates in load carrying in postbuckling. The method works adequately for the case of uniaxial compression, and it is not recommended when there is biaxial loading or compression combined with shear Kassapoglou, 2013. Based on the average stress $s_{st}$ experienced by the stringers and the first critical skin stress, $s_{cr}$ between adjacent stringers of spacing b, the following relation has been pro- posed by Marguerre for determination of $w_e$:

The above effective width concept is widely and effectively applied as an adequate reliable tool for prediction of ultimate loads of metal flat stiffened panels. When appropriately modified and adapted it might lend itself as an appropriate approach for determination of ultimate load capacities of axially compressed laminated composite stringer-stiffened curved panels as well Pevzner et al., 2008.

The effective width method simplifies the complex and non-uniform stress distribution in a buckled panel, replacing it with an equivalent and uniform stress acting over a reduced “effective width” of the skin adjacent to the stiffeners.

3.1Effective width for metallic structures¶

3.2Effective width for composite plates¶

An example on how the effective width changes with the loading fraction and material properties can be seen in this notebook. An illustration on how the internal load changes over the skin width can be found in this other notebook

These are also available in this web version of the documentation, see: Kassapoglou, Fig. 7.10 and Kassapoglou, Fig. 7.12.

3.3Effective width for composite shells¶

The Technion Effective Width (TEW) method Pevzner et al., 2008 is an engineering approximation for analyzing the postbuckling behavior of curved, laminated composite structures. The TEW method extends the effective width concept to curved, anisotropic, laminated composite panels by reformulating an equivalent column model to account for the unique bending, torsional, and coupled instability modes of composite structures.

The TEW analysis process is summarized as follows:

1. First Buckling Calculation: Determining the initial local buckling of the skin between stringers using semi-empirical or approximate analytical solutions.

2. Iterative Convergence of Effective Width: Once the load exceeds the initial buckling load, an iterative algorithm calculates the effective width of the skin contributing to the load-carrying capacity. This process continues until the stress redistribution between the buckled skin and the stiffener reaches equilibrium.

3. Global Stability Analysis: Evaluating the global column stability based on the flexural, torsional, and warping rigidities of the equivalent skin-stringer cross-section to determine the ultimate collapse load.

An example of the TEW method is presented in this notebook.

Which is part of this web version of the documentation, see: TEW method.

The TEW paper reveals that the predictive fidelity of the TEW method is closely linked to the panel’s stiffness.

• Heavily Stiffened Panels (Cases I, V): The TEW method tends to overpredict the collapse load. This is likely because the model does not capture localized failure modes like skin-stiffener debonding that can occur before global buckling.

• Lightly Stiffened Panels (Cases II, IV): The TEW method tends to conservatively underpredict the collapse load. The small margin between initial skin buckling and global collapse in these flexible structures may lead the algorithm to predict failure prematurely.

• Asymmetric Stiffeners (Cases IV, V): The use of J-form stringers introduces coupled bending-torsion modes, which significantly complicates the buckling behavior and increases sensitivity to manufacturing imperfections.