User:Eas4200c.f08/Plate buckling

Rectangular plate under in-plane distributed load

Simply-supported boundary conditions
Refs: MIT OCW 2.081J / 16.230J Plates and Shells Spring 2007 Lecture notes, pp.54-61. Sun, Mechanics of Aircraft Structures, Wiley, 2006, Section 7.7. Timoshenko & Gere, Theory of elastic stability, McGraw-Hill, 1961, Sec.9.2, p.351.

For a simply-supported rectangular plate under in-plane distributed compressive load, see figure Lecture notes, p.57. Plate dimensions: Length $$\displaystyle a$$ along $$\displaystyle x$$ axis, length $$\displaystyle b$$ along $$\displaystyle y$$ axis, thickness $$\displaystyle h$$. Young's modulus $$\displaystyle E$$, Poisson's ratio $$\displaystyle \nu$$. Total resultant compressive load $$\displaystyle P_x$$ (force) pointed in the $$\displaystyle x$$ direction, distributed along the $$\displaystyle y$$ axis.

Buckling mode shapes $$\displaystyle \psi_{mn} (x,y)$$
 * {| style="width:100%" border="0"

$$  \displaystyle \psi (x,y) =   c_{mn} \sin \left(     \frac      {m \pi x}      {a}   \right) \sin \left(     \frac      {n \pi y}      {b}   \right) \, \ {\rm for} \ m,n = 1, 2, 3, \ldots $$ where $$\displaystyle c_{mn}$$ is a constant, $$\displaystyle (m,n)$$ are the number of half wave-lengths along $$(x,y)\;$$, respectively. The transverse displacement (out of plane) $$\displaystyle u_z (x,y)$$ is a linear combination of the buckling mode shapes:
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (1a)
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle u_z (x,y) =  \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} \psi_{mn} (x,y) $$
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (1b)
 * }

Critical buckling load $$\displaystyle (P_x)_{cr}$$


 * {| style="width:100%" border="0"

$$  \displaystyle (P_x)_{cr} =  k_c \frac {\pi ^2 D}  {b} $$
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (2)
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle k_c (m, a/b) :=  \left(      \frac      {m b}      {a}      +      \frac      {a}      {m b}   \right)^2 $$
 * style="width:95%" |
 * style="width:95%" |
 * style= | (3)
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle D  := \frac {E h^3} {12 (1 - \nu^2)} $$ The factor $$\displaystyle D$$ is called the plate bending stiffness, and is the equivalent of the bending stiffness $$\displaystyle EI$$ in beam bending. Note that the dimension of the second area moment of inertia $$\displaystyle I$$, denoted by $$\displaystyle [I]$$, is $$\displaystyle L^4$$, i.e., length to the power 4:
 * style="width:95%" |
 * style="width:95%" |
 * style= | (4)
 * }



\displaystyle [I] = L^4 $$ Since $$\displaystyle \sigma = E \epsilon$$, we have

\displaystyle [\sigma] =  [E] [\epsilon] =  [E] =  \frac {F} {L^2} $$ i.e., force per length squared, as $$\displaystyle [\epsilon] = 1$$. It follows that the dimension of the beam bending stiffness coefficient

\displaystyle [EI] =   \frac {F} {L^2} L^4 =  FL^2 $$ whereas the dimension of plate bending stiffness coefficient $$\displaystyle D$$ is

\displaystyle [D] =  [E] [h^3] =  \frac {F} {L^2} L^3 =  FL $$ since the dimension of the Poisson's ratio is $$\displaystyle [\nu] = 1$$.

In addition, with $$\displaystyle [k_c] = 1$$, we have

\displaystyle [(P_x)_{cr}] =  [k_c] \frac {[D]} {[b]} =  1   \cdot \frac {FL} {L} =  F $$ which is consistent with the definition of $$\displaystyle P_x$$, the total in-plane compressive load.

