Topic: program MATLAB

CASE 1
Problem Statement:
Damped free vibrations can be modelled by considering a block of mass m that is
attached to a spring and a dashpot as shown.

From Newton’s second law of motion, the displacement x of the mass as a function of
time can be determined by solving the differential equation:

where k is the spring constant, and c is the damping coefficient of the dashpot. If the mass
is displaced from its equilibrium position and then released, it will start oscillating back
and forth. The nature of the oscillations depends on the size of the mass and the values of
k and c.
a)    MATLAB Script file:
clear all;
clf;
%This is a Script file that calls for m,k,c and x
m=input(‘Please enter the mass of the Block(m) ’);            %expects the mass of the Block from the user
k=input(‘Please enter the spring constant(k) ’);              %expects the spring constant from the user
c=input(‘Please enter the damping coefficient(c) ’);          %expects the damping coefficient from the user
T0=input(‘Please enter the starting t ’);                     %expects the starting time from the user
Tn=input(‘Please enter the ending t ’);                       %expects the ending time from the user
x(1,1)=input(‘Please enter the x1 value at the starting t ’); %expects the initial value of x
x(2,1)=input(‘Please enter the x2 value at the starting t ’); %expects the velocity @ t=0
fname=input(‘Please enter the function file please:’,’s’);
[tout,xout]=ode45(fname,[T0 Tn],[x(1,1) x(2,1)]);
fprintf(‘The Displacement is %f at time %f ’,xout,Tn);
fprintf(‘The Velocity is %f at time %f ’,tout,Tn);

b)    M-File to calculate x and v
Since the equation in the question is a Second Degree Linear Differential equation we will use the ODE45 MATLAB code. The

function can be called by using the following equation:
[tout,xout]=ode45(fname,[T0 Tn],[x(1,1) x(2,1)]);

tout = Velocity of the Block at a certain time
xout = Displacement of the Block at a certain time
The Vanderpol mentioned in Chapter 13 is used to solve the Second Degree Linear Differential equation into 2 First Degree

Differential equations.
% this funtion uses the vanderpol function
function yp = vanderpol(t,x)
yp=[x(2);(-0.1)*(3*x(2)+(28*x(1)))];

c)    M-File to create the animation of the mass of motion
The “movie” function is used to create an animation for the Block.
L=4;           %length of the Block
W=2;           %width of the Block
n=length(xout) %no. time we get xout
for i=1:n
axis equal;
axis([-3 18 -3 18]);
hold on;
X(1) = xout(i);
Y(1) = xout(i);
X(2) = X(1) + L;
Y(2) = xout(i);
X(3) = X(2)+W;
Y(3) = Y(2)+W;
X(4) = X(1)+W;
Y(4) = Y(1)+W;
X(5) = X(1);
Y(5) = Y(1);
plot(X,Y);
M(i)=getframe;
clf;
end
h = picture;
movie(h,M,3,n/2,[100 50 0 0]);

d)    Demonstration of the 2 cases
I.&nbsp;&nbsp; &nbsp;C = 3N-s/m for 0<=t<=20s

II.&nbsp;&nbsp; &nbsp;C=50N-s/m for 0<=t<=10s

e)&nbsp;&nbsp; &nbsp;Description of the Program
A script file was generated which would be used to calculate the velocity and displacement of the Block. However the question

is a second order differential equation which can be simplified to two first degree differential equations. The “movie”

command was used to generate the animation of the oscillation of the block. In Part D, there are 2 scenarios. In the first

case damping coefficient is 3N-s/m and in case 2 the damping coefficient is 50N-s/m. Since the damping coefficient in case to

is greater than case 1 hence the block will oscillate slower as compared to case 1.
CASE 2
Since I was a FEM student last semester I would like to consider the FEM analysis of 2 beam elements.

L1=12in&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;L2=18in
E=30x106psi&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;Area=0.5in x 0.5in
?=45°&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;downward vertical fore @ node 2=1000lb
This problem can be solved with the MATLAB. The geometric parameters, material properties, and applied forces are input.
The first step in the solution procedure is to discretize the domain, i.e. select the number of elements. The element

stiffness matrix for each of the beam elements is of the form.

The second step is to transform each element stiffness matrix in local coordinates to the global coordinate system. This

transformation is accomplished by:

