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

clear
clc

np=3;%輸入節(jié)點數(shù)
ne=3;%輸入單元數(shù)

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

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

kk=zeros（2*np2*np）;%分配整體剛度矩陣內(nèi)存空間
f?=zeros（2*np1）;%分配載荷內(nèi)存空間
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;%保存整體剛度矩陣

%邊界條件處理，采用對角元乘大數(shù)法
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;%高斯消元法求節(jié)點位移向量
Force?=?StiffnessSave?*?u;%求節(jié)點支反力

Elementforce=zeros（ne1）;%存儲每個桿單元的內(nèi)力
Elementstress=zeros（ne1）;%存儲每個桿單元的應(yīng)力
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;
????%求單元內(nèi)力
????Elementforce（lop?1）?=E*A/L*（（u（2*j-1）-u（2*i-1））*c+（u（2*j）-u（2*i））*s）;?
????%求單元應(yīng)力
????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