%%?薄板算例（直接法——平面應力問題）
close?all
clear?all
clc

%%?基本參數(shù)定義
numele=4;
numnode=5;
elenode=3;
L=1;
t=1;%單元厚度
elenodecorre=[1?2?3;1?3?4;2?5?3;3?5?4];
E=206;%合金鋼，單位Gpa
mu=0.25;%合金鋼，0.24-0.28
q=20;?%均布拉應力
%%?施加載荷前作圖
%inicoornode=[0?0;1?0;1/2?1/2;0?1;1?1];
ininodevector=[0?0?1?0?1/2?1/2?0?1?1?1];
figure
axis?off
axis?equal
hold?on
for?inicoori=1:numnode
????plot（ininodevector（2*inicoori-1）ininodevector（2*inicoori）‘bo‘）;
end
for?inicoorj=1:numele
????hold?on
????line（[ininodevector（2*elenodecorre（inicoorj1）-1）?...
????????ininodevector（2*elenodecorre（inicoorj2）-1）]...
????????[ininodevector（elenodecorre（inicoorj1）*2）?...
????????ininodevector（elenodecorre（inicoorj2）*2）]）;
????hold?on
????line（[ininodevector（2*elenodecorre（inicoorj2）-1）?...
????????ininodevector（2*elenodecorre（inicoorj3）-1）]...
????????[ininodevector（elenodecorre（inicoorj2）*2）?...
????????ininodevector（elenodecorre（inicoorj3）*2）]）;
????hold?on
????line（[ininodevector（2*elenodecorre（inicoorj3）-1）?...
????????ininodevector（2*elenodecorre（inicoorj1）-1）]...
????????[ininodevector（elenodecorre（inicoorj3）*2）?...
????????ininodevector（elenodecorre（inicoorj1）*2）]）;
end
%hold?on

%%??
sysmatK=zeros（1010）;
u=zeros（81）;
delt=zeros（101）;
f=zeros（101）;


for?elei=1:numele
????%?initical?coordinates?of?the?nodes?of?each?element
????x1=ininodevector（elenodecorre（elei1）*2-1）;
????y1=ininodevector（elenodecorre（elei1）*2）;
????x2=ininodevector（elenodecorre（elei2）*2-1）;
????y2=ininodevector（elenodecorre（elei2）*2）;
????x3=ininodevector（elenodecorre（elei3）*2-1）;
????y3=ininodevector（elenodecorre（elei3）*2）;????
????%?interpolation?function?（shape?function）
????area=0.5*（x2*y3+x1*y2+x3*y1-x1*y3-x2*y1-x3*y2）;
????
????%?對于每個節(jié)點進行計算
????%N1=0.5*（x2*y3-x3*y2+（y2-y3）*x+（x3-x2）*y）/area;
????%N2=0.5*（x3*y1-x1*y3+（y3-y1）*x+（x1-x3）*y）/area;
????%N3=0.5*（x1*y2-x2*y1+（y1-y2）*x+（x2-x1）*y）/area;
????b1=y2-y3;
????b2=y3-y1;
????b3=y1-y2;
????c1=x3-x2;
????c2=x1-x3;
????c3=x2-x1;
????B=[b1?0?b2?0?b3?0;0?c1?0?c2?0?c3;...
????????;c1?b1?c2?b2?c3?b3]/（2*area）;%由形函數(shù)偏導而來的單元應變矩陣
????D=[1?mu?0;mu?1?0;0?0?（1-mu）/2]*E/（1-mu^2）;%狀態(tài)彈性矩陣
????K=t*B‘*D*B;?%單元剛度矩陣
????A=area*B‘*t;%節(jié)點應力與單元邊上應力之間的關系矩陣
????S=D*B;%單元應力矩陣（單元位移向量與應力向量的轉化關系）
????inditemp=elenodecorre（elei:）;
????indicator=[2*inditemp（1）-1??2*inditemp（1）?2*inditemp（2）-1?...?
????????2*inditemp（2）?2*inditemp（3）-1??2*inditemp（3）];
????sysmatK（indicatorindicator）=K+sysmatK（indicatorindicator）;?%系統(tǒng)剛度矩陣組裝

end

f（3）=q*L*t;
f（9）=q*L*t;
soluK=sysmatK（[3?4?5?6?9?10][3?4?5?6?9?10]）;
soluf=f（[3?4?5?6?9?10]）;
soludelt=inv（soluK）*soluf;
delt（[3?4?5?6?9?10]）=soludelt;
delt
%%?施加載荷后作圖
finnodevector=ininodevector‘+delt

hold?on
for?fincoori=1:numnode
????plot（finnodevector（2*fincoori-1）finnodevector（2*fincoori）‘ro‘）;
end
for?fincoorj=1:numele
????hold?on
????line（[finnodevector（2*elenodecorre（fincoorj1）-1）?...
????????finnodevector（2*elenodecorre（fincoorj2）-1）]...

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----
?????文件????????2862??2019-01-05?09:48??有限元平面應力及桁架機構matlab代碼\Truss.m
?????文件????????3830??2019-01-03?16:46??有限元平面應力及桁架機構matlab代碼\thin_plate.m
?????文件??????????79??2019-01-11?18:54??有限元平面應力及桁架機構matlab代碼\文件說明.txt
?????目錄???????????0??2019-01-11?18:52??有限元平面應力及桁架機構matlab代碼\