推导梁单元的几何刚度矩阵
线性弹性稳定问题，所谓“线性”指的是:①杆的轴向力或板的张力由线性弹性分析决定；②在屈曲引起的无限小位移过程中，轴向力或张力保持不变。对于板来说，就是由线性弹性平面应力分析求得张力，而且在达到屈曲时，张力保持不变。至于非线性屈曲或屈曲后的性态，这将是非线性大位移问题。
▲图1
首先分析单元轴向力在由于横向虚位移引起的虚变形上所作的虚功。图1所示为杆件单元上的任一微段
dx,
AB是其失稳之前的位置,
A'B'为轴压力达到临界值时可能出现的分支平衡位置。设微段由平衡位置
A'B'发生无限小的横向虚位移至
A''B'',
ds1和
ds2分别为虚位移发生前、后微段的长度。微段的虚应变可表示为
\delta \epsilon = \frac{ds2-ds1}{ds1} \quad \cdots (1)
由弧长公式
\begin{split}
ds1 &= \sqrt{1+(\frac{dv}{dx})^2}\approx dx+\frac{1}{2}(\frac{dv}{dx})^2dx\\
ds2 &= \sqrt{1+[\frac{d(v+\delta v)}{dx}]^2}\approx dx+\frac{1}{2}(\frac{dv}{dx}+\frac{d\delta v}{dx})^2dx\\
\end{split} \quad \cdots (2)
(2)代入(1)，并略去高阶微量,得由横向虚位移产生的轴向虚应变为
\delta \epsilon = \frac{dv}{dx} \frac{d\delta v}{dx} \quad \cdots (3)
▲图2
如图2所示，单元轴向荷载所作的虚功为
\delta W_e = F_P \Delta = F_P\int_{x=0}^l\delta \epsilon dx \quad \cdots (4)
式中单元的轴向荷载
F_P以受拉为正。(3)代入(4)可得
\delta W_e = F_P\int_{x=0}^l \frac{dv}{dx} \frac{d\delta v}{dx} dx \quad \cdots (5)
由节点位移插值得挠度
v = \mathbf N \mathbf q^e =[N_1,N_2,N_3,N_4]\begin{Bmatrix}
v_1 \\
\theta_1 \\
v_2 \\
\theta_2 \\
\end{Bmatrix} \quad \cdots (6)
于是
\begin{split}
\frac{dv}{dx} &= \frac{d \mathbf N }{dx}\mathbf q^e \\
\frac{d\delta v}{dx} &= \frac{d \mathbf N }{dx}\delta \mathbf q^e \\
\end{split} \quad \cdots (7)
(7)代入(5)可得
\delta W_e = \delta \mathbf q^e F_P\int_{x=0}^l (\frac{d \mathbf N }{dx})^T \frac{d \mathbf N }{dx} dx \mathbf q^e \quad \cdots (8)
有关形函数，应变矩阵等参见有限元 | 基于虚功原理推导梁单元刚度矩阵
存在于单元中的应力
\sigma 在上述虚应变中所作的虚功为
\begin{split}
\delta W_i & = \int_V \delta \epsilon^T \sigma dV \\
& = \int_V \delta \mathbf q^{eT} \mathbf B^T E \mathbf B \mathbf q^e dV \\
& = \delta \mathbf q^{eT} \int_V \mathbf B^T E \mathbf B dV \mathbf q^e
\end{split} \quad \cdots (9)
由虚功原理
\delta W_e = \delta W_i 可得
\int_V \mathbf B^T E \mathbf B dV \mathbf q^e = F_P\int_{x=0}^l (\frac{d \mathbf N }{dx})^T \frac{d \mathbf N }{dx} dx \mathbf q^e \quad \cdots (10)
记
\mathbf k^e = \int_V \mathbf B^T E \mathbf B dV
\mathbf k_{\sigma}^e = F_P\int_{x=0}^l (\frac{d \mathbf N }{dx})^T \frac{d \mathbf N }{dx} dx
式中，
\mathbf k_{\sigma}^e 与材料物理常数
E无关，只与单元的几何尺寸有关，因此称为几何刚度矩阵。经积分计算后可得
\mathbf k_{\sigma}^e = F_P
\begin{bmatrix}
\frac{6}{5l} & \frac{1}{10} & \frac{-6}{5l} & \frac{1}{10} \\
\frac{1}{10} & \frac{2l}{15} & \frac{-1}{10} & \frac{-l}{30} \\
\frac{-6}{5l} & \frac{-1}{10} & \frac{6}{5l} & \frac{-1}{10}\\
\frac{1}{10} & \frac{-l}{30} & \frac{-1}{10} & \frac{2l}{15}
\end{bmatrix}
(10)可写成
(\mathbf k^e - \mathbf k_{\sigma}^e)\mathbf q^e = 0\quad \cdots (11)
对于一个结构，有
(\mathbf K - \mathbf K_{\sigma})\mathbf q = 0\quad \cdots (12)
该齐次方程组有非零解，则
det(\mathbf K - \mathbf K_{\sigma}) = 0\quad \cdots (13)
(13)可求结构的临界荷载。
