有限元方法求结构的临界荷载的最后一步该怎么算？
用有限元方法求结构的临界荷载时，仅通过能量变分或者虚功原理得到
很多参考资料在举例时往往忽略了具体求解过程。本文重点就是如何对(1)进行数值求解。
在(1)中， 称为梁单元的几何刚度矩阵或初应力矩阵，其与轴向力 成正比，与材料无关。以上推导中轴向力取拉为正，轴向力为压力时削弱了梁的抗弯刚度。如果结构的初应力矩阵是按照某个参考载荷 求出的，那么，结构在临界失稳状态时的屈曲载荷为 ，对应的初应力矩阵也随之改变为入 ，则式(1)为
结构失稳的临界载荷对应于结构进入随遇平衡状态，此时在外载荷不变的情况下，结构可由原来的平衡位置转入邻近的平衡位置，因此齐次方程组(2)有非零解，则
在以上推导过程中杆件轴力是以受拉为正的。一般在作刚架稳定性分析时习惯上将轴向力取受压为正。此时，式(2)应改写为
例子1
用有限元法来求解一端固支、一端简支，长度为 的杆的屈曲荷载(划分2个单元)。
▲图1
图1a的有限元模型如图1b所示。采用先处理法，整体结点未知位移列阵如下
得到2个单元弹性刚度矩阵分别为
单元1,2几何刚度矩阵
整体弹性刚度矩阵为
整体几何刚度矩阵为
令 ,则 可得
代入(3)得
展开是关于 的3次方程，解得最小的 ，故临界荷载为
用sympy求解如下
import marimo as mo
import sympy as sp
# marimo类似于jupyter
# 定义符号
μ, l = sp.symbols('μ, l', real = True, positive = True )
l2 = l**2
# 定义矩阵
M = sp.Matrix([[ 24-μ*12/5, 0, 6*l- l*μ/10],
[ 0, 8*l2 - 4*l2*μ/15, 2*l2+l2*μ/30],
[6*l - l*μ/10, 2*l2+l2*μ/30, 4*l2 - 2*l2*μ/15] ])
# 计算行列式
det_M = M.det()
# 解关于μ的三次方程
sol = sp.solve(det_M, μ)
# 数值解
sol_numeric = [s.evalf() for s in sol]
# 格式化输出,f字符串包裹的是markdown格式
mo.md(f"""
**解的结果:**
$$ |M| = {sp.latex(det_M)} $$
$$ μ = {sp.latex(sol_numeric)} $$
""")
例子2
用有限单元法计算图2a所示刚架的临界荷载。
▲图2
图2a的有限元模型如图2b所示。采用先处理法，整体结点未知位移列阵如下
得到3个单元弹性刚度矩阵分别为
由于单元1没有轴向力作用，因此，单元1几何刚度矩阵忽略不计。单元2,3几何刚度矩阵
整体弹性刚度矩阵为
整体几何刚度矩阵为
令 ,则 可得
代入(3)得
展开是关于 的3次方程，解得最小的 ，故临界荷载为
import marimo as mo
import sympy as sp
# marimo类似于jupyter
# 定义符号
μ, l = sp.symbols('μ, l', real = True, positive = True )
l2 = l**2
# 定义矩阵
M = sp.Matrix([[12*l2-μ*4*l2, -24*l+6*l*μ, 4*l2+ l2*μ],
[-24*l+6*l*μ, 192-288*μ, 0 ],
[4*l2+ l2*μ, 0, 16*l2 - 8*l2*μ] ])
# 计算行列式
det_M = M.det()
# 解关于μ的三次方程
sol = sp.solve(det_M, μ)
# 数值解
sol_numeric = [s.evalf() for s in sol]
# 格式化输出,f字符串包裹的是markdown格式
mo.md(f"""
**解的结果:**
$$ |M| = {sp.latex(det_M)} $$
$$ μ = {sp.latex(sol_numeric)} $$
""")
有限元模型有3个自由度，这个方法最终求关于 的3次方程的根。也就是说，有 个自由度，就要解关于 的 次方程的根。
