Solution procedure of piston problem using Matlab

Problem Statement

We are going to use the mechanism (figure below) of a 2-bar link and a piston. The problem we will look at is that, constant moment (torque), M, is applied to the first link (left bar, bar 1) and we want to know the displacement of piston, x, as a function of time. No damping or friction exists between the piston and the walls in this problem.

Free Body Diagram and Equation of Motion

We start with examining the kinematics of the problem first. Current problem has certain kinematic restraints and we need those to setup the equations of motion of the system. When we setup the equations of motion, we need FreeBody Diagram (FBD). Generally, when solving dynamic problems (or static problems), you have to draw the FBD first since it will give you a clear view of the problem. With correct FBD's, you can write down the equations for both force balance and moment balance more easily and you will end up with differential equations, which are equations of motion of the system. So, don't overlook the importance of FBD.

The procedure for developing the equations of motion including the kinematics of the system is shown

First we will examine the kinematics of the problem. Since the mass is restricted by the linkages, the following equations are obtained.
\begin{displaymath}l_1 \mbox{sin}\theta = l_2 \mbox{sin}\phi\end{displaymath} (1)
 
\begin{displaymath}x= l_1 \mbox{cos}\theta + l_2 \mbox{cos}\phi\end{displaymath} (2)
From these relations we can write some more relations for velocities in terms of the three degrees of freedom.
\begin{displaymath}l_1 \dot{\theta}\mbox{cos}\theta = l_2 \dot{\phi}\mbox{cos}\phi\end{displaymath} (3)
 
\begin{displaymath}\dot{x}= -\left(l_1 \dot{\theta}\mbox{cos}\theta + l_2 \dot{\phi}\mbox{cos}\phi\right)\end{displaymath} (4)
However, we do not need the expressions about $\dot{\theta},\; \dot{\phi},\;\mbox{and} \; \dot{x}$ in the current problem since the external force is a constant torque acting on the first linkage. We need only one equation of motion for $\theta$, using the kinematic relationships to eliminate the other variables.

The equation of motion for $\theta$ can be written as follows
\begin{displaymath}I_1\ddot{\theta}=M\end{displaymath} (5)
where I1 is the moment of inertia for linkage 1 (a given quantity). With a very simple equation of motion and using the ODE45 solver in MATLAB we can determine all the variables in the problem. Note that from the equation of motion we can guess the response of the system. Rewriting the EOM as
\begin{displaymath}\ddot{\theta}=\frac{M}{I_1}\end{displaymath} (6)
hilights the fact that the angular acceleration is a constant. Therefore, the angular velocity will increase linearly and the position will change as t2.
 
 

After we get the equations of motion, we need to solve them with given initial conditions. There are several ways to solve the ordinary differential equations, however, only the integration method (using ODE45 solver in Matlab) is explained below.

Closed Form Solution

Due to the complexity of this problem, we will not find a closed form solution. This is a case where using a numerical method is best. 

Integrating EOM with Matlab

How do we solve Ordinary Differential Equations (ODE) using Matlab? Do we have to write numerous codes to do that? The answer is NO. Life is not that miserable. There are several ways in Matlab to solve ODE's, but we are going to discuss only about typical ODE solver in Matlab, ODE45, which is rigorous enough for our purpose.

So, How Do We Use ODE45 Solver?

 ODE45 uses a 4th and 5th order Runge-Kutta method to integrate the equations of motion that we already got from the previous step. Don't worry, you don't have to fully understand 4th and 5th order Runge- Kutta method to use ODE45 solver, but it is recommended you to read about the method to get the full benefit out of using ODE45.

ODE45 has a typical form to use (Matlab version 4.2c1 only since version 5.0 has a little bit different format). 

[T, Y] = ODE45(F, T0, Tfinal, Y0, TOL)
T in the left bracket is for time, Y is our main variables in vector form. For example, if we want to know the displacement of a piston and use "x" as a displacement variable, "x" will be the main variable in the equations of motion and that "x" correspond to Y in the left bracket. Now, F in the right side is a name of the subprogram for ODE45 solver (will be explained in detail a little bit later). T0 is the starting point of calculation time, usually 0. Tfinal is final calculation time. Y0 is the initial conditions for ODE variable, Y. TOL is the tolerance for the error calculation in ODE45 solver (default value is 1e-6), lower it is, more accurate the answer will be, but slower the calculation is.
Now, the tricky part is how to use a subprogram right. Example will be the best way to learn how to do it right.

First, we are going to look at the main program and then the subprogram. The following is a Matlab source code for main program. The first few lines are for constant input values like mass of the piston, etc. The numbers given are all arbitrary. 

m = 1;                      %  mass of the piston
l1 = 1;                     %  length of the bar1 (left bar)
l2 = 2;                     %  length of the bar2 (right bar)
M = 5;                      %  constant moment at bar1
I1 = 3;                     %  inertia of bar1

t_final = 10;               %  time duration for calculation
                            %  initial time is always assumed as 0 
                            %  in this problem
As you guess, "%" is the comment symbol, so anything after "%" is negelcted for program execution.
Now, for initial conditions for main variable, theta in this case, the rotation angle of bar 1. 
theta = 0;                  %  initial displacement of the bar
theta_dot = 0;              %  initial angular velocity of the bar1
theta_ddot = M/I1;          %  initial angular acceleration of the bar1
In this problem, we apply the constant torque on the left link (bar 1), so, the angular acceleration of the bar 1 will be constant and the magnitude is M/I1. Furthermore, we need to specify tolerance and make initial conditions in a vector (array) form, which is just a syntax, meaning you don't have any other choices. 
tol = 1e-20;                %  tolerance for ODE45 solver
                            %  smaller it is, more accurate the solution
                            %  will be
