题：用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

%构建节点荷载向量（等效节点荷载）

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;

%求解节点位移向量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’);