固定端梁的振型和固有频率计算-matlab

function FF_end_BEAM(~) % Fixed-Fixed (Clamped-Clamped) end beam evaluations % Mode shapes and natural frequency calculations % HELP: This script is to evaluate mode shapes and corresponding natural % frequencies of the fixed-fixed end beam by a user specified mechanical % properties and geometry size of the beam. % Prepare the followings: % - Material properties of the beam, viz. density (Ro), Young's modulus (E) % - Specify a cross section of the beam, viz. square,rectangular, circular % - Geometry parameters of the beam, viz. Length, width, thickness % % Fixed-Fixed (Clamped-Clamped) beam calculations % E Young's modulus % Ro density % Ix cross sectional moment of inertia (along the X axis) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % by Sulaymon L. ESHKABILOV, Ph.D % October, 2011; Revised/Updated- August, 2024 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PART I. fprintf('What is the X-section of your beam? \n') fprintf('If it is a Circle, Enter 1. If square, Enter 2. \n') fprintf('If it is a Rectangle, Enter 3: \n') fprintf('If your beam"s X-section is not listed here, Enter 4. \n') fprintf('To see example # 1 (retangular X-section), hit [ENTER] or Enter 0: \n') fprintf('To see example # 2 (retangular X-section), Enter 5: \n') fprintf('To use specific beam"s X-section including, E, Ix, Ro, enter 6 \n') CS=input(' Enter your choice: '); if isempty(CS) || CS==0 fprintf('Example #1. Rectangular X-section Aluminum beam \n') fprintf('Length=0.321 [m], Width=0.05 [m], Thickness=0.006 [m]; \n') fprintf('E=69.9*1e9 [Pa]; Ro=2770 [kg/m^3] \n') fprintf(' \n') L=.321; W=.05; Th=.006; A=W*Th; Ix=(1/12)*W*Th^3; E=69.9e+9; Ro=2770; elseif CS==1 R=input('Enter Radius of the X-section in [m]: '); L=input('Enter Length in [m]: '); Ix=(1/4)*pi*R^4; A=pi*R^2; YA=Material_Selection; if YA==1 E=input('Enter Young"s modulus in [Pa]: '); Ro=input('Enter material density in [kg/m^3]: '); else Material_Properties; E=input('Enter Young"s modulus in [Pa]: '); Ro=input('Enter material density in [kg/m^3]: '); end elseif CS==2 W=input('Enter Width of the X-section in [m]: '); L =input('Enter Length in [m]: '); Ix=(1/12)*W^4; A=W^2; YA=Material_Selection; if YA==1 E=input('Enter Young"s modulus in [Pa]: '); Ro=input('Enter material density in [kg/m^3]: '); else Material_Properties; E=input('Enter Young"s modulus in [Pa]: '); Ro=input('Enter material density in [kg/m^3]: '); end elseif CS==3 W=input('Enter Width of the X-section in [m]: '); Th=input('Enter Thickness of the X-section in [m]: '); L=input('Enter Length in [m]: '); Ix=(1/12)*W*Th^3; A=W*Th; YA=Material_Selection; if YA==1 E=input('Enter Young"s modulus in [Pa]: '); Ro=input('Enter material density in [kg/m^3]: '); else Material_Properties; E=input('Enter Young"s modulus in [Pa]: '); Ro=input('Enter material density in [kg/m^3]: '); end elseif CS==4 fprintf('Note: you need to compute Ix (area moment of inertia along x axis) \n') fprintf('and X-sectional area \n') L=input('Enter Lengthin [m]: '); Ix=('Enter Ix in [m^4]: '); A=('Enter X-sectional area in [m^2]: '); YA=Material_Selection; if YA==1 E=input('Enter Young"s modulus in [Pa]: '); Ro=input('Enter material density in [kg/m^3]: '); else Material_Properties; E=input('Enter Young"s modulus in [Pa]: '); Ro=input('Enter material density in [kg/m^3]: '); end elseif CS==5 fprintf('Display Example #2 \n') fprintf('This is a rectangular X-sectional steel beam. \n') fprintf('Length=0.45 m; Width=0.04 m; Thickness=0.003 m; \n') L=.45; W=.04; Th=.003; A=W*Th; Ix=(1/12)*W*Th^3; E=2.1*1e11; Ro=7.85*1e3; else F=warndlg('It is not clear what your choice of X-section of a beam is. Re-run the script and enter your beam"s X-section !!!','!! Warning !!'); waitfor(F) fprintf('Type in the Command Window after >> FF_end_BEAM \n') pause(3) return end %% PART II. fprintf('How many modes and mode shapes would you like to evaluate ? \n') HMMS=input('Enter the number of modes and mode shapes to compute: '); if HMMS>=7 fprintf(' \n') warning('NOTE: Up to 6 mode shapes (plots) are displayed via the script. Yet, using evaluated data (Xnx) of the script, more mode shapes can be plotted'); fprintf(' \n') end Nm=4*HMMS; jj=4; while jj<=Nm betaNL(jj)=fzero(@(betaNL)cosh(betaNL)*cos(betaNL)-1,jj); jj=jj+4; end index=betaNL~=0; betaNLall=betaNL(index); betaN=(betaNLall/L)'; k=1; wn=zeros(1,length(betaN)); fn=ones(1,length(wn)); while k<=length(betaN) wn(k)=betaN(k)^2*sqrt((E*Ix)/(Ro*A)); fn(k)=wn(k)/(2*pi); fprintf('Mode shape # %d corresponds to nat. freq (fn): %3.1f [Hz] \n', k, fn(k)) k=k+1; end x=linspace(0, L, 180); xl=x./L; sigmaN=(cosh(betaNLall)-cos(betaNLall))./(sinh(betaNLall)-sin(betaNLall)); Tc='(cosh(betaN(ii).*x(jj))-cos(betaN(ii).*x(jj)))-sigmaN(ii).*(sinh(betaN(ii).*x(jj))-sin(betaN(ii)*x(jj)))'; Xnx=zeros(length(betaN),length(x)); for ii=1:length(betaN) for jj=1:length(x) Xnx(ii,jj)=eval(Tc); end end % Plot mode shapes; fprintf('Upto 6 mode shapes can be displayed \n'); fprintf(' Yet, using evaluated data (Xn(x)), more mode shapes can be plotted \n'); MMS=HMMS; if MMS==1 plot(xl,Xnx(1,:), 'b-', 'LineWidth', 2) title('Mode shape of the Fixed-fixed (Clamped-Clamped)') legend('Mode #1', 'location', 'best'); xlabel('x/L'); ylabel('Mode shape X_n(x)'); grid hold off; shg elseif MMS==2 plot(xl,Xnx(1,:), 'b-', 'LineWidth', 2); hold on plot(xl,Xnx(2,:), 'r-', 'LineWidth', 2); grid title('Mode shapes of the Fixed-fixed (Clamped-Clamped)') legend('Mode #1', 'Mode #2', 'location', 'best') xlabel('x/L'); ylabel('Mode shape X_n(x)') hold off; shg elseif MMS==3 plot(xl,Xnx(1,:), 'b-', 'LineWidth', 2); hold on plot(xl,Xnx(2,:), 'r-', 'LineWidth', 2) plot(xl,Xnx(3,:), 'm-'); grid title('Mode shapes of the Fixed-fixed (Clamped-Clamped)') legend('Mode #1', 'Mode #2', 'Mode #3', 'location', 'best') xlabel('x/L'); ylabel('Mode shape X_n(x)') hold off; shg elseif MMS==4 plot(xl,Xnx(1,:), 'b-', 'LineWidth', 2); hold on plot(xl,Xnx(2,:), 'r-', 'LineWidth', 2) plot(xl,Xnx(3,:), 'm-', 'LineWidth', 2) plot(xl,Xnx(4,:), 'k-', 'LineWidth', 2); grid title('Mode shapes of the Fixed-fixed (Clamped-Clamped)') legend('Mode #1', 'Mode #2', 'Mode #3', 'Mode #4', 'location', 'best') xlabel('x/L'); ylabel('Mode shape X_n(x)') hold off; shg elseif MMS==5 plot(xl,Xnx(1,:), 'b-', 'LineWidth', 2); hold on plot(xl,Xnx(2,:), 'r-', 'LineWidth', 2) plot(xl,Xnx(3,:), 'm-', 'LineWidth', 2) plot(xl,Xnx(4,:), 'k-', 'LineWidth', 2) plot(xl,Xnx(5,:), 'g-', 'LineWidth', 2) grid