Finite element method in relation to structural dynamics during earthquakes: Abstract: The basic procedure of the finite element method, with application to simple vibration problems, is presented in this chapter. The element stiffness, mass (both consistent and lumped mass) and forced vibration are derived for truss elements, shafts and beam elements. The transformation of the above matrices with respect to the local coordinate system is now transformed to the global system. The equations of motion of the complete system of finite element and the incorporation of boundary conditions are discussed. Relevant computer programs in MATHEMATICA and MATLAB are presented for truss, torsion and beam elements. Although techniques presented are applicable to two- or-three dimensional systems, only the one-dimensional element is treated in this pages.
Finite
element method in relation to structural dynamics during earthquakes
Abstract: The basic
procedure of the finite element method, with application to simple
vibration problems, is presented in this chapter. The element stiffness, mass
(both consistent and lumped mass) and forced vibration are derived for truss
elements, shafts and beam elements. The transformation of the above matrices
with respect to the local coordinate system is now transformed to the global
system. The equations of motion of the complete system of finite element and
the incorporation of boundary conditions are discussed. Relevant computer
programs in MATHEMATICA and MATLAB are presented for truss, torsion and beam
elements. Although techniques presented are applicable to two- or-three
dimensional systems, only the one-dimensional element is treated in this
pages.
Key words: discrete element natural
frequency, modes, Rayleigh Ritz method, boundary conditions.
Introduction
The finite element method
is a powerful numerical method that is used to provide approximations to
continuous systems. The disciplines in which the finite element method can be
applied include stress analysis, heat transfer, electromagnetic analysis, fluid
flow and vibrations. Application of the finite element method to continuous
systems requires the systems to be divided into a finite number of discrete
elements. Interpolations for the dependent variables are assumed across
each element and are chosen to ensure appropriate inter-element continuity.
The interpolating functions are developed in terms of unknown values of the
dependent variables at discrete points, called nodes which are located for a
one-dimensional system at element boundaries. The defined interpolations are
used to provide approximation to the dependent variables across the system.
Lagrange’s equations are then applied for vibration problems resulting in a set
of differential equations for the dependent variables at nodes. Total structure
is obtained as assembly of elements.
The term finite element was coined by Prof Clough. Boundary
conditions for continuous systems are classified as being of two types:
geometric boundary conditions are those which must be
satisfied according to geometric constraints. For example, u(x =
0) = 0 u′(x
= 0) = 0 at a fixed end of a cantilever;
• natural
boundary conditions are those that must be satisfied as a result of free and
moment balances. For example, EI (∂^{2}u/∂x^{2}) (x = L) = 0 is a
moment boundary condition at the free end of a cantilever.
The chosen interpolating function
must satisfy geometric boundary conditions.
The type of problem that comes
from the realm of structural dynamics is restricted to the calculation of natural
frequencies and the corresponding mode shapes of free undamped vibration of
common structural components and forms. This requires the development of mass
matrix, which will represent the effect of dynamic loading (proportional to the
square of frequency) which is set up during vibration. Vibration problems are
eigenvalue problems in which eigenvalues represent the square of the natural
frequencies and eigenvectors define the shape of the structure when vibrating
at a particular frequency.
Dynamic analysis
Basically four different types of
problems can be distinguished in the field of dynamics, free vibration, steady
state vibration, transient response to known excitations and response to random
excitations.
Torsional vibration of a
shaft
The shaft element shown in Fig. 14.1 has two nodes at its two
ends. The unknown displacements at each end are the angles of twist φ_{1}, φ_{2}. The
displacement function, which is the angle of twist, is given by
Let us consider undamped free vibration case. The inertia
force is Iφ^{˙˙} acting
on the shaft where I is the mass moment of inertia/unit length. Hence total
potential energy can be written as
If [M]_{e} is a lumped mass matrix (it
is assumed that mass is lumped at nodes)
Even though consistent mass is
more accurate, lumped mass gives better results because both stiffness and mass
are over-estimated, thus resulting in the correct answer.
The next step is the assembly of stiffness matrix. Assume that the shaft is idealized into a number of elements as shown in Fig. 14.2 and that each element has two nodes and one degree of freedom at each node. The elements of the stiffness matrix will go into the corresponding position of global stiffness matrix as
When all elements of the element stiffness matrices and mass
matrices are assembled, and the boundary conditions incorporated, the final
equations of free vibration is as follows
Equation 14.26 is a typical eigenvalue problem.
Example 14.1
A uniform shaft having one end
fixed and the other end free is divided into two equal elements as shown in
Fig. 14.3. Determine the natural frequency for torsional vibration of the
shaft.
Axial vibration of rods
The total degrees of freedom for a bar element are the axial
displacements at the ends of the element instead of the angle of twist for
torsional vibration. Let u_{1}, u_{2} be the
displacements at the left end and right end of an element (see Fig. 14.1).
Axial displacement at any section is written as
Hence the stiffness matrix of an element for axial vibration
of a rod is the same as torsional vibration of the rod except that GJ is
to be replaced by EA.
Hence I in the torsional
vibration must be replaced by ρA for axial vibration of the rod.
Example 14.2
For the rod shown in Fig 14.3 determine the lowest frequency
for longitudinal vibration of the rod by the finite element method by
considering two elements. Solution
In the above method, admissible
functions are used as basis functions in a Rayleigh�'Ritz approximation to
solution of an eigenvalue problem. Sometimes the Rayleigh�'Ritz method is
difficult to apply to vibration problems. The assumed modes method, introduced
in the next section, is based on application of Lagrange’s equations and leads
to the same approximation as that for the same set of interpolating functions
as the Rayleigh�'Ritz method. In the next section we will use the finite element
method for the longitudinal vibration by the assumed modes method.
14.5 Assumed
modes method
Example 14.3
Consider a forced vibration of a
longitudinal bar as shown in Fig. 14.5. The displacement u is a function
of spatial coordinate x and time t, i.e. u(x, t).
Program
14.1: MATLAB program for the assumed modes method
%example 14.3
% assumed modes method to
determine natural frequencies m mode %shapes, and forced response of tapered
bar with attached mass and spring clc;
close all; digits(5); x=sym(‘x’);
%parameters
e=200*10^9; % youngs modulus in
N/sq.m rho=7800; % mass density in kgm/cu.m l=4.0; % span of the beam in m
m=10; % concentrated mass in kgm
k=3*10^7; % spring constant %functions
a=0.001*(1-.005*x);% area in sq.m
u(1)=sin(pi*x/(2*l)); % assumed modes u(2)=sin(3*pi*x/(2*l)); u(3)=sin(5*pi*x/(2*l));
%mass and stiffness matrices for
i=1:3
for j=1:i c1=subs(u(i),x,l)*subs(u(j),x,l);
Mint=rho*a*u(i)*u(j);
Kint=e*a*diff(u(i),x)*diff(u(j),x);
M(i,j)=int(Mint,x,0,l)+m*c1;
K(i,j)=int(Kint,x,0,l)+k*c1;
M(j,i)=M(i,j);
K(j,i)=K(i,j); end
end
disp(‘ stiffness matrix’)
K=vpa(K)
disp(‘ mass matrix’) M=vpa(M)
K1=double(K); M1=double(M); C=inv(K1)*M1; [V,D]=eig(C);
for i=1:3 w(i)=1/sqrt(D(i,i));
end
disp(‘ natural frequencies in
rad/sec’) w=vpa(w);disp(w’)
%Normalize mode shape vectors
E=V’*M*V;
for j=1:3 for i=1:3
P(i,j)=V(i,j)/sqrt(E(j,j)); end
end
disp(‘ modal matrix’)
P=vpa(P);disp(P) %mode shapes xx=linspace(0,l,37); P1=single(P);
for k=1:37 for i=1:3
v(i,k)=0; for j=1:3
v(i,k)=v(i,k)+P1(j,i)*subs(u(j),x,xx(k)); end
end end
plot(xx,v(1,:),’k’,xx,v(2,:),’*k’,xx,v(3,:),’�'k’);
xlabel(‘ x (m)’)
ylabel(‘ w(x) m’)
legend(‘ 1st mode’, ‘2nd mode’, ‘3rd mode’)
Output for example 14.3 stiffness
matrix
K =
[ .91319e8, -.29250e8, .30139e8]
[ -.29250e8, .57987e9, -.26250e8] [ .30139e8, -.26250e8, .15570e10]
mass matrix M =
[ 25.381, -9.9368, 9.9930] [
-9.9368, 25.437, -9.9368] [ 9.9930, -9.9368, 25.442]
natural frequencies in rad/sec
1895.7 5062.1 8984.0
modal matrix
[ .19713, -.94408e-1, -.53048e-1]
[ -.26221e-2, -.20576, .89500e-1] [ .79332e-3, .28353e-1, .22280]
Truss element
Consider a truss element oriented as shown in Fig. 14.7 in the
global coordinate system.
14.7.1 Element
stiffness and mass matrices
Since the truss element is a bar element subjected to axial
forces, the stiffness matrix is given by
The local displacements <w_{i} w_{j}>
can be written in terms of global displacements as
where {q}_{l} represents local
displacements and {q}_{g} represents global
displacements.
By using contra-gradient law
{F}_{g} = [T]^{T}{F}_{l --- -- }14.50_{}
{F}_{g} = [T]^{T}[k]_{l}{q}_{l
}_{ --- -- }14.51
_{ }
14.7.2 Assembly
The elements of the stiffness matrix and mass matrix will
assemble in the proper location of structural stiffness matrix as
3 Application of
boundary conditions
• If any of
the degrees of freedom is constrained, the row and column corresponding to that
degree of freedom are deleted from the assembled stiffness or mass matrix.
• Add
springs of very high stiffness at the constrained degree of freedom.
• Use
Lagrangian multiplier method to incorporate the constraints.
4 Solve as an eigenvalue problem_{}
and the corresponding mode shape is also obtained.
Example 14.4
Use the finite element method to obtain the lowest natural
frequency for the
truss shown in Fig. 14.8. The data are: L_{1} =
1.2 m; L_{2} = 2.68 m; L_{3} = 2.4 m; L_{4}
= 1.2 m; θ = 63.43� o ; sin θ = 0.894; cos θ = 0.447
E = 2 �- 10^{11}
N/m^{2}; A = 0.04 m^{2}; ρ = 7600 kg/m^{3}.
Assembling all the elements and eliminating the 1st, 2nd, 7th
and 8th (boundary conditions) we get resulting stiffness and mass matrices as
Program
14.2: MATLAB program for free vibration of trusses
Trussdyn
% solution of plane
truss free vibration by finite element method
clc;
K=zeros(12,12);
M=zeros(10,10);
e=200e9;
a=0.04;
l1=1.2;
rho=7600;
l2=2.683;
l3=2.4;
l4=1.2;
l5=2.4;
% calculate element
lengths
% calculate element
stiffness
k1=PlaneTrussElementStiffness(e,a,l1,0);
m1=PlaneTrussElementMass(rho,a,l1,0);
k2=PlaneTrussElementStiffness(e,a,l2,63.43);
m2=PlaneTrussElementMass(rho,a,l2,63.43);
k3=PlaneTrussElementStiffness(e,a,l3,90);
m3=PlaneTrussElementMass(rho,a,l3,90);
k4=PlaneTrussElementStiffness(e,a,l4,0);
m4=PlaneTrussElementMass(rho,a,l4,0);
% assemble element
stiffness to global stiffness
K=PlaneTrussAssemble(K,k1,4,3);
K=PlaneTrussAssemble(K,k2,1,3);
K=PlaneTrussAssemble(K,k3,2,3);
K=PlaneTrussAssemble(K,k4,1,2);
M=PlaneTrussAssemble(M,m1,4,3);
M=PlaneTrussAssemble(M,m2,1,3);
M=PlaneTrussAssemble(M,m3,2,3);
M=PlaneTrussAssemble(M,m4,1,2);
format long;
ks=1;
% apply boundary conditions for
stiffness matrix and mass K(9,1)=ks;
K(10,2)=ks;
K(11,7)=ks;
K(12,8)=ks;
M(9,1)=ks;
M(10,2)=ks;
M(11,7)=ks;
M(12,8)=ks;
K(1,9)=ks;
K(2,10)=ks;
K(7,11)=ks;
K(8,12)=ks;
M(1,9)=ks;
M(2,10)=ks;
M(7,11)=ks;
M(8,12)=ks;
invk=inv(K);
km=invk*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 k=j;
Qh=ee(j); else
end end Ql=Qh;
Qh=max(ee)+0.001;
om1(i)=ee(k);
omega(i)=sqrt(ee(k)); for m=1:ms
p1(m,i)=evec(m,k); end
end
%Normalizing the mode shape
L=p1'*m*p1;
%develop modal matrix for i=1:ms
for j=1:ms p(i,j)=p1(i,j)/L(j,j);
end end
disp(‘ Natural frequencies in
rad/sec’) disp(omega’)
disp(‘ normalized modal vector ‘)
disp(p)
function y =
PlaneTrussElementStiffness(E,A,L, theta) %PlaneTrussElementStiffness This
function returns the element
%stiffness
matrix for a plane truss
%element
with modulus of elasticity E,
%cross-sectional
area A, length L, and
%angle
theta (in degrees).
%The size
of the element stiffness
%matrix is
4 x 4.
x = theta*pi/180; C = cos(x);
S = sin(x);
y =
E*A/L*[C*C C*S -C*C -C*S ; C*S S*S -C*S -S*S ; -C*C -C*S C*C C*S ; -C*S -S*S
C*S S*S];
function y =
PlaneTrussElementMass(rho,A,L, theta) %PlaneTrussElementStiffness This function
returns the mass
%matrix
for a plane truss
%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.
x = theta*pi/180; C = cos(x);
S = sin(x);
% for consistent mass use the following
y = (rho*A*L/6)*[2 0 1 0 ; 0 2 0
1 ;1 0 2 0 ; 0 1 0 2]; %for lumped mass use the following %y=(rho*A*L/2)*[1 0 0
0;0 1 0 0;0 0 1 0;0 0 0 1];
function y = PlaneTrussAssemble(K,k,i,j)
%PlaneTrussAssemble This function assembles the element
stiffness
%matrix k
of the plane truss element 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 m=1:4 ii=lm(m); for n=1:4
jj=lm(n);
K(ii,jj)=K(ii,jj)+k(m,n); end
end
y = K;
OUTPUT
Natural frequencies in rad/sec 1.0e+003 *
0.00100000000000
0.00100000000000
1.25105969735243
3.23474518944287
4.53189428385080
Beam element
Consider a uniform beam element of length L and
cross-sectional area A and mass density ρ as shown in Fig. 14.9. The
modulus of elasticity of the material is E and I is the second
moment of area. The unknown displacement of the element are deflections and
rotations at the two ends, in total four degrees of freedom for each element or
two degrees of freedom/node. The displacement function is represented by the
equation given by
Program 14.3: MATHEMATICA program for
evaluation of stiffness matrix, and mass matrix of a beam element
Example 14.5
Find the fundamental frequency of
a simply supported uniform beam shown in Fig. 14.10.
Example 14.6
Set up the system of equation
governing free vibration in its own plane of the portal frame shown in Fig.
14.11.
A = 1.851
87 �- 10^{�'5};
I = 2.857 85 �- 10^{�'11};
E = 210 GPa;
ρ = 25 613.5 kg/m^{3}; L
= 0.2413 m.
Program 14.4: MATLAB program
to find the natural frequency of beams or rigid frames
FRAMEDYN
% dynamics of plane frame by
finite element method clc;
nj=7;
ne=6;
neq=3*nj;
K=zeros(neq,neq);
M=zeros(neq,neq);
%give nodi and nodj of each member
nodi=[1 2 3 4 5 6];
nodj=[2 3 4 5 6 7];
%give the values of e,a,i angle
and lengths of members e=210e9;
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=[90 90 0 0
-90 -90];
l=[.12065 .12065 .12065 .12065
.12065 .12065];
%give
density of the material
rho=25613.5;
%number of constraint degrees of
freedom nbou=6;
%the numbers of constrained
degrees of freedom nb=[1 2 3 19 20 21];
%
form 6 x 6 element stiffness and mass matrix and
assemble wilson method for n=1:ne
i=nodi(n);
j=nodj(n);
k=PlaneFrameElementStiffness(e,a(n),mi(n),l(n),angle(n));
m=PlaneFrameElementMass(rho,a(n),l(n),angle(n));
K=PlaneFrameAssemble(K,k,i,j);
M=PlaneFrameAssemble(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 k=j;
Qh=ee(j); else
end end Ql=Qh;
Qh=max(ee)+0.001;
om1(i)=ee(k);
omega(i)=sqrt(ee(k)); for m=1:ms
p1(m,i)=evec(m,k); end
end
%Normalizing the mode shape
L=p1'*m*p1;
%develop modal matrix for i=1:ms
for j=1:ms p(i,j)=p1(i,j)/L(j,j); end
end
disp(‘ Natural frequencies in
rad/sec’) disp(omega’)
disp(‘ normalized modal vector ’)
disp(p)
function y =
PlaneFrameElementStiffness(E,A,I,L,theta) %PlaneFrameElementStiffness This
function returns the element
%stiffness
matrix for a plane frame
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.
x = theta*pi/180; C = cos(x);
S = sin(x);
w1 = A*C*C + 12*I*S*S/(L*L);
w2 = A*S*S + 12*I*C*C/(L*L);
w3 = (A-12*I/(L*L))*C*S;
w4 = 6*I*S/L;
w5 = 6*I*C/L;
y = E/L*[w1
w3 -w4 -w1 -w3 -w4 ; w3 w2 w5 -w3 -w2 w5 ; -w4 w5 4*I w4 -w5 2*I ; -w1 -w3 w4
w1 w3 w4 ;
-w3 -w2 -w5 w3 w2 -w5 ; -w4 w5 2*I w4 -w5 4*I];
function y =
PlaneFrameElementMass(rho,A,l, theta) %PlaneFrameElementMass This function
returns the mass
%matrix
for a plane frame
%element
with mass density rho,
%cross-sectional
area A, length L, and
%angle
theta (in degrees).
%The size
of the element stiffness
%matrix is
6 x 6.
x = theta*pi/180; C = cos(x);
S = sin(x);
% for consistent mass use the following
%mass matrix of frame element consistent matrix
f33 _ 00=[0 . 333,0,0,0 .
167,0,0;0,0 . 37143,0 . 05238*l,0,0 . 12857, -
0.03095*l;0,.05238*l,.00952*l*l,0,0.03095*l,-.00714*l*l;
0.167,0,0,0.333,0,0;0,0.12857,0.03095*l,0,0.37143,-0.05238*l;0,-0.03095*l,-0.00714*l*l,0,-0.05238*l,0.00952*l*l];
t=[C,S,0,0,0,0;-S,C,0,0,0,0;0,0,1,0,0,0;0,0,0,C,S,0;0,0,0,-S,C,0;0,0,0,0,0,1];
n=t’*f33_00*t;
%lumped mass
%n=[0.5,0,0,0,0,0;0,.5,0,0,0,0;0,0,0,0,0,0;0,0,0,.5,0,0;0,0,0,0,.5,0;0,0,0,0,0,0];
y=rho*A*l*n;
function y =
PlaneFrameAssemble(K,k,i,j)
%PlaneFrameAssemble This function
assembles the element stiffness
%matrix k
of the plane frame element 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)=3*i-2; lm(2)=3*i-1;
lm(3)=3*i; lm(4)=3*j-2; lm(5)=3*j-1; lm(6)=3*j; for l=1:6
ii=lm(l); for n=1:6
jj=lm(n);
K(ii,jj)=K(ii,jj)+k(l,n); end
end
y = K;
OUTPUT
Natual frequencies in rad/sec 1.0e+004 *
0.00010000000000
0.01957884271941
0.07771746384284
0.12745125003361
0.13873461021422
0.31352858884033
Forced vibration of a beam
Example 14.7
Use a two element model for the
beam to determine the steady state response of the system shown in Fig. 14.13.
clc;
digits(5);
L=8; % length in m
rho=7600; %mass density in
kg/cu.m e=2e11; %youngs modulus of the material i=1.6*10^-6; % moment of
inertia in m^4 a=3.6*10^-3 % area in m^2
m=20; % mass of hanging block in
kg k=3*10^4; % stiffness of discrete spring in N/m s=L/2; %element length
f0=2500; % excitation amplitude
in N om=80; %excitation frequency in rad/s disp(‘global mass matrix’)
M=rho*a*s/420*[4*s^2,13*s,-3*s^2,0,0;13*s,312,0,-13*s,0;...
-3*s^2,0,8*s^2,-3*s^2,0;0,-13*s,-3*s^2,4*s^2,0;...
0,0,0,0,420*m/(rho*a*s)]; K=e*i/s^3*[4*s^2,-6*s,2*s^2,0,0;-6*s,24+k*s^3/(e*i),0,6*s,-k*s^3/(e*i);...
2*s^2,0,8*s^2,2*s^2,0;0,6*s,2*s^2,4*s^2,0;...
0,-k*s^3/(e*i),0,0,k*s^3/(e*i)];
M1=vpa(M);disp(M1) K1=vpa(K);disp(K1)
%natural frequencies
W2=eig(inv(K)*M); for i=1:5
w(i)=1/sqrt(W2(i)); end
disp(‘ natural frequencies’)
w=vpa(w’)
%force vector disp(‘ force
vector’)
f=f0*s*[-s/48;13/16;-s/8;-5*s/192;0]
%use of undetermined coefficients z=-om^2*M+K;
W=inv(z)*f;
x=linspace(0,L,21); for k=1:21
if x(k) <s xi=x(k)/s;
y(k)=(xi-2*xi^2+xi^3)*W(1)+(3*xi^2-2*xi^3)*W(2);
y(k)=y(k)+(-xi^2+xi^3)*W(4);
else xi=(x(k)-s)/s;
y(k)=(1-3*xi^2+2*xi^3)*W(2)+(xi^2-2*xi^2+xi^3)*W(3);
y(k)=y(k)+(-xi^2+xi^3)*W(4); end
end plot(x,y,‘-k’); xlabel(‘
x(m)’);
ylabel(‘w(x) (m)’);
W=vpa(W);
disp(‘ steady state amplitude in m’);disp(W)
OUTPUT
global mass matrix
[ 16.677, 13.550, -12.507, 0.,
0.] [ 13.550, 81.298, 0., -13.550, 0.] [ -12.507, 0., 33.353, -12.507, 0.] [
0., -13.550, -12.507, 16.677, 0.] [ 0., 0., 0., 0., 20.]
[ .32000e6, -.12000e6, .16000e6, 0., 0.]
[ -.12000e6, .15000e6, 0.,
.12000e6, -30000.] [ .16000e6, 0., .64000e6, .16000e6, 0.]
[ 0., .12000e6, .16000e6,
.32000e6, 0.] [ 0., -30000., 0., 0., 30000.]
natural frequencies w =
15.162
42.623
74.044
186.76
339.31
force vector f =
1.0e+003
*
-0.8333 8.1250 -5.0000 -1.0417 0
steady state amplitude in m -.43073e-1
-.10930e-1
.25384e-1 -.22861e-1
.33458e-2
Vibration of a Timoshenko beam
We have already seen that
Euler�'Bernoulli’s theory predicts the frequencies for a shallow beam with
adequate precision. But with the increasing depth of the beam, the effect of
transverse shear deformation and rotary inertia become more important. Many
varieties of Timoshenko beam elements have been proposed. It has been observed
that the element described below is adequate for practical use.
The deflection function is given by
where k is Timoshenko’s
shear constant = 5/6 for rectangular section and 9/10 for circular sections.
Example 14.8
Find the fundamental frequency of
a Timoshenko beam by using the program ‘Timoshenkovib’ given that
E = 210
GPa; G = 70 GPa; ρ = 7800 kg/m^{3}, A =
1.85187e�'5;
I =
2.85785e�'11; P = 200 N; k = 5/6.
The lowest frequency obtained from the program is 11.7 rad/s.
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));
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);
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;
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;
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
Summary
In this chapter, free and forced vibrations of beams, trusses
and frames are discussed. In the next chapter, we will apply other numerical
methods such as differential quadrature and transformation methods to find the
natural frequencies of structures.