In this demo, we will find how the speed and 1-D velocity distributions predicted by kinetic gas theory can be found without any "derivation", using the simulation below. First, we will import a module containing the classes for this simulation, CollidingParticles, along with some modules for plotting:
#NAME: 1D Velocity Distributions
#DESCRIPTION: Finding the 1D velocity distributions predicted by kinetic gas theory through simulation of hard spheres.
from pycav.mechanics import *
from vpython import *
import numpy as np
import matplotlib.pyplot as plt
Next, we will run the simulation itself. First, we create a canvas which the simulation will display in. We then create a 1x1x1 container, then create a system composed of 200 hard sphere particles inside that container. The simulation will take a minute or two to finish. Once it is done, execute the last two cells. If it takes too long, try running the simulation with the visualisation off, by intialising system with visualize = False, or running the simulation for less time, or having less particles.
container = Container(1)
system = System(collides = True, interacts = False, visualize = True, container = container)
avg_speed = 1
system.create_particles_in_container(number = 100, speed = avg_speed, radius = 0.03)
dt = 1E-2
while system.steps <= 1000:
rate(150)
system.simulate(dt)
--------------------------------------------------------------------------- KeyboardInterrupt Traceback (most recent call last) <ipython-input-2-ee98fc628c2a> in <module>() 6 while system.steps <= 1000: 7 rate(150) ----> 8 system.simulate(dt) /Applications/anaconda/lib/python3.5/site-packages/pycav/mechanics.py in simulate(self, dt) 854 # Collision detection between particles using domains if has container. 855 if self.collides: --> 856 self._collision_detection_with_domains() 857 858 # Collision detection with walls of container if has one. /Applications/anaconda/lib/python3.5/site-packages/pycav/mechanics.py in _collision_detection_with_domains(self) 1064 # Should be fine as if particles are moving so quickly, simulation is inacurrate anyway.) 1065 if self.steps % 3 == 0: -> 1066 self._assign_particles_to_domains() 1067 1068 # The actual collision detection /Applications/anaconda/lib/python3.5/site-packages/pycav/mechanics.py in _assign_particles_to_domains(self) 923 domain.particles = [] 924 for particle in self.particles: --> 925 if domain.contains(particle) is True: 926 domain.particles.append(particle) 927 /Applications/anaconda/lib/python3.5/site-packages/pycav/mechanics.py in contains(self, particle) 378 Particle which is being checked to see if inside container or not 379 """ --> 380 relative_pos = particle.pos - self._origin_pos 381 for i in range(0, 3): 382 if relative_pos[i] > self.dimension or relative_pos[i] < 0: /Applications/anaconda/lib/python3.5/site-packages/pycav/mechanics.py in _origin_pos(self) 304 """ 305 l2 = self.dimension / 2 --> 306 return (self.pos - np.array([l2, l2, l2])) 307 308 @property KeyboardInterrupt:
We expect a 1D velocity distribution centred on 0, with the following form.
$$ \\f(v_{x}) = \sqrt{\frac{m}{2\pi k_{B}T}} e^{\frac{-mv_{x}^{2}}{2k_{B}T}}$$Assuming the relation
$$ \\<v> = \sqrt{\frac{8k_{B}T}{m}} $$,where v denotes the speed in this case, we can find the expected distribution for the system's 1D speed distribution, and compare it with the actual measured one.
%matplotlib notebook
#alpha is a constant = 2kT/m
alpha = ((avg_speed)**2) * np.pi / (4)
@np.vectorize
def f(v):
return ((1/(alpha*np.pi))**(0.5)) * np.exp(-(v**2)/alpha)
plt.hist(system.one_d_velocities, 20, normed = True, facecolor='green')
v = np.linspace(np.amin(system.one_d_velocities),np.amax(system.one_d_velocities),50)
l = plt.plot(v, f(v), 'r--', linewidth=1)
plt.ylabel('Probability density')
plt.xlabel('velocity')
plt.title('1D velocity distribution')
plt.show()
We expect a 1D speed distribution centred on 0, with the following form.
$$ \\f(v_{x}) = \frac{m}{2\pi k_{B}T}^{\frac{3}{2}} e^{\frac{-mv^{2}}{2k_{B}T}}$$Assuming the relation
$$ \\<v> = \sqrt{\frac{8k_{B}T}{m}} $$again, we can find the expected distribution for the system's speed distribution, and compare it with the actual measured one.
%matplotlib notebook
#alpha is a constant = 2kT/m
alpha = ((avg_speed)**2) * np.pi / (4)
@np.vectorize
def f(v):
return ((1/(alpha*np.pi))**(1.5)) * (4*np.pi*((v)**2)) * np.exp(-(v**2)/alpha)
values, bins, __ = plt.hist(system.speeds, 20, normed = True, facecolor='green')
v = np.linspace(np.amin(system.speeds),np.amax(system.speeds),20)
l = plt.plot(v, f(v), 'r--', linewidth=1)
plt.ylabel('Probability density')
plt.xlabel('velocity')
plt.title('3D speed distribution')
plt.show()