Piston animation using Hamiltonian mechanics

pytelg

Rate me:

0.00/5 (No votes)

14 Feb 2011CPOL5 min read

13.6K

Piston animation using Hamiltonian mechanics

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

The animation task consists of two smaller subtasks. 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 the 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 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 a good article about Hamiltonian mechanics on Wikipedia so if you haven’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. The 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:

Piston’s mass
Wheel’s mass
Rod’s mass
Rod’s length
Rod’s length error
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$
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 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 are 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 a devoted Mathematica user. I entered the equation from above and solved it numerically. The results presents the 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 implicitly called within trigonometric functions. This makes solving this equation difficult. There is a 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 4^th 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 6^th 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.viewport.UpdateLayout();
 this.renderCanvas.SetValue( WidthProperty, width );
 this.renderCanvas.SetValue( HeightProperty, height );
 this.renderCanvas.UpdateLayout();
 this.viewbox.UpdateLayout();
 }

The source code is downlodable from here and executable from here.

Filed under: C#, CodeProject, WPF Tagged: Hamilton, Mechanics, WPF

This article was originally posted at http://cadcammodelling.wordpress.com/2011/02/14/piston-animation-using-hamiltonian-mechanics

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Written By

pytelg

Student

Poland

I'm student of CAD/CAM design at Warsaw University of Technology, currently at the last 5th course year and successfully graduated from BSc studies in computer science in February 2010 from the same university. I'm interested in windows programming technolgies, geometric modelling, programming of numerically controlled devices, virtual reality and computer graphics. You can follow me on facebook : http://www.facebook.com/profile.php?id=1249817870 or linkedIn: http://pl.linkedin.com/in/pytelg.

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

Piston animation using Hamiltonian mechanics

The Derivation of Motion Equation

Quick Verification of Obtained Motion Equation in Mathematica

Comment on Implicit Differential Equations

How I Solved this Equation in my Program?

Source Code Considerations

Source Code that Solves the Equation Below

3D Adjustable Camera

License

Comments and Discussions