題：用MATLAB編制有限元程序

? 鋼軌長度60米，扣件間距0.6米，鋼軌參數：密度7800Kg/m3，彈性模量2.1×1011N/m2，慣性矩3.09×10-5m4，截面積為77.45cm2。軌下膠墊剛度為20KN/mm，膠墊阻尼損耗因子為0.25，計算（1）跨中受到200KN豎直向下的力時，鋼軌的內力分布；（2）跨中受到200KN豎直向下的力時，鋼軌的內力分布。

%材料參數&幾何參數

E=2.1*10^8;A=77.45*10^-4;I=3.09*10^-5;L=0.6;q=0.592;rub=20000

%平面梁單元剛度矩陣ke

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]

%將單元剛度矩陣集結為總體剛度矩陣

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

%再計入膠墊剛度系數

for ?i=2:3:302

? ? ? ? K(i,i)=K(i,i)+rub;

end

return

%構建節(jié)點荷載向量（等效節(jié)點荷載）

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é)點位移向量U，U為一個303×1的列陣

U=pinv(K)*F;

fjd=-rub*[U(2:3:302,1)];%計算各膠墊（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');