Critical buckling stress $$\displaystyle (\sigma_{xx})_{cr}$$


 * {| style="width:100%" border="0"

$$  \displaystyle (\sigma_{xx})_{cr} =  \frac {(P_x)_{cr}} {bh} =  k_c \frac {\pi ^2 D}  {b^2 h} $$
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (7)
 * }

Since the factor $$\displaystyle k_c$$ depends on the number of half wave-lengths $$\displaystyle m$$ and the aspect ratio $$\displaystyle a/b$$, for a given aspect ratio $$\displaystyle a/b$$ of the rectangular plate, the minimum critical load $$\displaystyle (\sigma_{xx})_{cr}^\star$$ and the associated number of half wave-lengths $$\displaystyle m^\star$$ are obtained by minimizing the function $$\displaystyle k_c$$, i.e.,
 * {| style="width:100%" border="0"

$$  \displaystyle (\sigma_{xx})_{cr}^\star =  k_c^\star \frac {\pi ^2 D}  {b^2 h} $$ with
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (8)
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle k_c^\star (a/b) =  \min_{m = 1, 2, \cdots} \ k_c (m, a/b) =  \min_{m = 1, 2, \cdots} \left(     \frac      {m b}      {a}      +      \frac      {a}      {m b}   \right)^2 =  \left(      \frac      {m^\star b}      {a}      +      \frac      {a}      {m^\star b}   \right)^2 $$ where
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (9)
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle m^\star (a/b) =  \underset{m=1,2, \cdots}{\rm argmin} \ k_c (m, a/b) =  \underset{m=1,2, \cdots}{\rm argmin} \left(     \frac      {m b}      {a}      +      \frac      {a}      {m b}   \right)^2 $$ i.e., $$\displaystyle m^\star (a/b)$$ is the minimizer of $$\displaystyle k_c (m,a/b)$$ for a given aspect ratio $$\displaystyle a/b$$.
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (10)
 * }

For a simply-supported column, see Lecture notes, p.77 and Sun [2006, Sec 7.2], the critical load is


 * {| style="width:100%" border="0"

$$  \displaystyle (P_{cr})_{column} =  \frac {\pi^2 EI} {L^2} $$ where $$\displaystyle L$$ is the length of the column.
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (11)
 * }

See also Lecture notes, pp.4-5, University of Auckland, New Zealand.

Clamped boundary conditions
For rectangular plates with clamped boundary conditions on all 4 edges under in-plane distributed load, use the Excel file for calculating critical stress: Be sure to save the excel file you downloaded somewhere since you may not be able to download it again within one day.

The above excel file is based on the table in the book by Young & Budynas, Roark's formula for stress and strain, McGraw-Hill, 7th edition, 2002, p.730, where the following formula and notation were used
 * {| style="width:100%" border="0"

$$  \displaystyle \sigma^\prime =  K   \frac{E}{1 - \nu^2} \left(     \frac{t}{b}   \right)^2 $$ where $$\displaystyle \sigma^\prime \equiv (\sigma_{xx})_{cr}$$, $$\displaystyle K$$ is a function different from $$\displaystyle k_c$$, and $$\displaystyle t \equiv h$$ the plate thickness.
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (12)
 * }

For simply-supported plates, comparing Eq.(12) to Eq.(7), one can deduce that
 * {| style="width:100%" border="0"

$$  \displaystyle K   = k_c \frac{\pi^2}{12} =   \left(      \frac      {m b}      {a}      +      \frac      {a}      {m b}   \right)^2 \frac{\pi^2}{12} $$ For clamped plates, we still have $$\displaystyle K = k_c \pi^2 / 12$$, but the coefficient $$\displaystyle k_c$$ is no longer given by Eq.(3).
 * style="width:95%" |
 * style="width:95%" |
 * style= | (13)
 * }

