%-----------------------------------------------------%
%该程序为平面桁架结构分析有限元MATLAB程序
%不能用于计算含倾斜支座的桁架结构
%输入的信息有节点数、节点信息、单元数、单元信息
%程序中采用的单位是国际标准单位
%边界条件处理采用的方法是:对角元素置1法
%返回的信息主要有:
% (1)每个单元的刚度矩阵
% (2)整体刚度矩阵
% (3)节点位移向量
% (4)节点载荷向量
% (5)每个单元的位移向量
% (6)每个单元的内力向量
% (7)每个单元的应力
%
%ver:2017/1/6
%-----------------------------------------------------%
clear
clc
np = 3 ; %输入节点数
ne = 2 ; %输入单元数
%node(节点信息):节点坐标,约束情况,节点载荷(不包含支反力)
%矩阵规模:np x 6
% node=[
% 0 0 1 1 0 0;
% 0.4 0 0 1 20000 0;
% 0.4 0.3 0 0 0 -30000;
% 0 0.3 1 1 0 0;
% ];
node=[0 0 1 1 0 0;
0.5 0.5*(3)^(1/2) 1 0 0 -10000;
1 0 1 1 0 0;];
% node=[0 0 1 1 0 0;
% 0.4 0 0 1 20000 0;
% 0.4 0.3 0 0 0 -25000;
% 0 0.3 1 1 0 0;];
%element(单元信息):左右节点ij,弹性模量E,截面积A
%矩阵规模:ne x 4
% element=[
% 1 2 300e9 0.0001;
% 2 3 300e9 0.0001;
% 1 3 300e9 0.0001;
% 4 3 300e9 0.0001;
% ] ;
element=[1 2 210e9 0.0001;2 3 210e9 0.0001];
% element=[1 2 295e9 0.0001;
% 2 3 295e9 0.0001;
% 1 3 295e9 0.0001;
% 4 3 295e9 0.0001;];
kk=zeros(2*np,2*np);%分配整体刚度矩阵内存空间
f =zeros(2*np,1);%分配载荷内存空间
elementstiffnesssave=zeros(4,4,ne);%分配空间存储每个单元的刚度矩阵
kksave=zeros(2*np*ne,2*np);%用于保存每集成一个单元后的整体刚度矩阵
%首先生成单元刚度矩阵
%然后集成整体刚度矩阵
for lop=1:ne
i=element(lop,1);
j=element(lop,2);
xi=node(i,1);
yi=node(i,2);
xj=node(j,1);
yj=node(j,2);
L=sqrt((xj-xi)*(xj-xi)+(yj-yi)*(yj-yi));
s=(yj-yi)/L;
c=(xj-xi)/L;
TT=[c*c,c*s,-c*c,-c*s;s*c,s*s,-s*c,-s*s;-c*c,-c*s,c*c,c*s;-c*s,-s*s,s*c,s*s]
ea=element(lop,3)*element(lop,4)/L;
ek = ea*[c*c,c*s,-c*c,-c*s;s*c,s*s,-s*c,-s*s;-c*c,-c*s,c*c,c*s;-c*s,-s*s,s*c,s*s];
ek %输出每个单元的刚度矩阵
elementstiffnesssave(:,:,lop)=ek;
kk(2*i-1:2*i, 2*i-1:2*i)=kk(2*i-1:2*i, 2*i-1:2*i)+ek(1:2, 1:2);
kk(2*i-1:2*i, 2*j-1:2*j)=kk(2*i-1:2*i, 2*j-1:2*j)+ek(1:2, 3:4);
kk(2*j-1:2*j, 2*i-1:2*i)=kk(2*j-1:2*j, 2*i-1:2*i)+ek(3:4, 1:2);
kk(2*j-1:2*j, 2*j-1:2*j)=kk(2*j-1:2*j, 2*j-1:2*j)+ek(3:4, 3:4);
kk %输出集成每个单元刚度矩阵后的整体刚度矩阵
kksave((lop-1)*2*np+1:lop*2*np, :)=kk;%保存每集成每个单元刚度矩阵后的整体刚度矩阵
end
%将保存每集成每个单元刚度矩阵后的整体刚度矩阵
% document=xlswrite('temp.xlsx',kksave,'A1:H32');
StiffnessSave=kk;%保存整体刚度矩阵
%边界条件处理,采用对角元素置1法
for lop=1:np
f(2*lop-1,1)=f(2*lop-1,1) + node(lop,5);
f(2*lop ,1)=f(2*lop ,1) + node(lop,6);
if node(lop,3) >= 1
kk(2*lop-1, :)=0;
kk(:, 2*lop-1)=0;
kk(2*lop-1, 2*lop-1) = 1;
f(2*lop-1,1)=0;
end
if node(lop,4) >= 1
kk(2*lop, :)=0;
kk(:, 2*lop)=0;
kk(2*lop, 2*lop) = 1;
f(2*lop,1)=0;
end
end
u=kk\f; %高斯消元法求节点位移向量
Force = StiffnessSave * u; %求节点支反力
Elementforce=zeros(ne,1);%存储每个杆单元的内力
Elementstress=zeros(ne,1);%存储每个杆单元的应力
for lop = 1:ne
i=element(lop,1);
j=element(lop,2);
E=element(lop,3);
A=element(lop,4);
xi=node(i,1);
yi=node(i,2);
xj=node(j,1);
yj=node(j,2);
L=sqrt((xj-xi)*(xj-xi)+(yj-yi)*(yj-yi));
s=(yj-yi)/L;
c=(xj-xi)/L;
%单元的位移列向量
emove = [u(2*i-1,1);u(2*i,1);u(2*j-1,1);u(2*j,1);];
%求单元内力
Elementforce(lop, 1) =E*A/L* [-c -s c s] * emove;
%求单元应力
Elementstress(lop, 1) =E/L*[-c -s c s] * emove;
end