Cad Cam Development

Piston animation using Hamiltonian mechanics

The short video that presents overall behaviour of piston animation software that is the topic of this post:

The animation task consists of two smaller subtask. The first: to derive and solve the correct motion equation so that a real-life simulation could be provided and the second one: to provide UI and some visual objects that will move accordingly to that equation. This article will rather focus on motion derivation than how to change object’s position on canvas in WPF but note that full source code is provided so that you can get insight in what’s happening within C#.

The equation of motion can be derived using many ways. The fundamental course in classical mechanics would greatly facilitate the problem of equations of motion but you can gain knowledge also from wikipedia. There are many well-written articles, so read and spend time on learning! You can use Lagrange mechanics, Euler Hamilton etc. This problem could be easily solved using Lagrange mechanics, but I’m using Hamilton mechanics because this is one is the most versatile and it would be most educational to analyze the Hamilton mechanics in application to relatively easy problem of piston animation. There is quite good article about Hamiltonian mechanics on wikipedia so if you heaven’t read Goldstein yet and want to know what Hamiltonian mechanics is, read this article carefully. This post displays some steps in performing motion simulation of piston. Full source code and executable link is available at the very bottom of the article. Using WPF interface you can easily manipulate the following variables:

  1. Piston’s mass
  2. Wheel’s mass
  3. Rod’s mass
  4. Rod’s length
  5. Rod’s length error
  6. The length of rod is error bound. Error propagates accordingly to probabilistic normal distribution function
    f(t)=\frac{1}{\sigma\sqrt{2\pi}}\int^x_{-\infty} e^{-\frac{(t-\mu)^2}{2\sigma^2}} dt\ where:\\ \sigma\ -\ error\ standard\ deviation,\ user\ adjustable\\ \mu-\ mean\ value=0

  7. Initial angular wheel’s velocity

The program I’ve written provides assistance in piston motion simulation. Due to the accurate physical approach I presented here (Hamiltonian mechanics) this kind of software can be very educative. That’s great that using some numerical methods you can simulate almost any physical movement, with almost any complexity on your computer. This could be helpful either in understanding in what’s happening behind mechanical mechanism or just to simulate before production phase. In this example I was surprised to find out some movement mysteries like the change of wheel’s mass does not propagate either on the motion of piston or the motion of rod. There many real-time updatable charts so that you can get the full understanding on how the piston works so analyze them carefully.

To get more understanding of code and the derivation of motion equation see the photo below.

Piston schema

The derivation of motion equation:

(y_{s}-y_{p})^2=L^2\\\\x_{s}=x_{m}\\\\y_{s}=\sqrt{L^2-(x_{m}-x_{p})^2}+y_{m}\\\\\sin \psi = \frac {x_{m}-x_{p}}{L}= \frac{R\sin \phi}{L} \\\\\cos \phi = \frac{y_{s}-y_{p}}{L}\\\\    \psi=arcsin \frac{R\sin\phi}{L}\\\\    y_{s}(\phi)=\sqrt{L^2-R^2\sin^2\phi}+R\cos\phi+y_{m}\\\\    T=\frac{I_{1}\Omega^2_{\phi}}{2}+\frac{I_{2}\Omega^2_{\psi}}{2}+\frac{mv^2}{2}\\\\    U=m_{p}gh_{p}+m_{r}gh_{r}\\\\    Lagrange\ function: \\\\    L(q,\dot{q},t)=T(q,\dot{q},t)-U(q,t)\\\\    U(\phi,t)=mg(\sqrt{L^2-R^2\sin^2\phi}+R\cos\phi+y_{m})+m_{r}g(R\cos\phi+y_{m})\\\\    v_{s}(\phi)=-R\sin\phi-\frac{R^2\cos\phi\sin\phi}{\sqrt{L^2-r^2\sin^2\phi}}\\\\    T(\phi,\dot{\phi},t) =\frac{\dot{q}^2}{2}(I_{1}+I_{2})+\frac{m}{2}(R^2\sin^2\phi+\frac{R^4\cos^2\phi\sin^2\phi}{L^2-R^2\sin^2\phi}+\frac{2R^3\sin^2\phi\cos\phi}{\sqrt{L^2-R^2\sin^2\phi}})\\\\    Momentum\\\\    p= \frac{\partial H}{\partial \dot{\phi}}=\dot{\phi}(I_{1}+I_{2})\\\\    \dot{\phi}=\frac{p}{I_{1}+I_{2}}\\\\    Hamilton\ function:\\\\    H(p,\phi,t)=p\dot{\phi}-L

