# Collisions

No preview image

### 1 collaborator Michael Novak (Author)

### Tags

(This model has yet to be categorized with any tags)
Visible to everyone | Changeable by the author
Model was written in NetLogo 6.0.4 • Viewed 114 times • Downloaded 4 times • Run 0 times Download this modelEmbed this model

## WHAT IS IT?

This model is one in a series of GasLab models. They use the same basic rules for simulating the behavior of gases. Each model integrates different features in order to highlight different aspects of gas behavior. This model is different from the other GasLab models in that the collision calculations take the circular shape and size of the particles into account, instead of modeling the particles as dimensionless points.

## HOW IT WORKS

The model determines the resulting motion of particles that collide, with no loss in their total momentum or total kinetic energy (an elastic collision).

To calculate the outcome of collision, it is necessary to calculate the exact time at which the edge of one particle (represented as a circle), would touch the edge of another particle (or the walls of a container) if the particles were allowed to continue with their current headings and speeds.

By performing such a calculation, one can determine when the next collision anywhere in the system would occur in time. From this determination, the model then advances the motion of all the particles using their current headings and speeds that far in time until this next collision point is reached. Exchange of kinetic energy and momentum between the two particles, according to conservation of kinetic energy and conservation of momentum along the collision axis (a line drawn between the centers of the two particles), is then calculated, and the particles are given new headings and speeds based on this outcome.

## HOW TO USE IT

INITIAL-NUMBER-PARTICLES determines the number of gas particles used with SETUP. If the world is too small or the particles are too large, the SETUP procedure of the particles will stop so as to prevent overlapping particles.

SMALLEST-PARTICLE-SIZE and LARGEST-PARTICLE-SIZE determines the range of particle sizes that will be created when SETUP is pressed. (Particles are also assigned a mass proportional to the area of the particle that is created.)

The SETUP button will set the initial conditions. The GO button will run the simulation.

Monitors:

• % FAST, % MEDIUM, % SLOW: the percentage of particles with different speeds: fast (red), medium (green), and slow (blue).
• AVERAGE SPEED: average speed of the particles.
• AVERAGE ENERGY: average kinetic energy of the particles.

Plots:

• SPEED COUNTS: plots the number of particles in each range of speed.
• SPEED HISTOGRAM: speed distribution of all the particles. The gray line is the average value, and the black line is the initial average.
• ENERGY HISTOGRAM: distribution of energies of all the particles, calculated as m*(v^2)/2. The gray line is the average value, and the black line is the initial average.

Initially, all the particles have the same speed but random directions. Therefore the first histogram plots of speed will show only one column. If all the particles have the same size (and therefore the same mass), then the first histogram plot of energy will also show one column. As the particles repeatedly collide, they exchange energy and head off in new directions, and the speeds are dispersed -- some particles get faster, some get slower. The histogram distribution changes accordingly.

## THINGS TO NOTICE

Run the model with different masses, by setting the MAX-PARTICLE-SIZE larger than the MIN-PARTICLE-SIZE. Mass is scaled linearly to the area of the particle, so that a particle that is twice the radius of another particle, has four the area and therefore four times the mass.

With many different mass particles colliding over time, different sized particles start to move at different speed ranges (in general). The smallest mass particles will be usually moving faster (red) than the average particle speed and the largest mass particles will be usually slower (blue) than the average particle speed. This emergent result is what happens in a gas that is a mixture of particles of different masses. At any given temperature, the higher mass particles are moving slower (such as Nitrogen gas: N2) then the lower mass particles (such as water vapor: H2O).

The particle histograms quickly converge on the classic Maxwell-Boltzmann distribution. What's special about these curves? Why is the shape of the energy curve not the same as the speed curve?

Look at the other GasLab models to compare the particle histograms. In those models, the particles are modeled as "point" particles, with an area of "zero". The results of those models should match this model, even with that simplification.

With particles of different sizes, you may notice some fast moving particles have lower energy than medium speed particles. How can the difference in the mass of the particles account for this?

## THINGS TO TRY

Setting all the particles to have a very slow speed (e.g. 0.001) and one particle to have a very fast speed helps show how kinetic energy is eventually transferred to all the particles through a series of collisions and would serve as a good model for energy exchange through conduction between hot and cold gases.

To see what the approximate mass of each particle is, type this in the command center:

``````ask particles [ set label precision mass 0 ]
``````

## EXTENDING THE MODEL

Collisions between boxes and circles could also be explored. Variations in size between particles could investigated or variations in the mass of some of the particle could be made to explore other factors that affect the outcome of collisions.

