Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Draw the mode shapes and get the natural frequencies of the cantilever beam(with a fo

  1. Jun 21, 2011 #1
    hi~I'm Rex.I have a problem,can anyone help me slove it.
    The problem is : Draw the mode shapes and get the natural frequencies of the cantilever beam(with a force in free end)
    The bar is shown in Fig.1(you can see in )
    and it elastic modulus of 200GPa, cross-sectional area of 0.001m^2, and density of 7860 kg/m^3. P=P0*sinΩ*t applied at point xi,wherei=1,2,3,4,5(five nodes) and respectively(P0=100 ,Ω=30 and t=10(s), t is time,P is force). the problem by commercial code such as matlab. I have a matlab code(Euler–Bernoulli beam),This is the matlab code(Euler-Bernoulli beam). http://www.mathworks.com/matlabcent...of-the-cantilever-beam-with-a-force-in-free-e You can see it there,and try to run in matlab.how can i modify it and transform to Cantilever beam with end load. If you interested this question,we can discuss it. My english is poor,i hope you understand the question which i asked.Nice to meet you.
     

    Attached Files:

    • 1_n.jpg
      1_n.jpg
      File size:
      5.6 KB
      Views:
      326
    Last edited: Jun 21, 2011
  2. jcsd
  3. Jun 21, 2011 #2
    Re: Draw the mode shapes and get the natural frequencies of the cantilever beam(with

    First off, your English isn't bad at all. I understood everything you were saying. I don't, however, understand the purpose of the 1,2,3,4 in boxes at the bottom of the diagram. I know the 1,2,3,4,5 at the top are for the nodes. I tried clicking on the link you provided and it said the link is no longer valid. Hopefully you can find the current link and then I'll be able to look at the code.

    If you wanted to approximate the cantilever beam as a mass-spring-damper system and get the natural frequencies of the system, here is a link that will help:

    http://faculty.uml.edu/pavitabile/22.403/web_downloads/Final_Project_Cantilever_101806.pdf
     
  4. Jun 21, 2011 #3
    Re: Draw the mode shapes and get the natural frequencies of the cantilever beam(with

    The purpose of the 1,2,3,4 in boxes at the bottom of the diagram is element.So the bar is with four elements(five nodes).We can try to put the force in this nodes,like in free end point(the fifth node).
     
  5. Jun 21, 2011 #4
    Re: Draw the mode shapes and get the natural frequencies of the cantilever beam(with

    Okay, I see what you mean now. Can I see your code or the code you want to modify?
     
  6. Jun 21, 2011 #5
    Re: Draw the mode shapes and get the natural frequencies of the cantilever beam(with

    this code is in Euler-Bernoulli beam. how to modify it and transform to a cantilever beam with a force in the fifth node(free end point).
    Main program:
    %--------------------------------------------------------------------------
    % Calculating the natural frequencies of a cantilever beam
    % length = 35 cm, width = 10 cm, thickness = 1cm
    %
    %
    % /│
    % /│=============================
    % /│ EI │
    % /│=============================
    % /│
    %
    %
    % PURPOSE :1.Get the natural frequencies of the Euler-Bernoulli beam.
    % 2.Draw the mode shapes of the Euler-Bernoulli beam.
    %
    % REFERENCE:1.K matrix assembly refer to DAVID S. BURNETT, "Finite
    % Element Analysis," 1987, AT&T Bell Laboratories, p363-370.
    % 2.M & K matrices assembly refer to J.N.REDDY,"An Introduction
    % to the Finite Element Method",1993,McGraw-Hill,p150-154,
    % p225-227,p237
    % Edit: RenYu-Chen
    % DATE : 12/9/2004
    %--------------------------------------------------------------------------
    clear all;clf;
    format long e;
    %
    % --------------------------------------------------------------------------
    % Input data of the cantilever beam
    % --------------------------------------------------------------------------
    % UNIT:M.K.S.
    E=200*10^9; % Young's modulus:N/m^2
    Wb=10/100; % beam width:m
    Tb=1/100; % beam thickness:m
    Lb=35/100; % beam length:m
    pb=7860; % beam density (per unit volume):kg/m^3
    I=Wb*Tb^3/12; % moment of inertia:m^4
    A=Wb*Tb; % cross section area of the beam:m^2
    %
    % --------------------------------------------------------------------------
    % Defining the elements and nodes
    % --------------------------------------------------------------------------
    % where we use 4 elements
    N_elem=4; % elements
    node=zeros(N_elem+1,2); % nodes
    for i=1:N_elem+1
    node(i,1)=i;
    node(i,2)=Lb/N_elem*(i-1);
    end
    NUM_NODE=length(node); % the size of nodes
    NUM_ELEM=length(node)-1; % the size of elements
    matrix_size=2*NUM_ELEM+2;
    %
    % --------------------------------------------------------------------------
    % Initializes K and F matrice
    % --------------------------------------------------------------------------
    %
    M=zeros(matrix_size,matrix_size);
    K=zeros(matrix_size,matrix_size);
    %
    % --------------------------------------------------------------------------
    % Assembling of K and M matrices
    % --------------------------------------------------------------------------
    %
    % (1). assembly of K matrix
    ELNO=0; % ELNO:the ith element
    for ii=1:2:matrix_size-3
    ELNO=ELNO+1;
    [KE]=k_elemu(node(ELNO,2),node(ELNO+1,2),Lb,E,I);
    [K]=mtxasmby(K,KE,ii); % matrix assembly by calling "mtxasmby.m"
    end % end of FOR loop -- assembly of K matrix
    %
    % (2). assembly of M matrix
    ELNO=0; % ELNO:the ith element
    for ii=1:2:matrix_size-3
    ELNO=ELNO+1;
    [ME]=m_elemu(node(ELNO,2),node(ELNO+1,2),Lb,pb,A);
    [M]=mtxasmby(M,ME,ii); % matrix assembly by call mtxasmby.m
    end % end of FOR loop -- assembly of M matrix
    %
    % --------------------------------------------------------------------------
    % Imposing the essential Boundary Conditions
    % --------------------------------------------------------------------------
    % In this case, one end of the Euler-Bernoulli beam is fixed.
    % The other end is free. Therefore, the K and M matrices can be deleted
    % 2 rows and 2 columns. In other words, we can delete the firstly and
    % secondly two rows and columns of the matrices.
    %
    % K_bc: apply BCs to the K_matrix
    % M_bc: be reduced orders as apply BCs to the K_matrix
    %
    K_bc=zeros(matrix_size-2,matrix_size-2);
    M_bc=zeros(matrix_size-2,matrix_size-2);
    %
    K_bc=K(3:matrix_size,3:matrix_size);
    M_bc=M(3:matrix_size,3:matrix_size);
    %
    % --------------------------------------------------------------------------
    % Calculate eigenvalues by the finite element formulation
    % --------------------------------------------------------------------------
    %
    ei=eig(K_bc,M_bc); % eigenvalues
    % sorted natural angular frequencies [rad/s]
    ef=sort(real(sqrt(ei)));
    % sorted natural angular frequencies [Hz]
    f=ef/(2*pi);
    %
    % --------------------------------------------------------------------------
    % Compare the first six modal testing and F.E.M values of frequencies with
    % theoretical values.
    % --------------------------------------------------------------------------
    %
    % Theoretical calculated values are expressed in xbar(i) = lambda(i)*L
    % variables and were copied here
    beta= [1.875104068706770e+000 ...
    4.694091132933031e+000 ...
    7.854757438070675e+000 ...
    1.099554073487850e+001 ...
    1.413716839104652e+001 ...
    1.727875953327674e+001 ];
    % first six frequencies of modal testing
    modal_text_f= [ 28.9 ...
    180 ...
    506 ...
    992 ...
    1640 ...
    2440 ];
    % Modal testing calculated values are expressed in xbar_ex(i) = lambda_ex(i)*L
    % variables and were copied here
    beta_ex= sqrt(modal_text_f*2*pi/sqrt(E*I/(pb*A*Lb^4)));
    % auxiliary constants
    c0=sqrt(E/pb);
    j=sqrt(I/A);
    b2=beta.*beta;
    b2_ex=beta_ex.*beta_ex;
    %
    om=b2*c0*j/(Lb^2)/(2*pi); % theoretical calculated angular frequencies [Hz]
    om_ex=b2_ex*c0*j/(Lb^2)/(2*pi); % modal test calculated angular frequencies [Hz]
    %
    % --------------------------------------------------------------------------
    % plot first six theoretical, modal testing and FEM natural frequencies
    % --------------------------------------------------------------------------
    %
    ix=1:6;
    figure(1); subplot(2,1,1)
    plot( ix,om(ix),'r-', ix, om_ex(ix),'b*-',ix,f(ix),'ko')
    title('Frequencies \omega','fontsize',18);
    xlabel('counter','fontsize',18);
    ylabel('angular frequencies','fontsize',18);
    legend('Theoretical', 'Modal texting','FEM');
    % --------------------------------------------------------------------------
    % compute and plot relative errors for the first 6 frequencies
    % --------------------------------------------------------------------------
    r=zeros(size(ix));
    for i=ix
    r(i)=100*(f(i)-om(i))/om(i);
    r_ex(i)=100*(abs(om_ex(i)-om(i)))/om(i);
    end
    figure(1); subplot(2,1,2)
    plot(ix,r_ex,'b*-',ix,r,'ro-')
    title('relative errors [%]','fontsize',18);
    xlabel('counter','fontsize',18);
    ylabel('[%]','fontsize',18);
    legend('relative errors for Modal testing','relative errors for FE');
    % print results for the first 6 frequencies
    disp(' ======================= relative errors [%] ======================= ')
    disp(' counter, Theoretical frequencies, Modal testing frequencies, FE frequencies, relative errors for Modal test[%], relative errors for FE[%]')
    disp([ix' om(ix)' om_ex(ix)' f(ix) r_ex' r'])
    %
    % --------------------------------------------------------------------------
    % Draw the first six Modal shape
    % --------------------------------------------------------------------------
    %
    % modes(1:6): the first six modes in Hertz
    % modeshapes(6, 10) : the vibration mode shapes of the first six modes
    %
    modes=zeros(6,1);
    modeshapes=zeros(6,N_elem);
    modes_ex=zeros(6,1);
    modeshapes_ex=zeros(6,N_elem);
    modes_fe=zeros(6,1);
    modeshapes_fe=zeros(6,N_elem);
    % loop over the six modes
    ix=1:6;
    beta_fe=sqrt(f(ix)'*2*pi/sqrt(E*I/(pb*A*Lb^4)));
    for i=ix;
    modes(i) = beta(i)^2 * sqrt(E*I/(pb*A*Lb^4));
    modes(i) = modes(i)/(2*pi);
    modes_ex(i) = beta_ex(i)^2 * sqrt(E*I/(pb*A*Lb^4));
    modes_ex(i) = modes_ex(i)/(2*pi);
    modes_fe(i) = beta_fe(i)^2 * sqrt(E*I/(pb*A*Lb^4));
    modes_fe(i) = modes_fe(i)/(2*pi);
    % coefficients for computing mode shape
    betaL=beta(i);
    betaL_ex=beta_ex(i);
    betaL_fe=beta_fe(i);
    %
    a1 = sin(betaL) - sinh(betaL);
    a2 = cos(betaL) + cosh(betaL);
    a1_ex = sin(betaL_ex) - sinh(betaL_ex);
    a2_ex = cos(betaL_ex) + cosh(betaL_ex);
    a1_fe = sin(betaL_fe) - sinh(betaL_fe);
    a2_fe = cos(betaL_fe) + cosh(betaL_fe);
    %
    x = 0;
    x_ex=0;
    x_fe=0;
    increment = Lb/N_elem;
    % compute modeshapes for j =1:5,
    for j =1:N_elem
    y = a1*(sin(x) - sinh(x));
    y = y +a2*(cos(x) - cosh(x));
    x = x + (beta(i)/Lb)*increment;
    modeshapes(i,j) = y;
    %
    y_ex = a1_ex*(sin(x_ex) - sinh(x_ex));
    y_ex = y_ex +a2_ex*(cos(x_ex) - cosh(x_ex));
    x_ex = x_ex + (beta_ex(i)/Lb)*increment;
    modeshapes_ex(i,j) = y_ex;
    %
    y_fe = a1_fe*(sin(x_fe) - sinh(x_fe));
    y_fe = y_fe +a2_fe*(cos(x_fe) - cosh(x_fe));
    x_fe = x_fe + (beta_fe(i)/Lb)*increment;
    modeshapes_fe(i,j) = y_fe;
    % increment the beam span by 5 intervals for plotting
    beamspan=(1:N_elem)*(1/N_elem)*Lb;
    %
    % plot the first six mode shape;
    for s=1:6
    subplot(2, 3, s)
    plot(figure(s))
    figure(i)
    ymax =max(abs(modeshapes(i,:)));
    ymax_ex =max(abs(modeshapes_ex(i,:)));
    ymax_fe =max(abs(modeshapes_fe(i,:)));
    plot( beamspan, modeshapes(i,:)/ymax,'r-', beamspan, modeshapes_ex(i,:)/ymax_ex,'b-.', beamspan, modeshapes_fe(i,:)/ymax_fe,'go');
    grid on;
    ylabel('modeshape amplitude','fontsize',5);
    xlabel('beam span [m]','fontsize',5);
    title(['Cantilever Beam, Mode ',num2str(i)],'fontsize',5);
    end
    end;
    end
    % end of porgram

    K Stiffness matrice is Subroutine ,should put this code in matlab work folder)

    % k_elemu
    function [KE]=K_elemu(prenode,postnode,Lb,E,I)
    % PURPOSE : This is a subprogram as for Stiffness matrice
    %
    % K=[ 12 6*lb -12 6*lb
    % 4*lb^2 -6*lb 2*lb^2
    % 12 -6**lb
    % symmetric 4*lb^2 ];
    %
    l=postnode-prenode;
    lb=l/Lb;
    % the length of present element per unit length, lb:length_bar
    KE(1,1)=12 ;
    KE(2,1)=-6*lb ; KE(2,2)=4*lb^2 ;
    KE(3,1)=-12 ; KE(3,2)=6*lb ; KE(3,3)=12 ;
    KE(4,1)=-6*lb ; KE(4,2)=2*lb^2 ; KE(4,3)=6*lb ; KE(4,4)=4*lb^2;
    KE=KE/lb^3*E*I/Lb^3;
    for i=2:4
    for j=1:i-1
    KE(j,i)=KE(i,j);
    end
    end

    M Mass matrice is Subroutine ,should put this code in matlab work folder)

    % m_elemu
    function [ME]=m_elemu(prenode,postnode,Lb,pb,A)
    % PURPOSE : This is a subprogram as for Mass matrice
    %
    % M=1/420[ 156 22*lb 54 -13*lb
    % 4*lb^2 13*lb -3*lb^2
    % 156 -22*lb
    % symmetric 4*lb^2];
    %
    l=postnode-prenode;
    lb=l/Lb;
    % the length of the present element per unit length, lb:length_bar
    ME(1,1)=156 ;
    ME(2,1)=22*lb ; ME(2,2)=4*lb^2 ;
    ME(3,1)=54 ; ME(3,2)=13*lb ; ME(3,3)=156 ;
    ME(4,1)=-13*lb ; ME(4,2)=-3*lb^2; ME(4,3)=-22*lb ; ME(4,4)=4*lb^2 ;
    ME=ME*lb/420*pb*A*Lb;
    for i=2:4
    for j=1:i-1
    ME(j,i)=ME(i,j);
    end
    end
     
  7. Jun 24, 2011 #6
    Re: Draw the mode shapes and get the natural frequencies of the cantilever beam(with

    Please help me

    I want to write matlab codes of multi-cracked beam vibrations.

    mustafaornek32@hotmail.com
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Draw the mode shapes and get the natural frequencies of the cantilever beam(with a fo
Loading...