H(p,\phi,t)=\frac{p^2}{2(I_{1}+I_{2})}+\frac{m}{2}(R^2\sin^2\phi+\frac{R^4\cos^2\phi\sin^2\phi}{L^2-R^2\sin^2\phi}+\frac{2R^3\sin^2\phi\cos\phi}{\sqrt{L^2-R^2\sin^2\phi}})+\\    +mg(\sqrt{L^2-R^2\sin^2\phi}+r\cos\phi+y_{m})+m_{r}g(R\cos\phi+y_{m})\\\\    \dot{\phi}=\frac{\partial H}{\partial p}=\frac{p}{I_{1}+I_{2}}\\    p=\dot{\phi}(I_{1}+I_{2})\\\\    \dot{p}=-\frac{\partial H}{\partial \phi} = -gmR\sin\phi-R\sin\phi-\frac{R^2\cos\phi\sin\phi}{\sqrt{L^2-R^2\sin^2\phi}}+\frac{m}{2}(2R^2\cos\phi\sin\phi+\\    +\frac{2R^4\cos\phi\sin\phi}{L^2-r^2\sin^2\phi}+\frac{2R^5\cos^2\phi\sin^3\phi}{(L^2-R^2\sin^2\phi)^\frac{3}{2}} +\frac{4R^3\cos^2\phi\sin\phi}{\sqrt{L^2-R^2\sin^2\phi}}-\frac{2R^3\sin^3\phi}{\sqrt{L^2-R^2\sin^2\phi}} )\\\\    \dot{p}=\ddot{q}(I_{1}+I_{2})\\\\    Implicit\ move\ equation:\\\\    \ddot{q}(I_{1}+I_{2})= -gmR\sin\phi-R\sin\phi-\frac{R^2\cos\phi\sin\phi}{\sqrt{L^2-R^2\sin^2\phi}}+\frac{m}{2}(2R^2\cos\phi\sin\phi+\\    +\frac{2R^4\cos\phi\sin\phi}{L^2-r^2\sin^2\phi}+\frac{2R^5\cos^2\phi\sin^3\phi}{(L^2-R^2\sin^2\phi)^\frac{3}{2}} +\frac{4R^3\cos^2\phi\sin\phi}{\sqrt{L^2-R^2\sin^2\phi}}-\frac{2R^3\sin^3\phi}{\sqrt{L^2-R^2\sin^2\phi}} )\\\\

Quick verification of obtained motion equation in Mathematica:

To verify the equation you can check it using Matlab or Mathematica. Personally I’m devoted Mathematica user. I entered the equation from above and solved it numerically. The results presents photo below. There is time on x axis and angle on y axis. You can clearly see that the angle is growing linearly according to time. That’s what we needed. Implementing this equation in your code should provide correct results. You could even try animating it in Mathematica environment. Click the photo below to enlarge and see the details: bigger chart and command I used for evaluation.

Numerical solving of piston's Hamilton move equation

Comment on implicit differential equations

As you have already noticed the move equation is not linear. The angle which you solve for is time dependent and implictly called within trigonometric functions. This makes solving this equation difficult. There is widely-known and solved problem of analytical solution of pendulum movement. That’s astonishing how complicated analytical solution it has compared to its obviousness and similarly it’s with the pendulum problem. The substantial help comes from numerical methods which make it a bit easier.

