有限元| 刚架临界荷载的新算法

针对前一篇有限元 | 有限元法计算刚架的临界荷载中计算刚架临界荷载的不足之处，提出新的算法。

梁单元新的弹性刚度矩阵和几何刚度矩阵

▲图1

图1所示为杆件单元上的任一微段    ,    是其失稳之前的位置,    为轴压力达到临界值时可能出现的分支平衡位置。设微段由平衡位置    发生无限小的横向虚位移至    。弹性稳定分析是二阶分析，虚位移应该在变形后。    和    分别为虚位移发生前、后微段的长度。微段的轴向虚应变可表示为

由弧长公式

(2)代入(1)，并略去高阶微量,得由横向虚位移产生的轴向虚应变为

在梁单元的4个节点自由度中，    和    的单位不同，现在用

将单位统一起来。其中    是单元长度。梁的挠度

其中

于是

两种坐标系的映射

由(9)(10)得

梁横截面虚应变

其中

▲图2

如图2所示，单元轴向荷载所作的虚功为

式中单元的轴向荷载    以受拉为正。(3)代入(4)可得

内力虚功

由虚功原理    可得

记

则有

上式是一个求广义特征值问题。

经积分计算后可得

算例

▲图3

已知简支梁的长度为    ，划分2个单元，求其临界荷载。

考虑边界条件之后，整体弹性刚度矩阵

整体几何刚度矩阵

由

得

其中    ，在MATLAB中求特征值与特征向量的代码如下

A = [2,3,1,0;
     3,12,0,-3;
     1,0,4,1;
     0,-3,1,2];

B = [4/3, 1, -1/3 0;
     1, 24, 0, -1;
     -1/3, 0, 8/3, -1/3;
     0, -1, -1/3, 4/3];

[X,D] = eig(A,B)
% D是对角矩阵，对角线上的元素是特征值。

D中最小的值是0.1243，则

结果与材料力学相同。特征值0.1243对应得特征向量为    则

便是屈曲模态。