Solving a Cantilever Beam Problem with MATLAB

[amsyabee]

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.

 

[timthereaper]

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

 

[amsyabee]

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

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).

 

[timthereaper]

Okay, I see what you mean now. Can I see your code or the code you want to modify?

 

[amsyabee]

Okay, I see what you mean now. Can I see your code or the code you want to modify?

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

 

[mustafaornek]

Please help me

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

mustafaornek32@hotmail.com