Develop a simulation of response of the slider-crank mechanism shown in the picture over a period of 2 seconds. The data for the problem is as follows:
Crank length = 0.2m, crank mass = 1kg.
Connecting rod length = 0.5m, connecting rod mass = 3kg.
Slider mass = 2kg.
Spring constant = 10000 N/m, unstretched length = 0.5m.
Viscous damping coefficient = 1000 Ns/m
Crank length = 0.2m, crank mass = 1kg.
Connecting rod length = 0.5m, connecting rod mass = 3kg.
Slider mass = 2kg.
Spring constant = 10000 N/m, unstretched length = 0.5m.
Viscous damping coefficient = 1000 Ns/m
Case 1: A
constant torque of 500 Nm (counter clockwise) is applied to the crank, starting
from rest.
The starting crank angle is 30 degrees above the horizontal
Constraint Equations:
For Joint P1
If we further manipulate the above equations and substitute
values of x1, y1 and x2, y2 from
equations (1), (2), (3), and (4) into equations (3), (4), (5) and (6), we get
the following constraint equations:
If we determine the generalized coordinates as qs
considering that x1=q1 , y1=q2, θ1=q3, x2=q3,
y2=q5, θ2=q6,
x3=q7, y3 = q8, and θ3=q9, we
can write the constraint equations as:
Constraint forces are:
     MATLAB Codes for Case 1 |
