Draw the mode shapes and get the natural frequencies of the cantilever beam(with a force in free end)

[EDIT: 20110621 11:15 CDT - merge comment into question, clarify - WDR]

hi, I'm Rex. I have a problem, can anyone help me solve 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 my blog. <http://www.facebook.com/profile.php?...type=1&theater#!/photo.php?fbid=224903064204253&set=o.21222876640&type=1&theater> ) and it has 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, where i=1,2,3,4,5(five nodes) and respectively (P0=100 ,Ω=30 and t=10(s), t is time, P is force).

[I need to solve] 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.de/matlabcentral/answers/9519-draw-the-mode-shapes-and-get-the-natural-frequencies-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. I can mail the code to you,tell me your mail please.

If you interested this question,we can discuss it. This is my mail:[e-mail]

My english is poor, i hope you understand the question which i asked. Nice to meet you.

Answer 1

[Chen Rex]

%--------------------------------------------------------------------------
%   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_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_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

this code is in Euler-Bernoulli beam. how to modify it and transform to a cantilever beam with a force in the five nodes point(free end point).

[Xiaojian]

I tried to run this code, but it looks something is missing, can anyone tell me why? where can i find the missing code?

"Error in ==> FEMbeam at 69

[K]=mtxasmby(K,KE,ii); % matrix assembly by calling "mtxasmby.m""

[Simon Streicher]

I found the following fix to work:

replace

%[K]=mtxasmby(K,KE,ii);  % matrix assembly by calling "mtxasmby.m"
K((ELNO*2-1):(ELNO+1)*2,(ELNO*2-1):(ELNO+1)*2)=K((ELNO*2-1):(ELNO+1)*2,(ELNO*2-1):(ELNO+1)*2)+KE;

and replace

%[M]=mtxasmby(M,ME,ii);  % matrix assembly by call mtxasmby.m
M((ELNO*2-1):(ELNO+1)*2,(ELNO*2-1):(ELNO+1)*2)=M((ELNO*2-1):(ELNO+1)*2,(ELNO*2-1):(ELNO+1)*2)+ME;

[AMIRRUL]

>> cantileverbeam

Error: File: cantileverbeam.m Line: 288 Column: 6

All functions in a script must be closed with an 'end'.

How to fix it? Can anyone help me?

Answer 2

[Sean de Wolski]

Aren't the eigenvectors the buckling modes? ~ Corresponding to the (ascending) sorted eigenvalues?

To change it to a cantilever, change your boundary conditions, possibly your discretization to include any important points, and your force vector.

[Chen Rex]

I don't know how to modify the code with essential Boundary Conditions?? Can you modify this code or write a new code??

tell me the force vector please!

M and K matrice is 10rows and 10columns.

MX"+KX=f(t),f(t) is P0=100sin30t.

The boundary conditions is like this web. http://emweb.unl.edu/Mechanics-Pages/Scott-Whitney/325hweb/Beams.htm

but how can i change my code??

[Sean de Wolski]

I can modify this code for you. I could also write a new code that's more consolidated. I'll do it for your simple cantilever beam once I have $10,000.00 American wired to a bank account that I'll have to set up just for this transaction.

Seriously, you're COPYING someone's well commented, well documented code and you're not going to take the time to learn how it works or why it works? As a civil/structural engineer myself, I wouldn't want to drive across any beam you design.

[Chen Rex]

Sorry!! I'm try to modify it and slove this problem.The code editor is RenYu-Chen. I just study how to tansform the code to cantilever beam with a force in free end point. I promise that i will not make money or trading for this code.But i'm not a engineer because i'm a student.I am trying to fix the vibration of higher learning.Help me and i am very grateful to you.

[Sean de Wolski]

Yes. So what is more useful for you as a student:

-Figuring out what RenYu did yourself: what they're doing, how/why they're doing it and what modifications you need to make it work for your problem.

or

-Having one of us do you homework assignment for you?

[Chen Rex]

Do not want to be too complicated, just a question. I just want to see my problem is how to solve other people's. I looked for books related to the discussion, but no information on the cantilever beam with a concentrated force of the code.I know you can modify this code.

If you do not want to help me I would not mind, thank you reply.

[Sean de Wolski]

I have no problem helping you. I will not do you homework for you. Helping would be something like: "Why does the author overwrite K_bc with a new K_bc? Doesn't that make the previous line useless?" or "What is the last double for-loop doing? How could I consolidate it?"

Helping is not: "Here's a bunch of code, now can you change it so it meets the specs of my homework assignment so I can get back to lounging?"

The first example of helping is the point of this community - people with more experience assisting people with less experience become more proficient. The second example is example is the point of contractors - pay people to do your work for you.

[Chen Rex]

I understand your intention. so,can you help me try to modify it??How to create force matrix in this code??the boundary conditions is?

What is the difference in Left point and right point??

What is the difference in Euler-Bernoulli beam code or cantilever beam code??

[Sean de Wolski]

If the code is well written there is no difference; the engine is the same just the parameters change. K_bc,M_bc, F.

F has to be a vector of forces corresponding to each DOF. Most will be zero. Figure out which DOF has your load applied and place it at that point in a column vector of zeros (zeros everywhere there is no force)

[Chen Rex]

so, the F vector is [0 0 0 0 0 0 0 0 P0]?? if the force put in free end point(fifth noide).one nod have two DOF.How to combine the matrix and M,K,F in the code???Euler-Bernoulli beam without considering the force matrix.Can you try to write three element code for me,and i try to do four element in myself??

[Sean de Wolski]

I would expect there to be three DOF per node: x,y,theta.

So for a vertical load at the 5th node if there are five nodes total:

F = zeros(15,1);

F(14) = P; %or negative P if it's down.

[Chen Rex]

Can not just consider the situation in two degrees of freedom:x,y?

F = zeros(15,1);What does it mean??

F(14) = P,Meaning is like this [0;0;0;0;0;0;0;0;0;0;0;0;0;P;0]??

Can i input

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

F=f(t) % force

t=time % time

and output mode shapes??

[Sean de Wolski]

Two DOF corresponds to a truss element. You need three for a beam/frame.

[Chen Rex]

Oh!In this case(cantiliver beam with a force),can you write the code of beam with 3 nodes(two elements).

[Chen Rex]

Are you offline??

[Walter Roberson]

Sean is a student and has, I am sure, numerous things to do.

[Sean de Wolski]

It's a gorgeous day to go water-skiing Walter !!

[Chen Rex]

wow!!It sounds pretty good. It's fun!! Nice to meet you Walter Roberson.

Do you and Sean de have time to help me solve this problem?

Answer 3

[Adjal yassine]

Salam alikom , hello I would like to tell you that the important thing is to understand to finite element method then it is easy to understand this program with Matlab , try to understand theory then the program, there is good book in this field , Matlab codes for finite element method and you will find many examples also there are many anlaytical methods for analysis the vibration of beam, one of this method is finite element method , what I can do with this method , we can give an apporiximate solution