Program 14.
6: MATLAB program to find the
frequency of a Timoshenko beam
TIMOSHENKOVIB
% dynamics of Timoshenko beam by finite element method
clc;
ne=5;
nj=ne+1;
neq=2*nj;
K=zeros(neq,neq);
M=zeros(neq,neq);
% give nodi and nodj of each member
for i=1:nj
nodi(i)=i;
nodj(i)=i+1;
end
%give the values of e,g,a,i and lengths of beam
e=210e9;
g=70e9;
a=[1.85187e-5 1.85187e-5 1.85187e-5 1.85187e-5 1.85187e-5 1.85187e-5] ;
mi=[2.85785e-11 2.85785e-11 2.85785e-11 2.85785e-11 2.85785e-11...
2.85785e-11];
angle=0;
L=10;
for i=1:ne
l(i)=L/ne;
end
% give density of the material
rho=7800.0;
%give axial load of the member
P=200;
%give timoshenko shear constant ko=5/6 for rect ko=9/10 for circular
ko=5/6;
% number of constraint degrees of freedom
nbou=2;
% the numbers of constrained degrees of freedom
nb=[1 2*nj-1];
% form 6 x 6 element stiffness and mass matrix and assemble wilson method
for n=1:ne
i=nodi(n);
j=nodj(n);
phi(n)=12.0*e*mi(n)/(ko*a(n)*g*l(n)^2);
k=TimoshenkoElementStiffness(e,a(n),mi(n),l(n),P,phi(n));
Finite element method in relation to structural dynamics 519
m=TimoshenkoElementMass(rho,a(n),mi(n),l(n),phi(n));
K=TimoshenkoAssemble(K,k,i,j);
M=TimoshenkoAssemble(M,m,i,j);
end
% apply boundary conditions using wilson method
for i=1:nbou
ii=nb(i);
for j=1:neq
K(ii,j)=0.0;
K(j,ii)=0.0;
M(ii,j)=0.0;
M(j,ii)=0.0;
end
K(ii,ii)=1;
M(ii,ii)=1;
end
% find inv(K)*M
invk=inv(K);
km=invk*M;
format long;
% find the eigen values and mode shapes of inv(K)*M
[ms,ns]=size(M);
% %eigen values and eigen vectors
[evec,ev]=eig(km);
for i=1:ms
ee(i)=1/ev(i,i);
end
Qh=max(ee)+0.001;
Ql=0;
for i=1:ms
for j=1:ms
if ee(j) > Ql & ee(j) < Qh
kl=j;
Qh=ee(j);
else
end
end
Ql=Qh;
Qh=max(ee)+0.001;
om1(i)=ee(kl);
omega(i)=sqrt(ee(kl));
for lm=1:ms
p1(lm,i)=evec(lm,kl);
520 Structural dynamics of earthquake engineering
end
end
%Normalizing the mode shape
LL=p1'*M*p1;
%develop modal matrix
for i=1:ms
for j=1:ms
p(i,j)=p1(i,j)/LL(j,j);
end
end
disp( Natural frequencies in rad/sec)
disp(omega)
disp( normalized modal vector )
disp(p)
function y = TimoshenkoElementStiffness(E,A,I,L,P,phi)
%TimoshenkoElementStiffness This function returns the element
% stiffness matrix for a Timoshenko beam element
% element with modulus of elasticity E,
% cross-sectional area A, moment of
% inertia I, length L, and angle
% theta (in degrees).
% The size of the element stiffness
% matrix is 6 x 6.
con=E*I/(1+phi);
w1 = 12*con/L^3+1.2*P/L;
w2 = 6*con/L^2+P/10;
w3 = con*(4+phi)/L+2*P*L/15;
w4 = con*(2-phi)/L-P*L/30;
y = [w1,w2,-w1,w2;w2,w3,-w2,w4;-w1,-w2,w1,-w2;w2,w4,-w2,w3];
function y = TimoshenkoElementMass(rho,A,I,l,phi)
%TimoshenkoElement Mass matrix This function returns the mass
% matrix for a Timoshenko beam
% element with mass density rho,
% cross-sectional area A, length L, and
% angle theta (in degrees).
% The size of the element stiffness
% matrix is 4 x 4.
a(1,1)=13/35+7*phi/10+phi^2/3;
a(1,2)=(11/210+11*phi/120+phi^2/24)*l;
a(1,3)=9/70+3*phi/10+phi^2/6;
a(1,4)=-(13/420+3*phi/40+phi^2/24)*l;
Finite element method in relation to structural dynamics 521
a(2,2)=(1/105+phi/60+phi^2/120)*l^2;
a(2,3)=(13/420+3*phi/40+phi^2/24)*l;
a(2,4)=-(1/140+phi/60+phi^2/120)*l^2;
a(3,3)=(13/35+7*phi/10+phi^2/3);
a(3,4)=-(11/210+11*phi/120+phi^2/24);
a(4,4)=(1/105+phi/60+phi^2/120)*l^2;
b(1,1)=1.2;
b(1,2)=(0.1-0.5*phi)*l;
b(1,3)=-1.2;
b(1,4)=(0.1-0.5*phi)*l;
b(2,2)=(2/15+phi/6+phi^2/3)*l^2;
b(2,3)=(-1/10+phi/2)*l;
b(2,4)=(-1/30-phi/6+phi^2/6)*l^2;
b(3,3)=1.2;
b(3,4)=(-0.1+0.5*phi)*l;
b(4,4)=(2/15+phi/6+phi^2/3)*l^2;
for i=2:4
for j=1:(i-1)
a(i,j)=a(j,i);
b(i,j)=b(j,i);
end
end
y=(rho*A*l)*a+(rho*I/((1+phi^2)*l))*b;
function y = TimoshenkoAssemble(K,k,i,j)
%Timoshenko beam This function assembles the element stiffness
% matrix k of the Timoshenko beam with nodes
% i and j into the global stiffness matrix K.
% This function returns the global stiffness
% matrix K after the element stiffness matrix
% k is assembled.
lm(1)=2*i-1;
lm(2)=2*i;
lm(3)=2*j-1;
lm(4)=2*j;
for l=1:4
ii=lm(l);
for n=1:4
jj=lm(n);
K(ii,jj)=K(ii,jj)+k(l,n);
end
end
y = K;
522 Structural dynamics of earthquake engineering
end
y = K;
OUTPUT
Natual frequencies in rad/sec
1.0e+002 *
0.01000000000000
0.11701648818926
0.23474939088448
0.35376511095837
0.47263579452189
0.61764685636320
0.77724091967848
0.93007476129177
1.12011725011258
1.36597976298605
1.66895296196556
1.66895296196556
