We are going to use the mechanism (figure below) of a disk with radius L rotating on a fixed axis through its center and a massless link of length L connected to a massless piston. The problem we will look at is that for a constant moment (torque), M₀, applied to the disk, find the differential equation of motion for theta and the expression for x as a function of theta. No damping or friction exists between the piston and the walls in this problem.   However, there is gas sealed in the cylinder which compresses from an original volume of 4LA where A is the cross-sectional area of the cylinder.  To find the pressure in the cylinder and therefore the force on the piston, we will use the ideal gas law, PV=nRT, where P is the pressure in the cylinder, V is the volume of the cylinder, n is the number of moles of gas, R is the ideal gas constant, and T is the temperature in Kelvin which we will assume remains constant.  Recall that P=F/A when setting up the equations.

Extra Problem: derive the equations of motion for this system by drawing free body diagrams for each of the three bodies in this problem.  Assume the disk has mass, m, and the cross-sectional area of the cylinder is A (this will not appear in your final equation of motion).  Also, so there is no confusion, assume that the cylinder is not in the same plane as the disk, i.e. the piston is free to travel the entire length of the cylinder.   The results here will be used later in a computer simulation.

MATLAB Program

Following the example on the Web, use ODE45 in MATLAB to simulate the piston problem.  The correct equation of motion was handed out in lecture. Use the following parameters in your program:

m_disk = 1 kg
L = 0.1 m
n = 0.05 moles
R = 8.314 J/(mol*K)
T = 300 K
tspan = 0:0.0005:.2;

Assume the disk starts at rest and at an initial angle of 0 radians.

Vary the moment M from a value of 5 Nm to 25 Nm to note the change in behavior of the system for different moments. Does the system behave as you would expect? Use MATLAB to make plots of x (the displacement of the piston) versus time and theta versus time for two different moments showing two different types of behavior. Explain the differences between the plots.

Use the following equation of state within your MATLAB subroutine:

state = [((1/abs(cos(theta(2))))*((nRT/(4*l-l*cos(theta(2))-abs(l*cos(theta(2)))) ...
  - nRT/(2*l))*(l*abs(cos(theta(2)))*sin(theta(2))+l*sin(theta(2)) ...
  *cos(theta(2)))) - M)/I; theta(1)];

Note: The elipses (...) act to allow a line of code to continue on the next  typed line.  Write this equation out by hand to see how the equation differs from (or is the same as) what you solved for in part a.

Also, x=l*(cos(theta)+abs(cos(theta)))