題：用MATLAB編制有限元程序
? 鋼軌長度60米，扣件間距0.6米，鋼軌參數(shù)：密度7800Kg/m3泣矛，彈性模量2.1×1011N/m2，慣性矩3.09×10-5m4禾蚕，截面積為77.45cm2您朽。軌下膠墊剛度為20KN/mm，膠墊阻尼損耗因子為0.25换淆，計(jì)算（1）跨中受到200KN豎直向下的力時(shí)哗总，鋼軌的內(nèi)力分布几颜；（2）跨中受到200KN豎直向下的力時(shí)，鋼軌的內(nèi)力分布魂奥。
%材料參數(shù)&幾何參數(shù)
E=2.1*10^8;A=77.45*10^-4;I=3.09*10^-5;L=0.6;q=0.592;rub=20000
%平面梁單元?jiǎng)偠染仃噆e
ke=[E*A/L 0 0 -E*A/L 0 0;
0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^36*E*I/L^2;
0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^22*E*I/L;
-E*A/L 0 0 E*A/L 0 0;
0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3-6*E*I/L^2;
0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^24*E*I/L]
%將單元?jiǎng)偠染仃嚰Y(jié)為總體剛度矩陣
K=zeros(303);
? for i=1:100
? ? ? ? ?for j=1:6
? ? ? ? ? ? ? ? for k=1:6
? ? ? ? ? ? ? ? ? ?K(3*(i-1)+j,3*(i-1)+k)=K(3*(i-1)+j,3*(i-1)+k)+ke(j,k);
? ? ? ? ? ? ? ? ?end
? ? ? ? ? ?end
? ?end
%再計(jì)入膠墊剛度系數(shù)
for ?i=2:3:302
? ? ? ? K(i,i)=K(i,i)+rub;
end
return
%構(gòu)建節(jié)點(diǎn)荷載向量（等效節(jié)點(diǎn)荷載）
F=zeros(303,1);
for i=2:3:302
? ? ? ? ?F(i,1)=-q*L;
end
return
F(2,1)=-q*L/2;F(302,1)=-q*L/2;
F(152,1)=-q*L-200;
F(3,1)=-q*L^2/12;
F(303,1)=q*L^2/12;
%求解節(jié)點(diǎn)位移向量U，U為一個(gè)303×1的列陣
U=pinv(K)*F;
fjd=-rub*[U(2:3:302,1)];%計(jì)算各膠墊（jd）的支反力
%繪制剪力圖
for i=1:100
? ? ? if i<51
? ? ? ?ffl=sum(fjd(1:i)); holdon
? ? ? ?plot([(i-1)*L,i*L],[0,0],'k-');hold on
? ? ? ?plot([(i-1)*L,(i-1)*L],[0,ffl-(i-1)*q*L]);hold on
? ? ? ? plot([i*L,i*L],[0,ffl-i*q*L]);hold on
? ? ? ? plot([(i-1)*L,i*L],[ffl-(i-1)*q*L,ffl-i*q*L]);holdon;
else
? ? ? ?ffl=sum(fjd(1:i))-200;holdon
? ? ? ?plot([(i-1)*L,i*L],[0,0],'k-');holdon
? ? ? ?plot([(i-1)*L,(i-1)*L],[0,ffl-(i-1)*q*L]);holdon
? ? ? ?plot([i*L,i*L],[0,ffl-i*q*L]);holdon
? ? ? ? plot([(i-1)*L,i*L],[ffl-(i-1)*q*L,ffl-i*q*L]);holdon;
? ?end
end
title('剪力圖');xlabel('長度/m');ylabel('剪力/KN');
%繪制彎矩圖
syms x;
for i=1:100
? ? ? if i<51
? ? ? ? M=-q/2*x^2;
?for ?j=1:i
? ? ? ?M=M+fjd(j)*(x-(j-1)*L);
end
fplot(-M,[(i-1)*L, i*L]); hold on;
else
? ? M=-q/2*x^2-200*(x-30);
for ?j=1:i
? ? ? M=M+fjd(j)*(x-(j-1)*L);
end
fplot(-M,[(i-1)*L, i*L]); hold on;
? ?end
end
plot([0,60],[0,0],'k-');
axis([0 60 -4510]);
title('彎矩圖');xlabel('長度/m');ylabel('彎矩/KN*m');