When the clamped plate is close to being square, the following (buckling) shape $$\displaystyle \psi (x,y)$$ can be used to estimate the buckling load with reasonable accuracy (Timoshenko & Gere [1961], Theory of elastic stability, McGraw-Hill, p.386.)


 * {| style="width:100%" border="0"

$$  \displaystyle \psi (x,y) =  \frac {c} {4}  \left(      1 - \cos \frac{2 \pi x}{a}   \right) \left(     1 - \cos \frac{2 \pi y}{b}   \right) $$
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (14)
 * }

Another calculator of critical stress (free for square plate), Poisson's ratio cannot be changed, not useful.

Simply-supported boundary conditions
Refs: Timoshenko & Gere, Theory of elastic stability, McGraw-Hill, 1961, Sec.9.7, p.379. Sun, Mechanics of Aircraft Structures, Wiley, 2006, Section 7.7.4.

Because of the simply-supported conditions, the buckling mode shapes $$\displaystyle \psi (x,y)$$ are given as in Eq.(1a).

Denote the shear force per unit length of the plate as $$\displaystyle N_{xy}$$, i.e.,
 * {| style="width:100%" border="0"

$$  \displaystyle N_{xy} =  \sigma_{xy} h $$ where $$\displaystyle h$$ is the plate thickness. See Fig.7.23 on p.263 in Sun [2006]. In the skin or spar webs of an airfoil section, $$\displaystyle N_{xy}$$ is the same as the shear flow $$\displaystyle q$$.
 * style="width:95%" |
 * style="width:95%" |
 * style= | (20)
 * }

Stored strain energy in the plate
 * {| style="width:100%" border="0"

$$  \displaystyle U  = \frac{D}{2} \frac{\pi^4 a b}{4} \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} C_{mn}^2 \left(     \frac{m^2}{a^2}      +      \frac{n^2}{b^2}   \right)^2 $$
 * style="width:95%" |
 * style="width:95%" |
 * style= | (21)
 * }

Work done by the external shear force per unit length $$\displaystyle N_{xy}$$
 * {| style="width:100%" border="0"

$$  \displaystyle W  = - 4 N_{xy} \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} \sum_{p=1}^{\infty} \sum_{q=1}^{\infty} C_{mn} C_{pq} \frac {mnpq} {(m^2 - p^2)(n^2 - q^2)} $$
 * style="width:95%" |
 * style="width:95%" |
 * style= | (22)
 * }

Equating the stored strain energy to the work done by external shear load, i.e.,

\displaystyle U = W $$ we obtains the critical buckling shear load $$\displaystyle (N_{xy})_{cr}$$
 * {| style="width:100%" border="0"

$$  \displaystyle (N_{xy})_{cr} =  -   \frac{ab D}{32} \frac {     \displaystyle \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} C_{mn}^2 \left(        \frac{m^2 \pi^2}{a^2}	 +         \frac{n^2 \pi^2}{b^2}      \right)^2 }  {      \displaystyle \sum_{m=1}^{\infty} \sum_{n=1}^{\infty} \sum_{p=1}^{\infty} \sum_{q=1}^{\infty} C_{mn} C_{pq} \frac{mnpq}{(m^2 - p^2)(n^2 - q^2)} } $$
 * style="width:95%" |
 * style="width:95%" |
 * style= | (23)
 * }

The next step is to find the coefficients $$\displaystyle C_{mn}$$ and $$\displaystyle C_{pq}$$ that minimize $$\displaystyle (N_{xy})_{cr}$$ to give $$\displaystyle (N_{xy})_{cr}^\star$$ i.e., the minimum critical load. To do so, take the derivatives of $$\displaystyle (N_{xy})_{cr}$$ with respect to each coefficient $$\displaystyle C_{mn}$$ and set these derivatives to zero; we obtain a set of homogeneous equations with zero right-hand side in terms of the unknown coefficients $$\displaystyle C_{mn}$$.