The transformation angle for element 1 is q = 0, and, thus, the transformation matrix is simply the identity matrix. Because

the rotation angle is measured counter-clockwise, the transformation angle for element 2 is q = -45°.
The third step is to assemble the global stiffness matrix that describes the entire structure by properly combining the

individual element stiffness matrices.
The fourth step is to apply the constraints and reduce the global stiffness matrix so that the specific problem of interest

can be solved. For this problem, the displacements at node 1 and node 3 are known, i.e. U1 = U2 = U3 = U7 = U8 = U9 = 0, and

the loads are applied to node 2 are specified, i.e. F4 = 0, F5 = -1000 lb, and F6 = 0. Thus, the system of equations can be

written as:

This problem has three unknowns, U4, U5, and U6, and, thus, requires three independent equations. Matrix algebra allows three

independent equations to be constructed by the removal of the rows and columns that correspond to the known displacements,

i.e. the global stiffness matrix reduces to:

The reduced Global stiffness matrix is as follows:

MATLAB program to solve this:
function beam2
%
% Step 0: Input constants
%
m = 2; % number of elements
b = [0.5 0.5]; % width (in)
h = [0.5 0.5]; % height (in)
a = b.*h; % cross-sectional area
I = b.*h.^3./12; % moment of inertia
l = [12 18]; % length of each element (in)
e = [30*10^6 30*10^6]; % modulus of elasticity (psi)
theta = [0 -45*pi/180]; % orientation angle
f = -1000; % force (lbs)
%
% Step 1: Construct element stiffness matrix
%
k = zeros(6,6,m);
for n = 1:m
k11 = a(n)*e(n)/l(n); k22 = 12*e(n)*I(n)/l(n)^3;
k23 = 6*e(n)*I(n)/l(n)^2; k33= 4*e(n)*I(n)/l(n);
k36 = 2*e(n)*I(n)/l(n);
k(:,:,n) = [k11 0 0 -k11 0 0;
0 k22 k23 0 -k22 k23;
0 k23 k33 0 -k23 k36;
-k11 0 0 k11 0 0;
0 -k22 -k23 0 k22 -k23;
0 k23 k36 0 -k23 k33];
end
%
% Step 2: Transform element stiffness matrices to global coordinates
%
shift = 0;
kbar = zeros(6,6,m); Ke = zeros(9,9,m);
for n = 1:m
c = cos(theta(n)); s = sin(theta(n));
lamda = [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];
kbar (:,:,n) = lamda’*k(:,:,n)*lamda;
%
% Step 3: Combine element stiffness matrices to form
global stiffness matrix
%
for i = 1:6
for j = 1:6
Ke(i+shift,j+shift,n) = kbar(i,j,n);
end
end
shift = shift + 3;
end
K = sum(Ke,3); Kr = K;
%
% Step 4: Reduce global stiffness matrix with constraints
%
Kr(:,9) = []; Kr(9,:) = []; Kr(:,8) = []; Kr(8,:) = [];
Kr(:,7) = []; Kr(7,:) = [];
Kr(:,3) = []; Kr(3,:) = []; Kr(:,2) = []; Kr(2,:) = [];
Kr(:,1) = []; Kr(1,:) = []
%
% Step 5: Solve for unknown displacements
%
Fr = [0; f; 0];
Ur = KrFr
%
% Step 6: Solve for forces
%
U = [0;0;0; Ur; 0;0;0];
F = K*U

The result obtained from MATLAB is U4 = -0.0016 in, U5 = -0.0064 in, and U6 = -0.0003 radians.

GET THIS ASSIGNMENT DONE FOR FREE NOW, JUST PLACE YOUR ORDER AND DO NOT MAKE ANY PAYMENT: WE WILL COMPLETE THE PAPER FIRST BEFORE YOU PAY AND YOU PAY ONLY AFTER YOU RECEIVE AND APPROVE YOUR COMPLETED PAPER! IF YOU ARE NOT IMPRESSED, WE WILL NOT ASK ANYTHING FROM YOU! WHAT ELSE WE WILL DELIVER IT WITHIN YOUR SPECIFIED TIME FROM 3 HOURS TO 10 DAYS DEADLINE.

&nbsp;

THIS IS ONLY A LIMITED OFFER! TAKE ADVANTAGE!!