Homework 4: Dynamic problems
Given:
Geometrical parameters: L = 1 m, bh = 1010 mm2 (fig. 1).
Material: Steel E = 200 GPa, = 7850 kg/m3
Number of elements: Choose 10 or 20?
b
h
L = 1m
Fig. 1. Calculation domain
Require:
1. Calculate consistent mass matrix.
2. Calculate the first five natural frequencies and their mode shapes.
Solution:
1. Calculate consistent mass matrix
The consistent mass matrix:
x2
M N
T
Ndx (1.1)
x1
Where:
Material density, = 7850 kg/m3
K Ke (1.2)
N Shape function
1 me 2
x1 x2
= -1 =1
e
Fig. 2. 1-D element
For the 2-node prismatic bar element moving along x, the stiffness shape functions are:
1 1
N (1.3)
2 2
Natural coordinate of element,
And the consistent mass matrix:
Ale 2 1
me 1 (1.4)
6 2
A Area of cross section
A bh (1.5)
b, h Cross-sectional dimensions
Length of element.
2. Calculating the first five natural frequencies
Dynamic equilibrium equation for the case of no damping:
Mu Ku F (1.6)
Where:
M Global mass matrix
n
M mie (1.7)
i1
K Global stiffness matrix
n
K k ie (1.8)
i1
kei Element stiffness matrix
AE 1 1
ke
le 1 1
(1.9)
F Global force matrix
n
F fie (1.10)
i1
fe Element force matrix:
fle 1
fe 1 (1.11)
2
3. Results
Applying the FEM theory, a MATLAB programs were established to solve the problem. The achieved results are
represented in following figures. The MATLAB codes are listed in the appendices.
Calculating parameters: n = 10 element, dtmax = 0.1981, dt = 0.1 < dtmax
Fig. 3. First five natural frequencies and their modal shapes
MATLAB code
%% Homework 4: Dynamic problem
clear all; clc;
E=200e9; % Elastic modulus
L=1000; % length of bar
b=10e-3; % width
h=10e-3; % height
A=b*h; % cross area
n=10; % number of elements
le=L/n; % length of each element
ro=7850; % mass density
M=zeros(n+1,n+1); % mass
K=zeros(n+1,n+1); % stiffness
edof=zeros(n,1); % element DOF
u0=zeros(n,1);
du0=zeros(n,1);
%% Time step Force
tfinal=5e-5; % final time to simulate
dtmax=L/sqrt(E/ro); % time increment (s)
dt=input('Choose the time step dt= ');
t=0:dt:tfinal;
nstep=size(t,2)-1; % number of steps
F=zeros(n+1,nstep+1); % force
F(1,:)=[];
%% Assebly mass, stiffnes matrice and load vector
me=ro*A*le/6*[2 1;1 2];
ke=A*E/le*[1 -1;-1 1];
for i=1: n
edof=[i i+1];
M(edof,edof)=M(edof,edof)+me;
K(edof,edof)=K(edof,edof)+ke;
end
MM=M; KK=K;
MM(1,:)=[];
MM(:,1)=[];
KK(1,:)=[];
KK(:,1)=[];
%% determine natural frequencies and modes
[V,lamda]=eig(KK,MM);
vv=zeros(1,n);
V=[vv; V];
omega=sqrt(lamda);
%% Print out the natural frequencies and modes
figure
n=size(V,2);
for i=1:n
V(:,i)=V(:,i)/V(size(V,1),i);
subplot(n/2,1,i)
plot(V(:,i))
ylabel('d')
title(['mode',num2str(i)])
hold on
grid on
end
