质量块简支梁.rar_alivem67_梁计算_模态 梁_简支梁_简支梁模态

%%%%%%%%%%%%%%%%%%简支梁%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%简支梁%%%%%%%%%%%%%%%%% clear all,close all clc L=5; %长度 b=0.05; %宽 h=0.01; %厚 单位 米 I=(b*h^3)/12; E=21*10^10; EI=E*I; A=b*h; %横截面积 rho=7850; %密度 y1=0.2*L;%第一条裂纹位置 y2=0.5*L;%第二条裂纹位置 r1=0; %裂纹深度 %裂纹模拟 f1=1.8624*r1^2-3.95*r1^3+16.375*r1^4-37.226*r1^5+76.81*r1^6 ... -126.9*r1^7+172*r1^8-143.97*r1^9+66.56*r1^10; c1=5.346*h*f1; %无量纲固有频率求解： r=[]; for i =1:50 a=(i-1)+eps; b=i*1; fhandle=@f_SS_1crack; [fa,Ta] = feval(fhandle,a,L,c1,y1,y2); [fb,Tb] = feval(fhandle,b,L,c1,y1,y2); options = optimset; if (fa > 0) == (fb < 0) || (fa < 0) == (fb > 0) r(i)=fzero(fhandle,[a b],options,L,c1,y1,y2); end end beta=nonzeros(r); omega=sqrt(EI*beta.^4/(rho*A*L^4));%圆频率 %系数求解 for N_modal=1:6 w=beta(N_modal); [f,T]=fhandle(w,L,c1,y1,y2); A=T(1:11,1:11); B=T(1:11,12); WS{N_modal}=-inv(A)*B; end %模态 p=100; %取点数 xx=linspace(0,L,p); for as=1:p if xx(as)<=y1, K1=as; end if xx(as)<=y2, K2=as; end %模态 end for jj=1:6 %模态阶数，1到6阶 %beam-1 const(:,jj)=WS{jj}; s=xx(1:K1); Y1(:,jj)=const(1,jj)*cosh(beta(jj)*(s/L))+const(2,jj)*sinh(beta(jj)*(s/L))+const(3,jj)*cos(beta(jj)*(s/L))+const(4,jj)*sin(beta(jj)*(s/L)); %beam-2 s=xx(K1+1:K2); Y2(:,jj)=const(5,jj)*cosh(beta(jj)*(s/L))+const(6,jj)*sinh(beta(jj)*(s/L))+const(7,jj)*cos(beta(jj)*(s/L))+const(8,jj)*sin(beta(jj)*(s/L)); %beam-3 s=xx(K2+1:p); Y3(:,jj)=const(9,jj)*cosh(beta(jj)*(s/L))+const(10,jj)*sinh(beta(jj)*(s/L))+const(11,jj)*cos(beta(jj)*(s/L))+sin(beta(jj)*(s/L)); Y(1:K1,jj)=Y1(:,jj); Y(K1+1:K2,jj)=Y2(:,jj); Y(K2+1:p,jj)=Y3(:,jj); end %画出模态形状： figure(1) plot(xx,Y(:,1),'k-'),hold on plot(xx,Y(:,2),'k-.'),hold on plot(xx,Y(:,3),'k--'),hold on plot(xx,Y(:,4),'k:') xlabel('x') ylabel('Y','rotation',0) legend('1st','2nd','3rd','4th','Location','Best','orientation','Horizontan') legend('boxoff') title('简支梁前四阶模态振型') figure(2); subplot(2,2,1); plot(xx,Y(:,1)); xlabel('x'); ylabel('Y','rotation',0); title('1st mode shape'); subplot(2,2,2); plot(xx,Y(:,2)); xlabel('x'); ylabel('Y','rotation',0); title('2nd mode shape'); subplot(2,2,3); plot(xx,Y(:,3)); xlabel('x'); ylabel('Y','rotation',0); title('3rd mode shape'); subplot(2,2,4); plot(xx,Y(:,4)); xlabel('x'); ylabel('Y','rotation',0); title('4th mode shape'); %损伤位置： figure(3) for i=1:98 YY(i+1)=Y(i,1)+Y(i+2,1)-2*Y(i+1,1); end plot(xx(2:99),YY(2:99)); xlabel('x') ylabel('Y','rotation',0); YYYY=Y(:,1); qqq=YYYY; w=2; N = length(qqq); d = zeros(1,w); dh = zeros(1,w); D = zeros(1,N-w); for i=1:(N-w) s = qqq(i:i+w); for j=1:w d(j) = sqrt(((s(j+1) - s(1))^2) + (((j+1)-1)^2)); dh(j) = sqrt(((s(j+1) - s(j))^2) + (((j+1)-j)^2)); end dmax =max(d); L = sum(dh); D(i) = log10(w)/(log10(w) + log10(dmax/L)); end figure(4); plot(D); xlswrite('ceshi.xls',D',2);