坚持实干、坚持一线、坚持积累、坚持思考，坚持创新。
接着博文《[力学][有限元][FEM]Basics of Buckling Analysis [曲屈分析基础]》继续介绍桁架单元(Truss Element)用于屈曲分析的方法。该部分内容也是 书本 《有限单元法：编程与软件应用》屈曲分析章节的部分内容节选。
11.3 屈曲分析1：2D桁架
图 11‑2是本章的算例结构，是一榀XZ平面内的桁架结构，结构几何信息、构件属性和约束信息与第二章中静力分析时相同，节点5、6、7受到-z方向的100kN集中力作用。本章将基于2D桁架单元对该结构进行图中荷载样式下的屈曲分析，并将基于MATLAB编程计算的结果与SAP2000、midas Gen分析结果进行对比。
图 11‑2 结构模型示意
以下直接给出小变形情况下桁架单元的几何刚度矩阵
(11.3‑1)
其中P为单元轴力，L为单元长度，单元的几何刚度矩阵与单元的内力有关。
11.3.1 MATLAB代码与注释
屈曲分析的Matlab代码与前面章节中静力分析的Matlab代码略有不同，这里给出进行屈曲分析所需的主要代码。
% Truss 2D Buckling Analysis
% Author ： JiDong Cui（崔济东），Xuelong Shen（沈雪龙）
% Website : www.jdcui.com
% 20170609
……
此前代码与静力分析代码相同，不再赘述。
% Loading
force(10)=-100000;
force(12)=-100000;
force(14)=-100000;
……
根据图 11‑2，为相应的自由度施加荷载，作为参考荷载，用于后续求解单元的几何刚度矩阵。
% Solute balance function
disp('Displacement')
displacement(activeDof) = stiffness( activeDof , activeDof)\force(activeDof)
% Element force
InForce = zeros( numEle , 1 );
for i=1:numEle;
noindex = EleNode(i,:);
deltax = xx( noindex(2)) - xx( noindex(1));
deltay = yy( noindex(2)) - yy( noindex(1));
L = sqrt( deltax*deltax + deltay*deltay );
% Element stiffness matrix
C = deltax / L;
S = deltay / L;
% Corresponding freedom of the element nodes
eleDof = [noindex(1)*2-1 noindex(1)*2 noindex(2)*2-1 noindex(2)*2];
% Element force
InForce(i) = (E/L)*[ -C -S C S ]*displacement(eleDof)*A;
end
disp('Element Stress')
InForce;
该段代码与普通静力分析的代码基本相同，主要用于计算参考荷载作用下单元的内力，用于后续计算单元的几何刚度矩阵，以下将介绍屈曲分析部分的主要代码。
% Geometric Stiffness Matrix
stiffnessM = zeros(numDOF);
stiffnessG = zeros(numDOF);
% Traverse all elements
for i=1:numEle;
% Index of the element nodes
noindex = EleNode(i,:);
% Element length
deltax = xx( noindex(2)) - xx( noindex(1));
deltay = yy( noindex(2)) - yy( noindex(1));
L = sqrt( deltax*deltax + deltay*deltay );
% Element geometric stiffness matrix
C = deltax / L;
S = deltay / L;
T=[ C,S, 0,0;
-S,C, 0,0;
0,0, C,S;
0,0,-S,C
];
eleKm = EA/L*[ C*C C*S -C*C -C*S;
C*S S*S -C*S -S*S;
-C*C -C*S C*C C*S;
-C*S -S*S C*S S*S];
eleKg0=InForce(i)/L*[0,0,0,0;
0,1,0,-1;
0,0,0,0;
0,-1,0,1];
eleKg=T'*eleKg0*T;
% Corresponding freedom of the element nodes
eleDof = [noindex(1)*2-1 noindex(1)*2 noindex(2)*2-1 noindex(2)*2];
% Integrate element geometric stiffness matrix to the global geometric stiffness matrix
stiffnessM( eleDof , eleDof) = stiffnessM( eleDof , eleDof) + eleKm;
stiffnessG( eleDof , eleDof) = stiffnessG( eleDof , eleDof) + eleKg;
end
本段代码主要作用：利用参考荷载作用下的单元内力，求解单元的几何刚度矩阵，进而采用“对号入座”的方式，将单元的几何刚度矩阵叠加，获得结构整体的几何刚度矩阵 。
% Solution
A=-stiffnessM(activeDof,activeDof)\stiffnessG(activeDof,activeDof);
[Vo,Do] = eig(A); % Get the eigenvalues and eigenvectors
fo=diag(Do);
posneg=1./fo.*abs(fo);
Vf=sortrows([Vo',abs(fo),posneg],size(Vo,1)+1);
V=Vf(:,1:size(Vo,1))';
f=Vf(:,size(Vo,1)+1).*Vf(:,size(Vo,1)+2);
fac=1./f
本段代码通过求解标准特征值问题，求得特征值和特征向量。Matlab内置函数eig([A])可直接求解特征值和特征向量，求得的特征值放在向量{Do}中，相应的特征向量放在矩阵[Vo]中。其后的代码用于处理求解得到的数据，其中代码“V=”用于将特征向量排序，代码“f=”用于将结构自振频率进行排序，代码“fac=”用于求得屈曲因子，即特征值的倒数。图 11‑3所示是该桁架结构的前三阶屈曲模态。
(a) Buckling Mode 1
(b) Buckling Mode 2
(c) Buckling Mode 3
- 注释 ( Comments )
( 如果您发现有错误，欢迎批评指正。邮箱：jidong_cui@163.com . 如果您喜欢这篇博文，请在上面给我 点个赞 吧！ 🙂 🙂
( If you found any mistakes in the post, please let me know. Email : jidong_cui@163.com. If you like this posts, please give me a “thumbs up“ rating on the above button! )
- 微信公众号 ( Wechat Subscription)
欢迎关注 “结构之旅” 微信公众号