Introduce the following definitions:
 * {| style="width:100%" border="0"

$$  \displaystyle \vartheta :=  \frac{a}{b} $$ which is the aspect ratio of the plate, and
 * style="width:95%" |
 * style="width:95%" |
 * style= | (24)
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle \lambda :=  -   \frac{\pi^2}{32 \vartheta} \frac{\pi^2 D}{b^2 h (\sigma_{xy})_{cr}} $$
 * style="width:95%" |
 * style="width:95%" |
 * style= | (25)
 * }

The first 5 of these equations read as follows
 * {| style="width:100%" border="0"

$$  \displaystyle \left[ \begin{array}{lllll} \frac{\lambda (1 + \vartheta^2)^2}{\vartheta^2} &	 \frac{4}{9} &	 0	 &	 0	 &	 0	 \\	 \frac{4}{9} &	 \frac{16 \lambda (1 + \vartheta^2)^2}{\vartheta^2} &	 - \frac{4}{5} &	 - \frac{4}{5} &	 \frac{36}{25} \\	 0	 &	 - \frac{4}{5} &	 \frac{\lambda (1 + 9 \vartheta^2)^2}{\vartheta^2} &	 0	 &	 0	 \\	 0	 &	 - \frac{4}{5} &	 0	 &	 \frac{\lambda (9 + \vartheta^2)^2}{\vartheta^2} &	 0	 \\	 0	 &	 \frac{36}{25} &	 0	 &	 0	 &	 \frac{\lambda (9 + 9 \vartheta^2)^2}{\vartheta^2} \end{array} \right] \left\{ \begin{array}{l} C_{11} \\	 C_{22} \\	 C_{13} \\	 C_{31} \\	 C_{33} \end{array} \right\} =  \left\{ \begin{array}{l} 0	 \\	 0	 \\	 0	 \\	 0	 \\	 0     \end{array} \right\} $$ Let $$\displaystyle \mathbf K _{5 \times 5}$$ be the 5x5 matrix shown in Eq.(26). Notice the nice symmetry in $$\displaystyle \mathbf K _{5 \times 5}$$: (1) It is symmetric, (2) On the diagonal coefficients, there is a factor 9 that appears with the index 3 in the coefficients $$\displaystyle C_{mn}$$, e.g., compare ($$\displaystyle K_{33}$$ for $$\displaystyle C_{13}$$) with ($$\displaystyle K_{44}$$ for $$\displaystyle C_{31}$$), and with ($$\displaystyle K_{55}$$ for $$\displaystyle C_{33}$$).
 * style="width:95%" |
 * style="width:95%" |
 * style= | (26)
 * }

Case 1: Truncate the above system in Eq.(26) so to use only the first two equations in Eq.(26) with only two unknowns $$\displaystyle C_{11}$$ and $$\displaystyle C_{22}$$, i.e.,
 * {| style="width:100%" border="0"

$$  \displaystyle \left[ \begin{array}{ll} \frac{\lambda (1 + \vartheta^2)^2}{\vartheta^2} &	 \frac{4}{9} \\	 \frac{4}{9} &	 \frac{16 \lambda (1 + \vartheta^2)^2}{\vartheta^2} \end{array} \right] \left\{ \begin{array}{l} C_{11} \\	 C_{22} \end{array} \right\} =  \left\{ \begin{array}{l} 0	 \\	 0     \end{array} \right\} $$ with the matrix now denoted by $$\displaystyle \mathbf K _{2 \times 2}$$. Set the determinant of $$\displaystyle \mathbf K _{2 \times 2}$$ to zero, i.e.,
 * style="width:95%" |
 * style="width:95%" |
 * style= | (27)
 * }

