%材料力学课程设计第一题所有问题的求解
%1.以下求解未知力矩
%参考数据q=15300,FB=4400;
clc;clear;
l1=1.1;l2=1.6;l3=3.1;l4=1.6;l5=2.1;FA=2680;
q=input('输入q的值(单位N/m):');
FB=input('输入FB的值(单位N):');
close all;
MC=-FA*l1;MG=-FB*l5;
w2=q*l2^3/12;w3=q*l3^3/12;w4=q*l4^3/12;
A=[2*(l2+l3),l3;l3,2*(l3+l4)];
C=[-6*(w2/2+w3/2)-MC*l2;-6*(w3/2+w4/2)-MG*l4];
RE=A\C;
fprintf('\n1.求解未知力矩的结果:\n')
MD=RE(1);MF=RE(2);
fprintf('MD=%.1fNm , MF=%.1fNm\n\n',MD,MF);
%以下求解支反力,需要调用M函数F_Liang.m
display('第一段梁求解结果:');
[FC,FD1]=F_Liang(FA,0,MC,MD,l2,q);
fprintf('FC=%.1fN , FD1=%.1fN\n\n',FC,FD1);
display('第二段梁求解结果:');
[FD2,FF1]=F_Liang(0,0,MD,MF,l3,q);
fprintf('FD2=%.1fN , FF1=%.1fN\n\n',FD2,FF1);
display('第三段梁求解结果:');
[FF2,FG]=F_Liang(0,FB,MF,MG,l4,q);
fprintf('FF2=%.1fN , FG=%.1fN\n\n',FF2,FG);
FD=FD1+FD2;FF=FF1+FF2;
display('综合得到多跨梁各支座的支反力如下:')
fprintf('FC=%.0fN\nFD=%.0fN\nFF=%.0fN\nFG=%.0fN\n\n\n',FC,FD,FF,FG);
%2.以下画出内力图
%求解剪力方程,画出剪力图
FQ1max=-FA;FQ1min=-FA;
FQ2max=FC-FA;FQ2min=FC-FA-l2*q;
FQ3max=FC-FA-l2*q+FD;FQ3min=FC-FA-l2*q+FD-l3*q;
FQ4max=FC-FA-l2*q+FD-l3*q+FF;FQ4min=FC-FA-l2*q+FD-l3*q+FF-l4*q;
FQ5min=FC-FA-l2*q+FD-l3*q+FF-l4*q+FG;FQ5max=FQ5min;
FQ1=-FA;
FQ2=subs(@(x)FQ2max-q*(x-l1));
FQ3=subs(@(x)FQ3max-q*(x-l1-l2));
FQ4=subs(@(x)FQ4max-q*(x-l1-l2-l3));
FQ5=FQ5max;
disp('2.第一段到第五段梁的剪力方程(单位:N)分别为:')
fprintf('FQ1=%s\n',char(vpa(FQ1,7)))
fprintf('FQ2=%s\n',char(vpa(FQ2,7)))
fprintf('FQ3=%s\n',char(vpa(FQ3,7)))
fprintf('FQ4=%s\n',char(vpa(FQ4,7)))
fprintf('FQ5=%s\n\n',char(vpa(FQ5,7)))
figure(1);
hold on;
ezplot(num2str(FQ1),[0,l1]);
ezplot(FQ2,[l1,l1+l2]);
ezplot(FQ3,[l1+l2,l1+l2+l3]);
ezplot(FQ4,[l1+l2+l3,l1+l2+l3+l4]);
ezplot(num2str(FQ5),[l1+l2+l3+l4,l1+l2+l3+l4+l5]);
xlim([0,l1+l2+l3+l4+l5]);
ylim([-5e4,5e4]);
plot([l1,l1],[FQ1min,FQ2max],'--');
plot([l1+l2,l1+l2],[FQ2min,FQ3max],'--');
plot([l1+l2+l3,l1+l2+l3],[FQ3min,FQ4max],'--');
plot([l1+l2+l3+l4,l1+l2+l3+l4],[FQ4min,FQ5max],'--');
plot([0,l1+l2+l3+l4+l5],[0,0],'k');
text(0.25,FQ1max-0.35e4,num2str(FQ1max));
text(l1-0.45,FQ2max+0.35e4,num2str(FQ2max));
text(l1+l2-0.45,FQ2min-0.35e4,num2str(FQ2min));
text(l1+l2-0.45,FQ3max+0.35e4,num2str(FQ3max));
text(l1+l2+l3-0.45,FQ3min-0.35e4,num2str(FQ3min));
text(l1+l2+l3-0.45,FQ4max+0.35e4,num2str(FQ4max));
text(l1+l2+l3+l4-0.45,FQ4min-0.35e4,num2str(FQ4min));
text(8,FQ5max+0.35e4,num2str(FQ5max));
xlabel('x /m');ylabel('FQ /N');title('剪力图');
%求解弯矩方程,画出弯矩图
syms x;
M1=int(num2str(-FA),0,x);M1r=subs(M1,x,l1);
M2=int(FQ2,l1,x)+M1r;M2r=subs(M2,x,l1+l2);
M3=int(FQ3,l1+l2,x)+M2r;M3r=subs(M3,x,l1+l2+l3);
M4=int(FQ4,l1+l2+l3,x)+M3r;M4r=subs(M4,x,l1+l2+l3+l4);
M5=int(num2str(FQ5max),l1+l2+l3+l4,x)+M4r;M5r=subs(M5,x,l1+l2+l3+l4+l5);
[x2,M2max]=fminbnd(eval(['@(x)',char(-M2)]),l1,l1+l2);M2max=-M2max;
[x3,M3max]=fminbnd(eval(['@(x)',char(-M3)]),l1+l2,l1+l2+l3);M3max=-M3max;
[x4,M4max]=fminbnd(eval(['@(x)',char(-M4)]),l1+l2+l3,l1+l2+l3+l4);M4max=-M4max;
disp('第一段到第五段梁的弯矩方程(单位:Nm)分别为:');
fprintf('M1=%s\n',char(vpa(expand(M1),7)));
fprintf('M2=%s\n',char(vpa(expand(M2),5)));
fprintf('M3=%s\n',char(vpa(expand(M3),6)));
fprintf('M4=%s\n',char(vpa(expand(M4),6)));
fprintf('M5=%s\n\n',char(vpa(expand(M5),7)));
fprintf('剪力图和弯矩图如图所示\n\n\n');
figure(2);
hold on;box off;
ezplot(M1,[0,l1]);
ezplot(M2,[l1,l1+l2]);
ezplot(M3,[l1+l2,l1+l2+l3]);
ezplot(M4,[l1+l2+l3,l1+l2+l3+l4]);
ezplot(M5,[l1+l2+l3+l4,l1+l2+l3+l4+l5]);
plot([0,l1+l2+l3+l4+l5],[0,0],'k');
xlim([0,l1+l2+l3+l4+l5]);
ylim([-1.2e4,1e4]);
plot([l1,l1],[0,M1r],'--');
plot([l1+l2,l1+l2],[0,M2r],'--');
plot([l1+l2+l3,l1+l2+l3],[0,M3r],'--');
plot([l1+l2+l3+l4,l1+l2+l3+l4],[0,M4r],'--');
%text(l1-0.45,FQ2max+0.35e4,num2str(FQ2max));
%text(l1-0.45,M1r-0.6e3,num2str(M1r));
%text(x2-0.3,0.6e3,'num2str(M2max)');
%text(x2,-0.7e3,'\bf\downarrow');
%text(l1+l2-0.3,M2r-0.6e3,num2str(M2r));
%text(x3-0.3,M3max+0.6e3,num2str(M3max));
%text(l1+l2+l3-0.3,M3r-0.6e3,num2str(M3r));
%text(x4-0.3,M4max+0.6e3,num2str(M4max));
%text(l1+l2+l3+l4-0.15,M4r-0.6e3,num2str(M4r));
title('弯矩图');xlabel('x /m');ylabel('Mz /Nm');
%3.以下画出最大弯曲正应力的变化曲线
fprintf('3.最大弯曲正应力的变化曲线如图所示\n\n\n');
h1=0.1;h2=0.12;h3=0.11;
b1=0.06;b2=0.08;b3=0.07;
t=0.005;E=210e9;
Wz=@(b,h)(b*h^3-(b-2*t)*(h-2*t)^3)/(6*h);
Wz1=Wz(b1,h1);
Wz2=Wz(b2,h2);
Wz3=Wz(b3,h3);
sigma1=abs(M1/Wz1);
sigma2=abs(M2/Wz2);
sigma3=abs(M3/Wz2);
sigma4=abs(M4/Wz2);
sigma5=abs(M5/Wz3);
figure(3);
hold on;
ezplot(sigma1,[0,l1]);
ezplot(sigma2,[l1,l1+l2]);
ezplot(sigma3,[l1+l2,l1+l2+l3]);
ezplot(sigma4,[l1+l2+l3,l1+l2+l3+l4]);
ezplot(sigma5,[l1+l2+l3+l4,l1+l2+l3+l4+l5]);
plot([0,l1+l2+l3+l4+l5],[0,0],'k');
xlim([0,l1+l2+l3+l4+l5]);
ylim([0,2.5e8]);
Y12=max(subs(sigma1,x,l1),subs(sigma2,x,l1));
Y23=max(subs(sigma4,x,l1+l2+l3+l4),subs(sigma5,x,l1+l2+l3+l4));
plot([l1,l1],[0,Y12],'--');
plot([l1+l2+l3+l4,l1+l2+l3+l4],[0,Y23],'--');
xlabel('x /m');ylabel('\sigma_{max} /Pa');title('最大弯曲正应力变化曲线');
%4.以下求最大挠度
%求各段挠度
syms C1 D1 C2 D2 C3 D3 C4 D4 C5 D5;
v2=int(int(M2*2/(Wz2*h2*E)))+C2*x+D2;
[C2,D2]=solve(subs(v2,x,l1),subs(v2,x,l1+l2));
v3=int(int(M3*2/(Wz2*h2*E)))+C3*x+D3;
[C3,D3]=solve(subs(v3,x,l1+l2),subs(v3,x,l1+l2+l3));
v4=int(int(M4*2/(Wz2*h2*E)))+C4*x+D4;
[C4,D4]=solve(subs(v4,x,l1+l2+l3),subs(v4,x,l1+l2+l3+l4));
theta2=subs(diff(v2,x));
theta2C=subs(theta2,x,l1);
theta4=subs(diff(v4,x));
theta4G=subs(theta4,x,l1+l2+l3+l4);
v1=int(int(M1*2/(Wz1*h1*E)))+C1*x+D1;
theta1=subs(diff(v1,x));
theta1C=subs(theta1,x,l1);
[C1,D1]=solve(theta1C-theta2C,subs(v1,x,l1));
v5=int(int(M5*2/(Wz3*h3*E)))+C5*x+D5;
theta5=subs(diff(v5,x));
theta5G=subs(theta5,x,l1+l2+l3+l4);
[C5,D5]=solve(theta5G-theta4G,subs(v5,x,l1+l2+l3+l4));
v1=subs(v1);
v2=subs(v2);
v3=subs(v3);
v4=subs(v4);
v5=subs(v5);
%画出挠曲线
figure(4);
hold on;box off;
ezplot(v1,[0,l1]);
ezplot(v2,[l1,l1+l2]);
ezplot(v3,[l1+l2,l1+l2+l3]);
ezplot(v4,[l1+l2+l3,l1+l2+l3+l4]);
ezplot(v5,[l1+l2+l3+l4,l1+l2+l3+l4+l5]);
plot([0,l1+l2+l3+l4+l5],[0,0],'k');
xlim([0,l1+l2+l3+l4+l5]);
ylim([-0.05,0.01]);
xlabel('x /m');ylabel('v /m');title('挠曲线图');
[v1x,v1min]=fminbnd(eval(['@(x)',char(v1)]),0,l1);
[v2x,v2max]=fminbnd(eval(['@(x)',char(-v2)]),l1,l1+l2);v2max=-v2max;
[v3x,v3min]=fminbnd(eval(['@(x)',char(v3)]),l1+l2,l1+l2+l3);
[v4x,v4max]=fminbnd(eval(['@(x)',char(-v4)]),l1+l2+l3,l1+l2+l3+l4);v4max=-v4max;
[v5x,v5min]=fminbnd(eval(['@(x)',char(v5)]),l1+l2+l3+l4,l1+l2+l3+l4+l5);
text(0.2,v1min-0.003,num2str(v1min));
text(v2x,v2max+0.003,num2str(v2max));
text(v3x,v3min-0.003,num2str(v3min));
text(v4x,v4max+0.003,num2str(v4max));
text(v5x-1,v5min-0.003,num2str(v5min));
fprintf('4.挠曲线图如图所示\n\n');
%用能量法求各段最大挠度
M_1=x;
M_21=-(l1+l2-v2x)/l2*(x-l1);M_22=(v2x-l1)/l2*(x-l1-l2);
M_31=-(l1+l2+l3-v3x)/l3*(x-l1-l2);M_32=(v3x-l1-l2)/l3*(x-l1-l2-l3);
M_41=-(l1+l2+l3+l4-v4x)/l4*(x-l1-l2-l3);M_42=(v4x-l1-l2-l3)/l4*(x-l1-l2-l3-l4);
M_5=-x+l1+l2+l3+l4+l5;
f11=int(2*M1*M_1/(E*Wz1*h1),x,0,l1);f10=-theta1C*l1;f1=f11+f10;
f21=int(2*M2*M_21/(E*Wz2*h2),x,l1,v2x);f22=int(2*M2*M_22/(E*Wz2*h2),x,v2x,l1+l2);f2=f21+f22;
f31=int(2*M3*M_31/(E*Wz2*h2),x,l1+l2,v3x);f32=int(2*M3*M_32/(E*Wz2*h2),x,v3x,l1+l2+l3);f3=f31+f32;
f41=int(2*M4*M_41/(E*Wz2*h2),x,l1+l2+l3,v4x);f42=int(2*M4*M_42/(E*Wz2*h2),x,v4x,l1+l2+l3+l4);f4=f41+f42;
f51=int(2*M5*M_5/(E*Wz3*h3),x,l1+l2+l3+l4,l1+l2+l3+l4+l5);f50=theta4G*l5;f5=f50+f51;
f1=double(subs(f1));f2=double(subs(f2));f3=double(subs(f3));f4=double(subs(f4));f5=double(subs(f5));
A=[abs(f1),abs(f2),abs(f3),abs(f4),abs(f5)];
fmax=max(A);
num=find(fmax==A);
fmax=eval(['f',num2str(num)]);
fmax_x=eval(['v',num2str(num),'x']);
fprintf('各段的最大挠度为:\nf1max= %.5fm\nf2max= %.6fm\nf3max= %.5fm\nf4max= %.5fm\nf5max= %.4fm\n\n',f1,f2,f3,f4,f5);
fprintf('综