The GUI based simulation was written in Matlab 2011a. A number of animated GIF files were made. In the table below, clicking on any image will start a gif animation in your browser.
These are relatively large animation GIF files and might take 1-2 seconds to load depending on network bandwidth. To run them locally, they can be downloaded as well.
|
| ||
|
Using low spring stiffness to reach 5000 RPM in one second. (5 MB) |
Using low spring stiffness, disk is heavier than above, was not able to reach 5000 RPM in 1 second. (600 KB) |
Same as above, but the torque was removed at 500 RPM. |
The Lagrangian energy method was used to derive the equations of motion. The physical model used to represent the connection between the piston and the crank arm uses the method of Knarnopp-Margolis. This method uses a virtual spring and virtual damper for the connection between the piston and the crank arm instead of a rigid connection. This leads to a simpler mathematical model.
The equation of motion will now be derived based on this method.
Let \(f(\theta )\) be the distance which is shown in the diagram above as \(x^{\prime }\), The kinematic constraint is the following\[ f\left ( \theta \right ) =R\sin \theta +\sqrt{L^{2}-R^{2}\cos ^{2}\theta }\] The extension in the spring is \[ \Delta =x-f\left ( \theta \right ) \] Hence the spring force is\begin{equation} F_{s}=k\left [ x-f\theta )\right ] \tag{1} \end{equation} The damping force is\begin{align} F_{d} & =b\dot{x}\nonumber \\ & =b\frac{d}{dt}\left ( x-f\left ( \theta \right ) \right ) \nonumber \\ & =b\left ( \dot{x}-\frac{df\left ( \theta \right ) }{d\theta }\dot{\theta }\right ) \tag{2} \end{align}
The system is a 2 DOF, with generalized coordinates \(\theta ,x\).
Let \(I=\frac{MR^{2}}{2}\) be the moment of inertia of the disk (called \(J_{fw}\) in the above diagram) around an axis through the disk center, where \(M\) is the mass of the disk (called \(m_{fw}\) in the diagram). Let \(m\) be the mass of the cylinder (called \(m_{p}\) in the problem diagram).
The kinetic energy of the system is\[ T=\frac{1}{2}I\dot{\theta }^{2}+\frac{1}{2}m\dot{x}^{2}\] The potential energy is\begin{align*} V & =\frac{1}{2}k\Delta ^{2}\\ & =\frac{1}{2}k\left ( x-f\left ( \theta \right ) \right ) ^{2} \end{align*}
Therefore the Lagrangian is\begin{align*} L & =T-V\\ & =\frac{1}{2}I\dot{\theta }^{2}+\frac{1}{2}m\dot{x}^{2}-\frac{1}{2}k\left ( x-f\left ( \theta \right ) \right ) ^{2} \end{align*}
For \(\theta \) we find\begin{align*} \frac{\partial L}{\partial \dot{\theta }} & =I\dot{\theta }\\ \frac{d}{dt}\frac{\partial L}{\partial \dot{\theta }} & =I\ddot{\theta } \end{align*}
and\begin{align*} \frac{dL}{d\theta } & =-k\left ( x-f\left ( \theta \right ) \right ) \left ( -\frac{df\left ( \theta \right ) }{d\theta }\right ) \\ & =k\left ( x-f\left ( \theta \right ) \right ) \frac{df\left ( \theta \right ) }{d\theta } \end{align*}
Therefore the EQM for \(\theta \) is\begin{equation} \frac{d}{dt}\frac{\partial L}{\partial \dot{\theta }}-\frac{dL}{d\theta }=Q_{\theta }\tag{3} \end{equation} To find the generalized force \(Q_{\theta }\), we make a virtual \(\delta \theta \) displacement while keeping \(x\) fixed and then find the virtual work done. This results in\begin{align*} \delta W & =\tau \delta \theta +F_{d}\left ( \delta f\left ( \theta \right ) \right ) \\ & =\tau \delta \theta +F_{d}\frac{df\left ( \theta \right ) }{d\theta }\delta \theta \\ & =\left ( \tau +F_{d}\frac{df\left ( \theta \right ) }{d\theta }\right ) \delta \theta \end{align*}
Hence from the above we see that \[ Q_{\theta }=\tau +F_{d}\frac{df\left ( \theta \right ) }{d\theta }\] Eq. (3) becomes\begin{equation} I\ddot{\theta }-k\left ( x-f\left ( \theta \right ) \right ) \frac{df\left ( \theta \right ) }{d\theta }=\tau +F_{d}\frac{df\left ( \theta \right ) }{d\theta }\tag{4} \end{equation} Now we evaluate \(\frac{df\left ( \theta \right ) }{d\theta }\)\begin{align} f\left ( \theta \right ) & =R\sin \theta +\sqrt{L^{2}-R^{2}\cos ^{2}\theta }\nonumber \\ \frac{df\left ( \theta \right ) }{d\theta } & =R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\tag{5} \end{align}
Substituting the above back into Eq. (2) gives\begin{equation} F_{d}=b\left ( \dot{x}-\left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \tag{6} \end{equation} Substituting Eqs. (6,5) into Eq. (4) gives \begin{multline*} I\ddot{\theta }-k\left ( x-f\left ( \theta \right ) \right ) \left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) =\\ \tau +\left ( b\left ( \dot{x}-\left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \right ) \left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \end{multline*} Simplifying the above gives the final EQM for \(\theta \)\begin{align} \ddot{\theta } & =\frac{\tau }{I}+\frac{kR}{I}\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \left ( \sqrt{L^{2}-R^{2}\cos ^{2}\theta }-R\sin \theta +x\right ) \nonumber \\ & +\frac{bR}{I}\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \left ( \dot{x}-R\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \tag{7} \end{align}
Now to find the EQM for \(x\)\begin{align*} \frac{\partial L}{\partial \dot{x}} & =m\dot{x}\\ \frac{d}{dt}\frac{\partial L}{\partial \dot{x}} & =m\ddot{x} \end{align*}
and\[ \frac{dL}{dx}=-k\left ( x-f\left ( \theta \right ) \right ) \] Therefore the EQM for \(x\) is\begin{equation} \frac{d}{dt}\frac{\partial L}{\partial \dot{x}}-\frac{dL}{dx}=Q_{x}\tag{8} \end{equation} To find \(Q_{x}\) we make a virtual \(\delta x\) displacement while keeping \(\theta \) fixed and find find the virtual work. Hence\[ \delta W=-F_{d}\delta x \] Therefore \(Q_{x}=F_{d}\) which was found in Eq. (2) above. Hence \[ Q_{x}=-b\left ( \dot{x}-\frac{df\left ( \theta \right ) }{d\theta }\dot{\theta }\right ) \] And Eq. (8) becomes\[ m\ddot{x}+k\left ( x-f\left ( \theta \right ) \right ) =-b\left ( \dot{x}-\frac{df\left ( \theta \right ) }{d\theta }\dot{\theta }\right ) \] Substituting the expressions for \(f\left ( \theta \right ) \) and \(\frac{df\left ( \theta \right ) }{d\theta }\) found earlier into the above results in
\[ m\ddot{x}+k\left ( x-\left ( R\sin \theta +\sqrt{L^{2}-R^{2}\cos ^{2}\theta }\right ) \right ) =-b\left ( \dot{x}-\left ( R\cos \theta +\frac{R^{2}\cos \theta \sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \] Therefore\begin{equation} \ddot{x}=\frac{k\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}{m}+\frac{kR\sin \theta }{m}-\frac{kx}{m}-\frac{b\dot{x}}{m}+\frac{bR\dot{\theta }\cos \theta }{m}+\frac{b\dot{\theta }R^{2}\cos \theta \sin \theta }{m\sqrt{L^{2}-R^{2}\cos ^{2}\theta }} \tag{9} \end{equation}
Now Eqs. (7) and (9) above will be cast into first order form.\[\begin{pmatrix} x_{1}\\ x_{2}\\ x_{3}\\ x_{4}\end{pmatrix} =\begin{pmatrix} \theta \\ x\\ \theta ^{\prime }\\ x^{\prime }\end{pmatrix} \overset{\frac{d}{dt}}{\rightarrow }\begin{pmatrix} \dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\\ \dot{x}_{4}\end{pmatrix} =\begin{pmatrix} x_{3}\\ x_{4}\\ \theta ^{\prime \prime }\\ x^{\prime \prime }\end{pmatrix} \] Hence\begin{align*} \begin{pmatrix} \dot{x}_{1}\\ \dot{x}_{2}\\ \dot{x}_{3}\\ \dot{x}_{4}\end{pmatrix} & =\begin{pmatrix} x_{3}\\ x_{4}\\ \frac{\tau }{I}+\frac{kR}{I}\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \left ( \sqrt{L^{2}-R^{2}\cos ^{2}\theta }-R\sin \theta +x\right ) \\ +\frac{bR}{I}\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \left ( \dot{x}-R\cos \theta \left ( 1+\frac{R\sin \theta }{\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\right ) \dot{\theta }\right ) \\ \frac{k\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}{m}+\frac{kR\sin \theta }{m}-\frac{kx}{m}-\frac{b\dot{x}}{m}+\frac{bR\dot{\theta }\cos \theta }{m}+\frac{b\dot{\theta }R^{2}\cos \theta \sin \theta }{m\sqrt{L^{2}-R^{2}\cos ^{2}\theta }}\end{pmatrix} \\ & =\begin{pmatrix} x_{3}\\ x_{4}\\ \frac{\tau }{I}+\frac{kR}{I}\cos x_{1}\left ( 1+\frac{R\sin x_{1}}{\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}\right ) \left ( \sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}-R\sin x_{1}+x_{2}\right ) \\ +\frac{bR}{I}\cos x_{1}\left ( 1+\frac{R\sin x_{1}}{\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}\right ) \left ( x_{4}-R\cos x_{1}\left ( 1+\frac{R\sin x_{1}}{\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}\right ) x_{3}\right ) \\ \frac{k\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}{m}+\frac{kR\sin x_{1}}{m}-\frac{kx_{2}}{m}-\frac{bx_{4}}{m}+\frac{bRx_{3}\cos x_{1}}{m}+\frac{bx_{3}R^{2}\cos x_{1}\sin x_{1}}{m\sqrt{L^{2}-R^{2}\cos ^{2}x_{1}}}\end{pmatrix} \end{align*}
The system was simulated on the computer. The first order equations above were solved using ode numerical solver. A GUI program was written to make it easier to do the simulation and observe the results. The following describes the simulation program user interface
We are asked to find the torque needed to reach \(5000\) RPM in one second. By trial and error, a torque of \(\tau =7\) \(NM\) was found to meet this condition, with \(f=10Hz\), and damping ratio \(\xi =0.2\), and using the same other parameters shown in the problem table. Using this torque, the \(5000\) RPM was reached at about \(0.5\) seconds. The force on the piston started to become lower after the torque was removed, as well as the piston acceleration.
The following diagram shows the final plots at the end of the one second run.
We are also asked to find the parameters which will cause the spring deflection to remain below fraction of a millimeter. It was found that the spring \(k\) must be made very large to simulate the effect of a rigid connection between the piston and the crank arm in order to keep the deflection very small. But when this was done, the simulation started to take very long time. Simulating \(1\) second on my computer would more than one hour.
To reduce vibration, the damping coefficient \(\xi \) was made larger than the \(0.2\) value in the table. But even when making \(\xi \) close to one, the simulation still took very long time to reach one second due to the large number of steps taken by the numerical solver. The spring stiffness was increased all the way to \(10^{7}\) N/m, and the force on the piston reached \(1.5MPa\) at one point. The value of \(f\) (spring frequency) to achieve this high spring stiffness was \(1000\) \(Hz.\) The damping ratio \(\xi \) was set to \(0.9\). The following is a plot of the simulation using the above parameters, showing that the deflection remained below half millimeter.
The acceleration on the piston was found numerically at each step by dividing the difference of the current speed of piston from the last time step speed, and dividing that over the difference in time as follows\[ a_{x}=\frac{\dot{x}\left ( t_{n}\right ) -\dot{x}\left ( t_{n-1}\right ) }{t_{n}-t_{n-1}}\] The force on the piston was found at each step as follows\[ F=-(F_{s}+F_{d}) \] Where \(F_{s}\) is the force in spring and \(F_{d}\) is the damping force. Expression for these forces were derived above.
The force on the piston depended on the torque used and the amount of spring stiffness and damping used. For example, in the first test case above when the spring was made very stiff, the largest force on the piston reached was \(1.5\times 10^{6}N\) and an acceleration of almost \(1000g\) was seen. After the torque was removed, the acceleration gradually dropped to \(500g\). The longer the simulation was run, these values became smaller. The simulation program updates both the acceleration and the force on the piston as the simulation is running.
The Knarnopp-Margolis model of the crank/piston connection was used to simulate the motion. This model is simpler mathematically than the rigid connection model. But it was hard to find the right parameters to approximate a rigid connection since making the spring stiffness large, causes the simulation to become very slow. When the spring stiffness was made very large to better approximate a rigid connection, and when the damping was made larger, the simulation would take about \(1\) hour for \(1\) second simulation time. Even when using a stiff ode solver such as ode15s the simulation was taking too long. Using a compiled language would eliminate this problem as it will be much faster, but this simulation was done using Matlab
The main advantage of the Knarnopp-Margolis model is that the implementation was simpler since the mathematical model was simpler. But more trial and error are needed to find optimal values for all the parameters of the problem which will cause the simulation time to become more reasonable, while approximating the rigid connection.
