ex1.m

这时一个基于梁元的有限元计算程序

clear
clc
E=210e6;
I=[60 5.69 5.69 5.69 39.5]*1e4;
L=[4.5 27-4.5 56.2-27 86.5-56.2 100-86.5];
P=[1.575 9.45 18.1 28 15.3 4.75]*1e5;
M(1)=0;
M(2)=P(1)*L(1)-P(2)*0.5*L(2);
M(3)=P(3)*0.5*L(2)-P(3)*0.5*L(3);
M(4)=P(4)*0.5*L(3)-P(5)*0.5*L(4);
M(5)=P(5)*0.5*L(4)-P(6)*L(5);
M(6)=0;
n=6;
k1=BeamElementStiffness(E,I(1),L(1));
k2=BeamElementStiffness(E,I(2),L(2));
k3=BeamElementStiffness(E,I(3),L(3));
k4=BeamElementStiffness(E,I(4),L(4));
k5=BeamElementStiffness(E,I(5),L(5));
K=zeros(2*n,2*n);
K=BeamAssemble(K,k1,1,2);
K=BeamAssemble(K,k2,2,3);
K=BeamAssemble(K,k3,3,4);
K=BeamAssemble(K,k4,4,5);
K=BeamAssemble(K,k5,5,6);
K1=K;
m=[3 5 7 9];
for i=1:length(m)
    ki=K1((m(i)+1):end,(m(i)+1):end);
    K1(m(i):(length(ki)+m(i)-1),m(i):(length(ki)+m(i)-1))=ki;
end
k=K1(1:8,1:8);
force=[P(1) M(1) M(2) M(3) M(4) M(5) P(6) M(6)]';
u=jacobi(k,force,[0 0 0 0 0 0 0 0]',0.01);
u=pinv(k)*force;
uk=[u(1) u(2) 10 u(3) 0 u(4) -10 u(5) 0 u(6) u(7) u(8)]';
F=K*uk;
u1=uk(1:4);
u2=uk(3:6);
u3=uk(5:8);
u4=uk(7:10);
u5=uk(9:12);

f1=BeamElementForces(k1,u1);
f2=BeamElementForces(k2,u2);
f3=BeamElementForces(k3,u3);
f4=BeamElementForces(k4,u4);
f5=BeamElementForces(k2,u5);
t=[0 4.5 27 56.2 86.5 100];
box off
subplot(1,2,1);
plot(t,[F(1) F(3) F(5) F(7) F(9) F(11)]');
title('节点受力图');
xlabel('节点位置     /m');
ylabel('F           /N');
subplot(1,2,2);plot(t,[F(2) F(4) F(6) F(8) F(10) F(12)]');
title('节点受力矩图');
xlabel('节点位置    /m');
ylabel('M          /N*m');