\displaystyle \det \mathbf K _{2 \times 2} =  0 $$ one finds
 * {| style="width:100%" border="0"

$$  \displaystyle \lambda =  \pm \frac{1}{9} \frac{\vartheta^2}{(1 + \vartheta^2)^2} $$ and thus by Eq.(25), the minimum critical buckling shear stress is
 * style="width:95%" |
 * style="width:95%" |
 * style= | (28)
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle (\sigma_{xy})_{cr}^{\star,2} =  \pm k_c^{\star,2} \frac{\pi^2 D}{b^2 h}  \, \ k_c^{\star,2} =  \frac{\pi^2}{32} \frac{9 (1 + \vartheta^2)^2}{\vartheta^3} $$ where the superscript 2 in $$\displaystyle (\sigma_{xy})_{cr}^{\star,2}$$ simply means that this minimum critical buckling stress was obtained with a 2x2 system, and the $$\displaystyle \pm$$ simply means that the critical buckling load does not depend on the direction of the shear stress (unlike the compressive buckling case). Compare Eq.(29) for $$\displaystyle (\sigma_{xy})_{cr}^{\star,2}$$ to Eq.(8) for $$\displaystyle (\sigma_{xx})_{cr}^{\star}$$.
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (29)
 * }

It was found that using a 2x2 system was not accurate enough; for a square plate, there was a 15% error for $$\displaystyle (\sigma_{xy})_{cr}^{\star,2}$$ when compared to a more accurate estimate when more unknowns were used. This error increases with for increasing aspect ratio $$\displaystyle \vartheta = a/b$$.

Case 2: Now use all 5 equations in Eq.(26) with all 5 unknowns $$\displaystyle \{ C_{11}, C_{22}, C_{13}, C_{31}, C_{33} \}$$. Setting the determinant of $$\displaystyle \mathbf K _{5 \times 5}$$ to zero, i.e.,

\displaystyle \det \mathbf K _{5 \times 5} =  0 $$ one finds
 * {| style="width:100%" border="0"

$$  \displaystyle \lambda =  \left[ \frac{\vartheta^4}{81 (1 + \vartheta^2)^4} \left\{ 1 	 + 	 \frac{81}{625} +	 \frac{81}{25} \left(	   \frac{1 + \vartheta^2}{1 + 9 \vartheta^2}	 \right)^2 +	 \frac{81}{25} \left(	   \frac{1 + \vartheta^2}{9 + \vartheta^2}	 \right)^2 \right\} \right]^{1/2} $$ Using the definition of $$\displaystyle \lambda$$ in Eq.(25), we obtain the minimum critical buckling stress
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (30)
 * }
 * {| style="width:100%" border="0"

$$  \displaystyle (\sigma_{xy})_{cr}^{\star,5} =  \pm k_c^{\star,5} \frac{\pi^2 D}{b^2 h}  \, \ k_c^{\star,5} =  \frac{\pi^2}{32} \frac{1}{\lambda \vartheta} $$
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style="width:92%; padding:10px; border:2px solid #8888aa" |
 * style= | (31)
 * }

Clamped boundary conditions
Use a numerical method such as the finite element method: I am looking for eager students who want to do a Highest-Honors thesis.

Refs: Young & Budynas, Roark's formula for stress and strain, McGraw-Hill, 7th edition, 2002, p.732 (a table with values of $$\displaystyle k_c$$ for $$\displaystyle \vartheta = 1, 2, \infty$$ and a reference to a 1931 technical report). Szilard, Theory and analysis of plates, Prentice Hall, New Jersey, 1974, p.704 (a graph and a reference to a 1964 book and a 1960 report).


 * {| style="width:100%" border="0"

$$  \displaystyle (\sigma_{xy})_{cr} =  K   \frac {E} {1 - \nu^2} \left(     \frac{h}{b}   \right)^2 $$
 * style="width:95%" |
 * style="width:95%" |
 * style= | (35)
 * }

It is also indicated in Young & Budynas [2002] that experimental results yielded a value of $$\displaystyle K = 4.1$$ for large $$\displaystyle \vartheta = a/b$$. The lower experimental value for $$\displaystyle K$$ is likely due to imperfections in the experimental plate specimens.