## NETLOGO FEATURES

Instead of advancing one tick at a time as in most models, the tick counter takes on fractional values, using the `tick-advance` primitive. (In the Interface tab, it is displayed as an integer, but if you make a monitor for `ticks` you'll see the exact value.)

## RELATED MODELS

Look at the other GasLab models to see collisions of "point" particles, that is, the particles are assumed to have an area or volume of zero.

## CREDITS AND REFERENCES

This model was developed as part of the GasLab curriculum (http://ccl.northwestern.edu/curriculum/gaslab/) and has also been incorporated into the Connected Chemistry curriculum (http://ccl.northwestern.edu/curriculum/ConnectedChemistry/)

## HOW TO CITE

If you mention this model or the NetLogo software in a publication, we ask that you include the citations below.

For the model itself:

Please cite the NetLogo software as: This model and associated activities and materials were created as part of the project: MODELING ACROSS THE CURRICULUM. The project gratefully acknowledges the support of the National Science Foundation, the National Institute of Health, and the Department of Education (IERI program) -- grant number REC #0115699.

Click to Run Model

```globals [
tick-delta        ;; how much simulation time will pass in this step
box-edge          ;; distance of box edge from origin
collisions        ;; list used to keep track of future collisions
particle1         ;; first particle currently colliding
particle2         ;; second particle currently colliding
init-avg-speed init-avg-energy             ;; initial averages
avg-speed avg-energy                       ;; current averages
fast medium slow                           ;; current counts
percent-slow percent-medium percent-fast   ;; percentage of current counts
]

breed [particles particle]

particles-own [
speed
mass
energy
]

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;; setup procedures ;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to setup
clear-all
set-default-shape particles "circle"
set box-edge max-pxcor - 1
ask patches with [(abs pxcor = box-edge or abs pycor = box-edge) and
abs pxcor <= box-edge and abs pycor <= box-edge]
[ set pcolor yellow ]
set avg-speed 1
make-particles
set particle1 nobody
set particle2 nobody
reset-ticks
set collisions []
update-variables
set init-avg-speed avg-speed
set init-avg-energy avg-energy
end

to make-particles
create-particles initial-number-particles [
set speed 1

set size smallest-particle-size
+ random-float (largest-particle-size - smallest-particle-size)
;; set the mass proportional to the area of the particle
set mass (size * size)
set energy kinetic-energy

recolor
]
;; When space is tight, placing the big particles first improves
;; our chances of eventually finding places for all of them.
foreach sort-by [ [a b] -> [ size ] of a > [ size ] of b ] particles [ the-particle ->
position-randomly
while [ overlapping? ] [ position-randomly ]
]
]
end

to-report overlapping?  ;; particle procedure
;; here, we use IN-RADIUS just for improved speed; the real testing
;; is done by DISTANCE
report any? other particles in-radius ((size + largest-particle-size) / 2)
with [distance myself < (size + [size] of myself) / 2]
end

to position-randomly  ;; particle procedure
;; place particle at random location inside the box
setxy one-of [1 -1] * random-float (box-edge - 0.5 - size / 2)
one-of [1 -1] * random-float (box-edge - 0.5 - size / 2)
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;; go procedures  ;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to go
choose-next-collision
ask particles [ jump speed * tick-delta ]
perform-next-collision
recalculate-particles-that-just-collided
if floor ticks > floor (ticks - tick-delta)
[
update-variables
update-plots
]
end

to update-variables
set medium count particles with [color = green]
set slow   count particles with [color = blue]
set fast   count particles with [color = red]
set percent-medium (medium / ( count particles )) * 100
set percent-slow (slow / (count particles)) * 100
set percent-fast (fast / (count particles)) * 100
set avg-speed  mean [speed] of particles
set avg-energy mean [energy] of particles
end

to recalculate-particles-that-just-collided
;; Since only collisions involving the particles that collided most recently could be affected,
;; we filter those out of collisions.  Then we recalculate all possible collisions for
;; the particles that collided last.  The ifelse statement is necessary because
;; particle2 can be either a particle or a string representing a wall.  If it is a
;; wall, we don't want to invalidate all collisions involving that wall (because the wall's
;; position wasn't affected, those collisions are still valid.
ifelse is-turtle? particle2
[
set collisions filter [ the-collision ->
item 1 the-collision != particle1 and
item 2 the-collision != particle1 and
item 1 the-collision != particle2 and
item 2 the-collision != particle2
] collisions
]
[
set collisions filter [ the-collision ->
item 1 the-collision != particle1 and
item 2 the-collision != particle1
] collisions
]
if particle1 != nobody [ ask particle1 [ check-for-wall-collision ] ]
if particle1 != nobody [ ask particle1 [ check-for-particle-collision ] ]
;; Slight errors in floating point math can cause a collision that just
;; happened to be calculated as happening again a very tiny amount of
;; time into the future, so we remove any collisions that involves
;; the same two particles (or particle and wall) as last time.
set collisions filter [ the-collision ->
item 1 the-collision != particle1 or
item 2 the-collision != particle2
] collisions
;; All done.
set particle1 nobody
set particle2 nobody
end

;; check-for-particle-collision is a particle procedure that determines the time it takes
;; to the collision between two particles (if one exists).  It solves for the time by representing
;; the equations of motion for distance, velocity, and time in a quadratic equation of the vector
;; components of the relative velocities and changes in position between the two particles and
;; solves for the time until the next collision

to check-for-particle-collision
let my-x xcor
let my-y ycor
let my-particle-size size
let my-x-speed speed * dx
let my-y-speed speed * dy
[
let dpx (xcor - my-x)   ;; relative distance between particles in the x direction
let dpy (ycor - my-y)    ;; relative distance between particles in the y direction
let x-speed (speed * dx) ;; speed of other particle in the x direction
let y-speed (speed * dy) ;; speed of other particle in the x direction
let dvx (x-speed - my-x-speed) ;; relative speed difference between particles in x direction
let dvy (y-speed - my-y-speed) ;; relative speed difference between particles in y direction
let sum-r (((my-particle-size) / 2 ) + (([size] of self) / 2 )) ;; sum of both particle radii

;; To figure out what the difference in position (P1) between two particles at a future
;; time (t) will be, one would need to know the current difference in position (P0) between the
;; two particles and the current difference in the velocity (V0) between the two particles.
;;
;; The equation that represents the relationship is:
;;   P1 = P0 + t * V0
;; we want find when in time (t), P1 would be equal to the sum of both the particle's radii
;; (sum-r).  When P1 is equal to is equal to sum-r, the particles will just be touching each
;; other at their edges (a single point of contact).
;;
;; Therefore we are looking for when:   sum-r =  P0 + t * V0
;;
;; This equation is not a simple linear equation, since P0 and V0 should both have x and y
;; components in their two dimensional vector representation (calculated as dpx, dpy, and
;; dvx, dvy).
;;
;; By squaring both sides of the equation, we get:
;;   (sum-r) * (sum-r) =  (P0 + t * V0) * (P0 + t * V0)
;; When expanded gives:
;;   (sum-r ^ 2) = (P0 ^ 2) + (t * PO * V0) + (t * PO * V0) + (t ^ 2 * VO ^ 2)
;; Which can be simplified to:
;;   0 = (P0 ^ 2) - (sum-r ^ 2) + (2 * PO * V0) * t + (VO ^ 2) * t ^ 2
;; Below, we will let p-squared represent:   (P0 ^ 2) - (sum-r ^ 2)
;; and pv represent: (2 * PO * V0)
;; and v-squared represent: (VO ^ 2)
;;
;;  then the equation will simplify to:     0 = p-squared + pv * t + v-squared * t^2

let p-squared   ((dpx * dpx) + (dpy * dpy)) - (sum-r ^ 2)   ;; p-squared represents difference
;; of the square of the radii and the square of the initial positions

let pv  (2 * ((dpx * dvx) + (dpy * dvy)))  ;; vector product of the position times the velocity
let v-squared  ((dvx * dvx) + (dvy * dvy)) ;; the square of the difference in speeds
;; represented as the sum of the squares of the x-component
;; and y-component of relative speeds between the two particles

;; p-squared, pv, and v-squared are coefficients in the quadratic equation shown above that
;; represents how distance between the particles and relative velocity are related to the time,
;; t, at which they will next collide (or when their edges will just be touching)

;; Any quadratic equation that is a function of time (t) can be represented as:
;;   a*t*t + b*t + c = 0,
;; where a, b, and c are the coefficients of the three different terms, and has solutions for t
;; that can be found by using the quadratic formula.  The quadratic formula states that if a is
;; not 0, then there are two solutions for t, either real or complex.
;; t is equal to (b +/- sqrt (b^2 - 4*a*c)) / 2*a
;; the portion of this equation that is under a square root is referred to here
;; as the determinant, d1.   d1 is equal to (b^2 - 4*a*c)
;; and:   a = v-squared, b = pv, and c = p-squared.
let d1 pv ^ 2 -  (4 * v-squared * p-squared)

;; the next test tells us that a collision will happen in the future if
;; the determinant, d1 is > 0,  since a positive determinant tells us that there is a
;; real solution for the quadratic equation.  Quadratic equations can have solutions
;; that are not real (they are square roots of negative numbers).  These are referred
;; to as imaginary numbers and for many real world systems that the equations represent
;; are not real world states the system can actually end up in.

;; Once we determine that a real solution exists, we want to take only one of the two
;; possible solutions to the quadratic equation, namely the smaller of the two the solutions:
;;  (b - sqrt (b^2 - 4*a*c)) / 2*a
;;  which is a solution that represents when the particles first touching on their edges.
;;  instead of (b + sqrt (b^2 - 4*a*c)) / 2*a
;;  which is a solution that represents a time after the particles have penetrated
;;  and are coming back out of each other and when they are just touching on their edges.

let time-to-collision  -1

if d1 > 0
[ set time-to-collision (- pv - sqrt d1) / (2 * v-squared) ]        ;; solution for time step

;; if time-to-collision is still -1 there is no collision in the future - no valid solution
;; note:  negative values for time-to-collision represent where particles would collide
;; if allowed to move backward in time.
;; if time-to-collision is greater than 1, then we continue to advance the motion
;; of the particles along their current trajectories.  They do not collide yet.

if time-to-collision > 0
[
;; time-to-collision is relative (ie, a collision will occur one second from now)
;; We need to store the absolute time (ie, a collision will occur at time 48.5 seconds.
;; So, we add clock to time-to-collision when we store it.
;; The entry we add is a three element list of the time to collision and the colliding pair.
set collisions fput (list (time-to-collision + ticks) self myself)
collisions
]
]
end

;; determines when a particle will hit any of the four walls

to check-for-wall-collision  ;; particle procedure
;; right & left walls
let x-speed (speed * dx)
if x-speed != 0
[ ;; solve for how long it will take particle to reach right wall
let right-interval (box-edge - 0.5 - xcor - size / 2) / x-speed
if right-interval > 0
[ assign-colliding-wall right-interval "right wall" ]
;; solve for time it will take particle to reach left wall
let left-interval ((- box-edge) + 0.5 - xcor + size / 2) / x-speed
if left-interval > 0
[ assign-colliding-wall left-interval "left wall" ] ]
;; top & bottom walls
let y-speed (speed * dy)
if y-speed != 0
[ ;; solve for time it will take particle to reach top wall
let top-interval (box-edge - 0.5 - ycor - size / 2) / y-speed
if top-interval > 0
[ assign-colliding-wall top-interval "top wall" ]
;; solve for time it will take particle to reach bottom wall
let bottom-interval ((- box-edge) + 0.5 - ycor + size / 2) / y-speed
if bottom-interval > 0
[ assign-colliding-wall bottom-interval "bottom wall" ] ]
end

to assign-colliding-wall [time-to-collision wall]  ;; particle procedure
;; this procedure is used by the check-for-wall-collision procedure
;; to assemble the correct particle-wall pair
;; time-to-collision is relative (ie, a collision will occur one second from now)
;; We need to store the absolute time (ie, a collision will occur at time 48.5 seconds.
;; So, we add clock to time-to-collision when we store it.
let colliding-pair (list (time-to-collision + ticks) self wall)
set collisions fput colliding-pair collisions
end

to choose-next-collision
if collisions = [] [ stop ]
;; Sort the list of projected collisions between all the particles into an ordered list.
;; Take the smallest time-step from the list (which represents the next collision that will
;; happen in time).  Use this time step as the tick-delta for all the particles to move through
let winner first collisions
foreach collisions [ the-collision ->
if first the-collision < first winner [ set winner the-collision ]
]
;; winner is now the collision that will occur next
let dt item 0 winner
;; If the next collision is more than 1 in the future,
;; only advance the simulation one tick, for smoother animation.
set tick-delta dt - ticks
if tick-delta > 1
[ set tick-delta 1
set particle1 nobody
set particle2 nobody
stop ]
set particle1 item 1 winner
set particle2 item 2 winner
end

to perform-next-collision
;; deal with 3 possible cases:
;; 1) no collision at all
if particle1 = nobody [ stop ]
;; 2) particle meets wall
if is-string? particle2
[ if particle2 = "left wall" or particle2 = "right wall"
stop ]
if particle2 = "top wall" or particle2 = "bottom wall"
stop ] ]
;; 3) particle meets particle
ask particle1 [ collide-with particle2 ]
end

to collide-with [other-particle]  ;; particle procedure
;;; PHASE 1: initial setup
;; for convenience, grab some quantities from other-particle
let mass2 [mass] of other-particle
let speed2 [speed] of other-particle
;; modified so that theta is heading toward other particle
let theta towards other-particle

;;; PHASE 2: convert velocities to theta-based vector representation
;; now convert my velocity from speed/heading representation to components
;; along theta and perpendicular to theta
let v1t (speed * cos (theta - heading))
let v1l (speed * sin (theta - heading))
;; do the same for other-particle
let v2t (speed2 * cos (theta - heading2))
let v2l (speed2 * sin (theta - heading2))

;;; PHASE 3: manipulate vectors to implement collision
;; compute the velocity of the system's center of mass along theta
let vcm (((mass * v1t) + (mass2 * v2t)) / (mass + mass2) )
;; now compute the new velocity for each particle along direction theta.
;; velocity perpendicular to theta is unaffected by a collision along theta,
;; so the next two lines actually implement the collision itself, in the
;; sense that the effects of the collision are exactly the following changes
;; in particle velocity.
set v1t (2 * vcm - v1t)
set v2t (2 * vcm - v2t)

;;; PHASE 4: convert back to normal speed/heading
;; now convert my velocity vector into my new speed and heading
set speed sqrt ((v1t * v1t) + (v1l * v1l))
;; if the magnitude of the velocity vector is 0, atan is undefined. but
;; speed will be 0, so heading is irrelevant anyway. therefore, in that
;; case we'll just leave it unmodified.
set energy kinetic-energy

if v1l != 0 or v1t != 0
[ set heading (theta - (atan v1l v1t)) ]
;; and do the same for other-particle
set speed sqrt ((v2t ^ 2) + (v2l ^ 2))
set energy kinetic-energy

if v2l != 0 or v2t != 0
[ set heading (theta - (atan v2l v2t)) ]
]

;; PHASE 5: recolor
;; since color is based on quantities that may have changed
recolor
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;; particle coloring procedures ;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to recolor ;; particle procedure
;;let avg-speed 1
;; avg-speed is assumed to be 0.5, since particles are assigned a random speed between 0 and 1
;; particle coloring procedures for visualizing speed with a color palette,
;; red are fast particles, blue slow, and green in between.
ifelse speed < (0.5 * avg-speed) ;; at lower than 50% the average speed
[ set color blue ]      ;; slow particles colored blue
[ ifelse speed > (1.5 * avg-speed) ;; above 50% higher the average speed
[ set color red ]        ;; fast particles colored red
[ set color green ] ]    ;; medium speed particles colored green
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;; time reversal procedure  ;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Here's a procedure that demonstrates time-reversing the model.
;; You can run it from the command center.  When it finishes,
;; the final particle positions may be slightly different because
;; the amount of time that passes after the reversal might not
;; be exactly the same as the amount that passed before; this
;; doesn't indicate a bug in the model.
;; For larger values of n, you will start to notice larger
;; discrepancies, eventually causing the behavior of the system
;; to diverge totally. Unless the model has some bug we don't know
;; about, this is due to accumulating tiny inaccuracies in the
;; floating point calculations.  Once these inaccuracies accumulate
;; to the point that a collision is missed or an extra collision
;; happens, after that the reversed model will diverge rapidly.

to test-time-reversal [n]
setup
while [ticks < n] [ go ]
let old-ticks ticks
reverse-time
while [ticks < 2 * old-ticks] [ go ]
ask particles [ set color white ]
end

to reverse-time
ask particles [ rt 180 ]
set collisions []
;; the last collision that happened before the model was paused
;; (if the model was paused immediately after a collision)
;; won't happen again after time is reversed because of the
;; "don't do the same collision twice in a row" rule.  We could
;; try to fool that rule by setting particle1 and
;; particle2 to nobody, but that might not always work,
;; because the vagaries of floating point math means that the
;; collision might be calculated to be slightly in the past
;; (the past that used to be the future!) and be skipped.
;; So to be sure, we force the collision to happen:
perform-next-collision
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;; reporters ;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

to-report init-particle-speed
report 1
end

to-report max-particle-mass
report max [mass] of particles
end

to-report kinetic-energy
report (0.5 * mass * speed * speed)
end

to draw-vert-line [ xval ]
plotxy xval plot-y-min
plot-pen-down
plotxy xval plot-y-max
plot-pen-up
end