How I solved this equation in my program?

Firstly, I substituted another variable for the first angle derivative, got first order nonlinear equation and solved for mass angular velocity. Secondly, according to the angular move analysis I solved another first order equation. all in one: I dispatched the second ODE into the system of two first order ODE equations and made it work. I used Runge Kutta 4th order to solve both of the equations of the system but the results could be more accurate. Analyze this page to get more understanding. The system I implemented is stable and does not take on energy, what is extremely important in solving this kind of work. Personally, after I completed this work, I would like to try to solve the equations with 6th order method like Runge-Kutta-Pohlberg in order to get smaller phase difference between two following steps. For this reason Runge-Kutty-Pohlberg would be the best method in solving the piston’s move equation if for instance Ducatti had asked me to do some pre-production simulation.

Ducatti 848 engine

Source code considerations

As I’ve written previously the implicit move equation is solved in two steps. Firstly I solved the angular velocity and secondly the angle using the ODE:

\omega=\frac{d\theta}{dt}\\\\v=\omega r\\\\    \int \mathrm{d}\theta=\int \frac{v}{r} \mathrm{d}t

Source code that solves the equation below.

 public void solve(ref double time)
            //solve angular velocity
            double k1 = getf( this.Wheel.Angle, time );
            double k2 = getf( this.Wheel.Angle + this.Step * k1 / 2.0, time + this.Step / 2.0 );
            double k3 = getf( this.Wheel.Angle + this.Step * k2 / 2.0, time + this.Step / 2.0 );
            double k4 = getf( this.Wheel.Angle + this.Step * k3, time + this.Step );
            double k = this.Wheel.AngularVelocity + this.Step * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;
            this.Wheel.AngularVelocity = k;

            //solve angle
            double a1 = getfAngle( this.Wheel.AngularVelocity, time );
            double a2 = getfAngle( this.Wheel.AngularVelocity + this.Step * a1 / 2.0, time + this.Step / 2.0 );
            double a3 = getfAngle( this.Wheel.AngularVelocity + this.Step * a2 / 2.0, time + this.Step / 2.0 );
            double a4 = getfAngle( this.Wheel.AngularVelocity + this.Step * k3, time + this.Step );
            double a = this.Wheel.Angle + (a1 + 2.0 * a2 + 2.0 * a3 + a4) / 6.0;
            this.Wheel.Angle = a;

            time += Step;

UML Diagram for piston animation

I’m using Charles Petzold library to provide visual objects I’m simulating. It’s very easy to write this kind of simple geometrical animation. I used Cylinder base class for all my simulating objects and it really comes in handy.. The UML diagram of my classes is above.

3D adjustable camera

If you click minimize/maximize you can see that the size of simulating objects is automatically adjusting. This means, that it will comply with the size of the screen that the software runs on. That’s very nice and helpful feature. Every professional CAD/CAM or computer graphics software has it. I implemented it by changing the camera position, but you play with FOV also. If you see the code below you’ll see that I use power function so that smaller screens can get smaller increase in camera position than bigger screens. This regulates the size of simulating objects on the screen, because it changes the projection matrix.

private void Viewbox_SizeChanged(object sender, SizeChangedEventArgs e)
 double aspect = e.NewSize.Width / e.NewSize.Height;
 camera.Position = new System.Windows.Media.Media3D.Point3D( camera.Position.X, camera.Position.Y, 7.0 + Math.Pow( aspect * 2.2, 2 ) );
 double width = this.renderCol.ActualWidth;
 double height = this.row1.ActualHeight;
 this.viewport.SetValue( WidthProperty, width );
 this.viewport.SetValue( HeightProperty, height );
 this.renderCanvas.SetValue( WidthProperty, width );
 this.renderCanvas.SetValue( HeightProperty, height );

Source code and executable download:
Source code is downlodable from here and executable from here


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s