function zdot = Project_1(t,z)
% Function to evaluate the derivatives for the SliderCrank Mechanism
% Separate position and velocity information
q = z( 1: 9); % position data
qdot = z( 10: 18); % velocity data
% Acceleration of gravity
g = 9.81; % Units - m/sec.^2
% Input torgue
tau = 500; % Units - N m
%Spring and Damper Constants
k = 10000; % Spring constant Units N/m
c = 1000; % Dashpot constant Units - Ns/m
% Rigid Body Data
m1 = 1; % Units - Kg
m2 = 3; % Units - Kg
m3 = 2; % Units - Kg
L1 = 0.2; % Units - m
L2 = 0.5; % Units - m
I1 = (1/12)*m1*L1^2; % Units - Kg m^2
I2 = (1/12) * m2 * L2^2; % Units - Kg m^2
I3 = 0; % Units - Kg m^2
% Build the Mass matrix
M = zeros( 9, 9);
M( 1, 1) = m1;
M( 2, 2) = m1;
M( 3, 3) = I1;
M( 4, 4) = m2;
M( 5, 5) = m2;
M( 6, 6) = I2;
M( 7, 7) = m3;
M( 8, 8) = m3;
M( 9, 9) = I3;
% Specify the generalized forces
Q = zeros( 9,1);
Q( 2) = -m1*g;
Q( 3) = tau;
Q( 5) = -m2*g;
Q( 7) = k*(0.5-q(7))-c*qdot(7);
Q( 8) = -m3*g;
% Specify the derivative of the constraint matrix
Phiq = zeros( 8,9);
% Body 1
Phiq( 1, 1) = 1;
Phiq( 1, 3) = L1/2*sin(q(3));
Phiq( 2, 2) = 1;
Phiq( 2, 3) = -L1/2*cos(q(3));
Phiq( 3, 3) = L1*sin(q(3));
Phiq( 3, 4) = 1;
Phiq( 3, 6) = L2/2*sin(q(6));
%**********************************************************************
% Body 2
Phiq( 4, 3) = -L1*cos(q(3));
Phiq( 4, 5) = 1;
Phiq( 4, 6) = -L2/2*cos(q(6));
Phiq( 5, 3) = L1*sin(q(3));
Phiq( 5, 6) = L2*sin(q(6));
Phiq( 5, 7) = 1;
Phiq( 6, 3) = -L1*cos(q(3));
Phiq( 6, 6) = -L2*cos(q(6));
Phiq( 6, 8) = 1;
%**********************************************************************
%Body 3
Phiq( 7, 8) = 1;
Phiq( 8, 9) = 1;
%**********************************************************************
% Right Side of Constraint Equation
Phiqdot = zeros( 8,1);
Phiqdot( 1) = (L1/2)*cos(q(3))*qdot(3)^2;
Phiqdot( 2) = (L1/2)*sin(q(3))*qdot(3)^2;
Phiqdot( 3) = L1*cos(q(3))*qdot(3)^2+L2/2*cos(q(6))*qdot(6)^2;
Phiqdot( 4) = L1*sin(q(3))*qdot(3)^2+L2/2*sin(q(6))*qdot(6)^2;
Phiqdot( 5) = L1*cos(q(3))*qdot(3)^2+L2*cos(q(6))*qdot(6)^2;
Phiqdot( 6) = L1*sin(q(3))*qdot(3)^2+L2*sin(q(6))*qdot(6)^2;
% Build Coefficient Matrix
C = [M Phiq';Phiq zeros( 8,8)];
D = inv(C);
% Build Right Vector
R = [ Q' -Phiqdot']';
% FInd Solution
ACC = D*R;;
qdd = ACC(1: 9);
% Determine the time derivative of the state zector
zdot = [qdot' qdd']';
close all;
clear all;
% Class Example
% Set the parameters
Tspan = 0.0:0.0001:2;
g = 9.81;
L1 = 0.2;
L2 = 0.5;
theta1 = pi/6;
theta2 = 2.940235-pi;
% Determine the initial conditions
q = zeros(9,1);
q(1) = 0.2/2*cos(theta1);
q(2) = 0.2/2*sin(theta1);
q(3) = theta1;
q( 4) = 2*(.1*cos(theta1))+0.5/2*cos(theta2);
q( 5) = -0.5/2*sin(theta2);
q(6) = theta2;
q(7) = 2*(0.1*cos(theta1))+0.5*cos(theta2);
qdot = zeros(9,1); % All initial velocities are zero
Z0 = [q' qdot']'; % Initial Condition Vector
options = odeset('RelTol',1.0e-9,'AbsTol',1.0e-6);
[Tout,Z] = ode45(@Project_1,Tspan,Z0,options);
% Body 1
x1 = Z(:,1);
y1 = Z(:,2);
Theta1 = Z(:,3);
figure(1);
plot(Tout, x1, 'k', Tout, y1,'k--');
grid on;
title('Crank Position');
xlabel('Time - seconds');
ylabel('Position - meters');
legend('Horizontal Position','Vertical Position');
figure(2);
plot(Tout, Theta1, 'k');
grid;
title('Angle of Crank');
xlabel('Time - seconds');
ylabel('Angle - radians');
x1dot = Z(:,1+ 9);
y1dot = Z(:,2+ 9);
Theta1dot = Z(:,3+ 9);
figure(3);
plot(Tout, x1dot, 'k', Tout, y1dot, 'k--', Tout, Theta1dot,'k-.');
grid;
title('Velocities of Crank');
xlabel('Time - seconds');
ylabel('Velocities - meters/s and angular velocity - rad/sec.');
legend('Horizontal Velocity','Vertical Velocity', 'Angular Velocity');
% Connecting Rod
x2 = Z(:,4);
y2 = Z(:,5);
Theta2 = Z(:,6);
figure(4);
plot(Tout, x2, 'k', Tout, y2,'k--',Tout,Theta2,'k-.');
grid;
title('Connecting Rod Position');
xlabel('Time - seconds');
ylabel('Position - meters and angle - rad');
legend('Horizontal Position','Vertical Position','Angular Position');
x2dot = Z(:,4+ 9);
y2dot = Z(:,5+ 9);
Theta2dot = Z(:,6+ 9);
figure(5);
plot(Tout, x2dot, 'k', Tout, y2dot, 'k--', Tout, Theta2dot,'k-.');
grid;
title('Velocities of Connection Rod');
xlabel('Time - seconds');
ylabel('Velocities - meters/s and angular velocity - rad/sec.');
legend('Horizontal Velocity','Vertical Velocity', 'Angular Velocity');
% Slider
x3 = Z(:,7);
y3 = Z(:,8);
Theta3 = Z(:,9);
figure(6);
plot(Tout, x3, 'k', Tout, y3,'k--',Tout,Theta3,'k-.');
grid;
title('Slider Position');
xlabel('Time - seconds');
ylabel('Position - meters and angle - rad');
legend('Horizontal Position','Vertical Position','Angular Position');
x3dot = Z(:,7+ 9);
y3dot = Z(:,8+ 9);
Theta3dot = Z(:,9+ 9);
figure(7);
plot(Tout, x3dot, 'k', Tout, y3dot, 'k--', Tout, Theta3dot,'k-.');
grid;
title('Velocities of Slider');
xlabel('Time - seconds');
ylabel('Velocities - meters/s and angular velocity - rad/sec.');
legend('Horizontal Velocity','Vertical Velocity', 'Angular Velocity');
Results of Case 1
You can find Case 2 of this problem, where the crank is being driven at a constan angular speed of 60rpm clockwise, from here...