x0 = [theta,theta_dot,theta_ddot];
Going to very high accuracy (very small tolerance, as you see above) for this problem gives nice looking plots. With tolerance of 1e-20 for this problem, it'll take about 2-4 minutes depending on what machine you're using while with tolerance of 1e-10 will give you the results right away no matter what machine you are using, but poor-looking plots. So, if you want to just check out whether the program runs and don't want to wait, reduce calculation time and take high tolerance (like 1e-10). BTW, with other problems, you usually will only need the default tolerance of 1e-6 to get some good-looking plots. After we specify all the inputs for ODE45 solver, we write 
[t,x] = ode45('pist_sol',0,t_final,x0,tol);
Here, "pist_sol" is the name of the subprogram and it is shown below. 
function xdot = pist_sol(t,x)

%
%   State variables
%

x_dot1 = x(2);
x_dot2 = x(3);
x_dot3 = 0;

xdot = [x_dot1; x_dot2; x_dot3];
The first line is for declaring subprogram. Array x is [x(1), x(2), x(3)], which is [displacement, velocity, acceleration]. And array xdot is for a derivative of array x, so xdot(1) is a velocity, xdot(2) is an angular acceleration, xdot(3) is a impulse, or point peak, the change in acceleration in time. Since in this problem, the angular acceleration is constant, the change in the angular acceleration is zero, i.e. xdot(3) = 0 all the time.

Now you got it all. You just write them down in the subprogram as above. It is that simple. Even for somewhat complicated problem, it is still not that difficult to write this relations (and they are called state variables) of the variables and their derivatives.

So, let's summarize things. First, you write some constant input values in your main program, like the mass of a piston, etc. Then, you need to write the initial conditions for your variables ( anything can be your variables, like displacement, rotation angle, velocity, etc. ). You write them in array format and you specify the tolerance as small or large as you want (remember the default is 1e-6 and it is good enough for most cases). After all that, you write "[t,x]..." in the above. That's all for the main program. For the subprogram, you declare that this is a subprogram. And then, you write all the derivatives of your variables in state variable format as above. Finally, you write those state variables in array format as the above, "xdot = [x_dot1 ...];". Give the name for the subprogram, the same name you used for the subprogram in line "[t,x] = (...)" in the main program and save it.
Plotting Results

After finishing the calculation, you plot them. Here are some of the plotting commands. They should be in the last part of the main program. 

figure(1); clf; orient tall;
subplot(2,1,1),plot(t,Theta);
title('Problem 3.175'); 
xlabel('Time (sec.)');  
ylabel('Theta');

subplot(2,1,2),plot(t,Phi);
xlabel('Time (sec.)');    
ylabel('Phi');
  
figure(2); clf; orient tall;
subplot(2,1,1),plot(t,X);
title('Problem 3.175');
xlabel('Time (sec.)');
ylabel('X (m)');

subplot(2,1,2),plot(Theta,X);  
xlabel('Theta (rad.)');
ylabel('X (m)');

mxx = max(X);       mnx = min(X);
mxt = max(Theta);   mnt = min(Theta);
axis([mnt,mxt,mnx,mxx])
In the first line, command "clf" clears the figure, "orient tall" gives the vertically tall plot, rather than the default landscape type plot. But, you will not see the difference on the screen between landscape one and vertical one. The difference happens when you print them with "print" command. For Unix users, you can make the plot as .ps file using print command: print -deps (your directory here, for example, /afs/engin.umich.edu/u/m/ a/marimba/Matlab/) filename. So, it will be like "print -deps /afs/engin.umich.edu/u/m/a/marimba /Matlab/example.ps". For PC and Macintosh version of Matlab, you can choose from the menu to print out plots. Now, I believe you get the ideas for next few lines, xlabel, ylabel and title commands. The subplot command divides a figure with one plot into the figure with several plots. Subplot(number of rows of plots, number of columns of plots, actual location of the current plot).

The command "max(X)" gives the maximum value of the array X. "min(X)" gives minimum value of the array X. The command "axis([x_lowest limit, x_highest limit, y_lowest limit, y_highest limit])" gives the limits for axis on the plot. Note that I used "Theta", "Phi" and "X" here. I just changed the names of the columns of array x, the original array we used in main program. You can see how or what I did in main program source code below.

NOTE that Matlab on Unix workstations is version of 5.0 and MatlabS on PC's and Macintosh's are 4.2 c1. If you run the program for version 4.2 c1 on Matlab 5.0, it'll still run and give you plots, but it'll also give you some error messages. Although, main programs are a little bit different from ver 4.2 c1 to ver 5.0, both can share the subprogram.

Here, the complete source programs of main program (PC,Mac), main program (Unix) and subprogram (PC,Mac,Unix).

Some results are shown in the following. Here, "Theta" is the rotation angle of first link, bar 1, "Phi" is the rotation angle of bar 2 and "X" is the displacement of piston and they are all in radians. As we saw in the constant input part (the very first part of the main program), the length of the bar 1 was 1 (m), the length of the bar 2 was 2 (m).

 

As you see in the first plot, "Theta" increases with time, but "Phi" fluctuates since the mechanism does not allow the second link to rotate, so it goes up and down. Also, you can notice that displacement "X" is getting faster as time goes, (a period of the sinusoidal movement of "X" decreases) since the moment at the bar 1 accelerates the piston.

NOTE that If you're trying to reproduce these results above, remember that the calculation will take some time (2-4 minutes) depending on what machine you're using. If the calculation time looks to long to you, then reduce the t_final to small value (originally it's 10 seconds), for example, 5 seconds.

If you need some more information about commands, try "help". For example, 

help plot
or the help command will give some information regarding any command. Also see the other Matlab help for this class here .

Return to Example Table