FEM-Course-Matlab/1.matlab平面纯弯梁有限元编程/固支梁/beam_main.m

49 lines
1.2 KiB
Mathematica
Raw Permalink Normal View History

2024-01-28 16:46:36 +00:00
clear all;close all;clc;
%<EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>mm N MPa
E=210e3;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD>
t=2;%<EFBFBD><EFBFBD><EFBFBD>߶<EFBFBD>
b=100;%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
I=1/12*100*t^3;%<EFBFBD><EFBFBD><EFBFBD>Ծ<EFBFBD>
L=1;%<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
TL=500;%<EFBFBD><EFBFBD><EFBFBD>ܳ<EFBFBD><EFBFBD><EFBFBD>
ele_num=TL/L;%<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
for i=1:ele_num
ke(:,:,i)=BeamElementStiffness(E,I,L);%%%<EFBFBD><EFBFBD>Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ƶ<EFBFBD>
end
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȫ<EFBFBD>ָնȾ<EFBFBD><EFBFBD><EFBFBD>
K=zeros((ele_num+1)*2);
for i=1:ele_num
K=BeamAssemble(K,ke(:,:,i),i,i+1); %ȫ<EFBFBD>ָնȾ<EFBFBD><EFBFBD><EFBFBD>
end
k=[K(3:(ele_num+1)*2-2,3:(ele_num+1)*2-2)]; %ȥ<EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD>ĸնȾ<EFBFBD><EFBFBD><EFBFBD>
f=zeros(2*(ele_num-1),1); %<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>f
f(ele_num-1)=-400; %<EFBFBD>м<EFBFBD><EFBFBD>ڵ<EFBFBD>ʩ<EFBFBD><EFBFBD>400N<EFBFBD><EFBFBD><EFBFBD>µļ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
u=k\f; %%%% inv(k)*f <EFBFBD><EFBFBD>˹<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
U=zeros((ele_num+1)*2,1);%<EFBFBD>ڵ<EFBFBD>λ<EFBFBD><EFBFBD>
U(3:ele_num*2)=u;%
for i=1:ele_num
Fe(i,:)=ke(:,:,i)*U((i-1)*2+1:(i-1)*2+4); %%%%<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>=<EFBFBD><EFBFBD>Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>*<EFBFBD>ڵ<EFBFBD>λ<EFBFBD><EFBFBD>
end
F=K*U; %%%%ȫ<EFBFBD>ָնȾ<EFBFBD><EFBFBD><EFBFBD>*<EFBFBD>ڵ<EFBFBD>λ<EFBFBD><EFBFBD>=<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܵ<EFBFBD><EFBFBD>IJ<EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
figure;
x = [0:L:500]';
M=Fe(:,2);
M(ele_num+1)=-Fe(end,4);
hold on;
title('Bending Moment Diagram');
plot(x,M);
y1 = zeros(1,ele_num+1);
plot(x,y1,'k')
xlabel('length(mm)')
ylabel('Moment(N.mm)')
title('moment diagram with constant thickness');
figure;
hold on;
title('Displacment Diagram');
for i=1:length(x)
disp_plot_constant(i)=U(2*(i-1)+1);
end
plot(x,disp_plot_constant);
save data disp_plot_constant