### Slider-Crank Mechanism-Case 2: A crank is being driven at a constant angular speed

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

﻿

Case 2: The crank is being driven at a constant angular speed of 60 rpm clockwise.
The starting crank angle is 30 degrees above the horizontal

#### Constaint Equations:

As derived in Case 1
In case 2, we need to find all initial values for the velocities by differentiating the constraint equation,

The last equation is added to case 2 because of the applied constant initial angular velocity to the crank.

﻿
The following m-file is used to find the initial values of velocities and accelerations of all of the three bodies.
``````

clear all;

L1=0.2;
L2=0.5;
g=9.81;
w=-2*pi

q= zeros(9,1);              % Initial Positions

q(3)=pi/6;
q(6) =-asin(0.2);
q(1)=L1/2*cos(q(3));
q(2) =L1/2*sin(q(3));
q(4) =L1*cos(q(3))+L2/2*cos(q(6));
q(5) =L1*sin(q(3))+L2/2*sin(q(6));
q(7) =L1*cos(q(3))+L2*cos(q(6));
q(8) =0;
q(9) =0

qdotc  =zeros(9,9); % The coefficient matrix of initial velocities

qdotc(1,1)=1;
qdotc(1,3)=q(2);%L1/2*sin(q(3));
qdotc(2,2)=1;
qdotc(2,3)=-q(1); %L1/2*cos(q(3));
qdotc(3,4)=1;
qdotc(3,3)=2*q(2);%L1*sin(q(3));
qdotc(3,6)=L2/2*sin(q(6));
qdotc(4,5)=1;
qdotc(4,3)=-2*q(1); %L1*cos(q(3));
qdotc(4,6)=-L2/2*cos(q(6));
qdotc(5,7)=1;
qdotc(5,3)=2*q(2); %L1*sin(q(3));
qdotc(5,6)=L2*sin(q(6));
qdotc(6,8)=1;
qdotc(6,3)=-2*q(1); %L1*cos(q(3));
qdotc(6,6)=-L2*cos(q(6));
qdotc(7,8)=1;
qdotc(8,9) =1;
qdotc(9,3)=1;

qdotR=zeros(9,1); % Right side of initial velocity matrix
qdotR(9)=w;

B = inv(qdotc);

qdot = inv(qdotc)*qdotR % Initial Velocities for the given w

qddc           = zeros(9,9); % Coeffecient matrix of acceleration
qddc(1,1) =1;
qddc(1,3) =L1/2*sin(q(3));
qddc(2,2) =1;
qddc(2,3) = -L1/2*cos(q(3));
qddc(3,3) =L1*sin(q(3));
qddc(3,4)=1;
qddc(3,6) =L2/2*sin(q(6));
qddc(4,3) =-L1*cos(q(3));
qddc(4,5) =1;
qddc(4,6) =-L2/2*cos(q(6));
qddc(5,3) =L1*sin(q(3));
qddc(5,6) =L2*sin(q(6));
qddc(5,7) =1;
qddc(6,3) =-L1*cos(q(3));
qddc(6,6) =-L2*cos(q(6));
qddc(6,8) =1;
qddc(7,8) =1;
qddc(8,9)=1;
qddc(9,3)=1;

qddR =zeros(9,1); % Right side of acceleration matrix
qddR(1) =-qdot(3)^2*L1/2*cos(q(3));
qddR(2) =-qdot(3)^2*L1/2*sin(q(3));
qddR(3) = -qdot(3)^2*L1*cos(q(3))-qdot(6)^2*L2/2*cos(q(6));
qddR(4) =-qdot(3)^2*L1*sin(q(3))-qdot(6)^2*L2/2*sin(q(6));
qddR(5) =-qdot(3)^2*L1*cos(q(3))-qdot(6)^2*L2*cos(q(6));
qddR(6) =-qdot(3)^2*L1*sin(q(3))-qdot(6)^2*L2*sin(q(6));

E=inv(qddc);

qdd =E*qddR  % Intial accelerations
```
```

The result of this file is:

Now we need to simulate the position and velocities of the three bodies using constraint equation and differentiating them.

﻿

#### MATLAB Codes for Case 2:

``````
function zdot   = Project_2(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         = 0;                    % 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( 9,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;
Phiq(9,3)   =1;
%**********************************************************************
% Right Side of Constraint Equation
Phiqdot     =  zeros( 9,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( 9,9)];
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
qdot(1) = 0.3142;
qdot(2)=-0.5441;
qdot(3)=-2*pi;
qdot(4)=0.7394;
qdot(5)=-0.5441;
qdot(6)=2.2214;
qdot(7)=0.8505;

Z0        =  [q' qdot']';       % Initial Condition Vector

options   =  odeset('RelTol',1.0e-9,'AbsTol',1.0e-6);
[Tout,Z]  =  ode45(@Project_2,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');

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 2

﻿Your can find the solution for Case 1 by clicking here...