Ideal Gas Simulator: Help Solving Collision Detection Issue

[JizzaDaMan]

This one's aimed primarily at the programmers in this forum.

I've written a programme in python (Set as a long term project from my A-level physics teacher) which will hopefully end up being a scientifically accurate ideal gas simulator in 2D - I'll set up a bunch of particles with random velocities, then let it run indefinitely. I'll write the kinetic energy of each particle to a file every 100 'time steps' (or data to that effect) that I can then plot a histogram from in excel to analyse the distributions of the velocities.

At the moment I'm at a stage where I have multiple particles bouncing around a box, and I'm just showing it graphically so I can see what's happening.

Unfortunately at the moment I have a pretty serious problem with my code; at seemingly random collisions between particles, total kinetic energy is not conserved. When all the particles are the same mass and without gravity (which I've tried implementing in the past to see if it works) the increase in energy that happens isn't actually noticeable to the eye (but still not scientifically accurate!), but when I change some of the particles' masses or add in gravity then it all just goes totally crazy.

I've spoken to my physics teacher about it and he seems to think that when I detect a collision, it's never going be the exact moment of impact, so I need to extrapolate backwards to determine the exact time of the collision and work things out from there. I've looked at source code of other programmes achieving the same thing however, and I haven't seen anything so complicated.

So yeah this problem is totally confusing me, so some help would be great. This is the collision detection routine which then calculates the velocities after the collision:

	def CollideBall(self, ball, timeframe = 1):				
		Impulse = self.VectorToBall(Vector(ball._getX() + ball.Velocity._getX() * timeframe, 
			ball._getY() + ball.Velocity._getY() * timeframe))
		
		if Impulse.Length() <= self.Size + ball.Size:
			Impulse.Normalize()
			Impact = Vector(self.Velocity._getX() - ball.Velocity._getX(), self.Velocity._getY() - ball.Velocity._getY())
			ImpactVelocity = Impact.DotProduct(Impulse)
			
			Impulse.MultiplyScalar(math.sqrt(abs(ImpactVelocity**2 * self.Mass * ball.Mass)))
			
			self.Velocity.X += Impulse.X / self.Mass
			self.Velocity.Y += Impulse.Y / self.Mass

			ball.Velocity.X -= Impulse.X / ball.Mass
			ball.Velocity.Y -= Impulse.Y / ball.Mass

I'm afraid I'm no professional programmer so my code's a bit messy, I hope you can follow it. All the variable names should be fairly self explanatory.

I implemented that by following this semi-tutorial thingy:
http://freespace.virgin.net/hugo.elias/models/m_snokr.htm

I had to change one little bit; notice that I squared the impact velocity on the 'Impulse.MultiplyScalar()' line, whereas the tutorial doesn't mention anything like that. I did that because before kinetic energy wasn't being conserved in any collisions at all.

The rest of my code is attached if you want to read it. In 'Physics.py' are the methods and classes that deal with the motion and collisions etc. The 'Bouncing_Ball.py' file just contains the code that initialises the particles and draws them. The 'graphics.py' module is just one I found on line which uses tKinter to draw basic shapes. (Unfortunately I had to convert them to text files in order to upload them).

Thanks so much for any help in advance :)

 

[DrClaude]

Is the total energy increasing as a function of time or is it stable on averge? If it is the latter, then there is nothing to worry about: I've never seen such a simulation where this is not the case. Discretization of time is the culprit.

 

[JizzaDaMan]

It is stable on average, but it's not just changing by the odd joule it's in the order of a kJ the first time. Once it happens once, then it begins to increase at an increasing rate. I won't say exponentially though because it is still random. When I change the mass of one molecule then the kinetic energy increases very quickly, and the rate is somewhere in between when I try to implement gravity.

 

[DrClaude]

What happens when you reduce the time step?

 

[JizzaDaMan]

It reduces the frequency of errors

 

[DrClaude]

I implemented that by following this semi-tutorial thingy:
http://freespace.virgin.net/hugo.elias/models/m_snokr.htm

I looked at that link, and the equations are incorrect. I don't understand why he takes the square root of the "impact speed", it doesn't make sense in terms of units. And $\sqrt{m_1 m_2}$ also doesn't make any sense, as you would rather expect the reduced mass $m_1 m_2 / (m_1 + m_2)$ to appear in there.

I looked on the internet, and found this page: http://introcs.cs.princeton.edu/java/assignments/collisions.html

which explains things correctly. In particular, scroll down to the section entitled "Collision resolution."

 

[JizzaDaMan]

ok thanks I'll have a look at that :)