%-----------------------------------------------------%
%該程序為平面桁架結構分析有限元MATLAB程序
%不能用于計算含傾斜支座的桁架結構
%輸入的信息有節點數、節點信息、單元數、單元信息
%邊界條件處理采用的方法是：對角元素乘大數法
%返回的信息主要有：
%???（1）每個單元的剛度矩陣
%???（2）整體剛度矩陣
%???（3）節點位移向量
%???（4）節點載荷向量
%???（5）每個單元的位移向量
%???（6）每個單元的內力向量
%???（7）每個單元的應力
%
%ver:2017/1/6
%-----------------------------------------------------%

clear
clc

np=3;%輸入節點數
ne=3;%輸入單元數

%node（節點信息）：節點坐標，約束情況，節點載荷
%矩陣規模：np?x?6
node=[
?????0?0?1?1?0?0;
????4?0?0?1?0?0;
????2?3?0?0?5000?-10000;
];

%element（單元信息）：左右節點ij，彈性模量E，截面積A
%矩陣規模：ne?x?4
element=[
????1?2?210e9?0.0001;
????1?3?210e9?0.0001;
????2?3?210e9?0.0001;
];?

kk=zeros（2*np2*np）;%分配整體剛度矩陣內存空間
f?=zeros（2*np1）;%分配載荷內存空間
Elementinfo=zeros（ne3）;%分配空間存儲單元sincosEA
elementstiffnesssave=zeros（44ne）;%分配空間存儲每個單元的剛度矩陣

%首先生成單元剛度矩陣
%然后集成整體剛度矩陣
for?lop=1:ne
????i=element（lop1）;
????j=element（lop2）;
????xi=node（i1）;
????yi=node（i2）;
????xj=node（j1）;
????yj=node（j2）;
????L=sqrt（（xj-xi）*（xj-xi）+（yj-yi）*（yj-yi））;
????s=（yj-yi）/L;
????c=（xj-xi）/L;
????ea=element（lop3）*element（lop4）/L;
????ek?=?ea*[c*cc*s-c*c-c*s;s*cs*s-s*c-s*s;-c*c-c*sc*cc*s;-c*s-s*ss*cs*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???%輸出集成每個單元剛度矩陣后的整體剛度矩陣
end

StiffnessSave=kk;%保存整體剛度矩陣

%邊界條件處理，采用對角元乘大數法
for?lop=1:np
????f（2*lop-11）=f（2*lop-11）?+?node（lop5）;
????f（2*lop??1）=f（2*lop??1）?+?node（lop6）;
????if?node（lop3）?>=?1
????????kk（2*lop-1?2*lop-1）?=?1e8*kk（2*lop-12*lop-1）;
????end
????if?node（lop4）?>=?1
????????kk（2*lop?2*lop）?=?1e8*kk（2*lop?2*lop）;
????end
end

u=kk\f;%高斯消元法求節點位移向量
Force?=?StiffnessSave?*?u;%求節點支反力

Elementforce=zeros（ne1）;%存儲每個桿單元的內力
Elementstress=zeros（ne1）;%存儲每個桿單元的應力
for?lop?=?1:ne
????i=element（lop1）;
????j=element（lop2）;
????E=element（lop3）;
????A=element（lop4）;
????xi=node（i1）;
????yi=node（i2）;
????xj=node（j1）;
????yj=node（j2）;
????L=sqrt（（xj-xi）*（xj-xi）+（yj-yi）*（yj-yi））;
????s=（yj-yi）/L;
????c=（xj-xi）/L;
????%求單元內力
????Elementforce（lop?1）?=E*A/L*（（u（2*j-1）-u（2*i-1））*c+（u（2*j）-u（2*i））*s）;?
????%求單元應力
????Elementstress（lop?1）?=E/L*（（u（2*j-1）-u（2*i-1））*c+（u（2*j）-u（2*i））*s）;
end

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????2871??2017-01-07?12:09??PlaneTrussfem.m

?????文件???????3780??2017-01-08?12:52??PlaneTrussfem_1.m

-----------?---------??----------?-----??----

?????????????????6651????????????????????2