Brownian motion - Thermolab

Een van de bekendste voorbeelden van botsende deeltjes in de natuur is Brownian motion. Fijn gemalen pollen in water lijken te dansen in willekeurige richting. Dit komt doordat de pollen worden geraakt door watermoleculen die in alle richtingen bewegen. Omdat de pollen veel zwaarder zijn dan watermoleculen, dus de beweging van de pollen is veel langzamer en minder “intens” dan die van de watermoleculen. Dit proces van willekeurige beweging door botsingen met kleinere deeltjes wordt Brownian motion genoemd en kunnen we simuleren op basis van ons (premature) botsingsmodel. Daarbij kunnen we ook gebruik maken van de zojuist geleerde manier van tracking van deeltjes, waarbij we een zowel het zware bolletjes als een enkel deeltje kunnen volgen.

Let op! We bestuderen hier nog geen thermische effecten, deze opdrachten zijn met name bedoeld om beter te begrijpen hoe het botsingsmodel in elkaar zit.

import numpy as np
import matplotlib.pyplot as plt

# Maken van de ParticleClass

class ParticleClass:
    # Het maken van het deeltje
    def __init__(self, m, v, r, R, c):
        self.m = m                         
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = np.array(R, dtype=float)  
        self.c = c

    # Het updaten van de positie, eventueel met zwaartekracht
    def update_position(self, a):
        self.r += self.v * dt + 1/2 * a * dt**2  
              
    # Harde wand
    def boxcollision(self):
        if abs(self.r[0]) + self.R > Box_length: 
            self.v[0] = -self.v[0]                                  # Omdraaien van de snelheid
            self.r[0] = np.sign(self.r[0]) * (Box_length - self.R)  # Zet terug net binnen box                 
        if abs(self.r[1]) + self.R > Box_length: 
            self.v[1] = -self.v[1]     
            self.r[1] = np.sign(self.r[1]) * (Box_length - self.R) 
            
    @property
    def momentum(self):
        return self.m * self.v
    
    @property
    def kin_energy(self):
        return 1/2 * self.m * np.dot(self.v, self.v)

# Aanmaken van de randvoorwaarden en initiele condities
Box_size_0 = 10
Box_length_0 = Box_size_0/2
Box_length = Box_length_0     # De grootte van de box kan wijzigen!

# Particles
dt = 0.1
particles = []
N = 40
v_0 = 1
a = -9.81

dt = 0.04

# Aanmaken van deeltjes
for i in range(N-1):
    vx = np.random.uniform(-v_0,v_0)
    vy = np.random.choice([-1, 1])*np.sqrt(v_0**2-vx**2)        
    pos = Box_length_0*np.random.uniform(-1,1,2)
    particles.append(ParticleClass(m=1.0, v=[vx, vy], r = pos, R=.5,c='blue')) 

#Aanmaken van de pollen
particles.append(ParticleClass(m=20.0, v=[0, 0], r = [0, 0], R=.5,c='red'))

Er is een doos vol met deeltjes op willekeurige positie aangemaakt. We willen kijken waar de deeltjes zijn terechtgekomen. Hieronder staat dit weergegeven.

# Inspecteren van beginsituatie
plt.figure()

plt.xlabel('x')
plt.ylabel('y')

plt.xlim(-Box_length_0,Box_length_0)
plt.ylim(-Box_length_0,Box_length_0)


for particle, particle_object in enumerate(particles):
    plt.plot(particle_object.r[0],particle_object.r[1],color=particle_object.c,marker='.')
    plt.arrow(particle_object.r[0],particle_object.r[1], 
        particle_object.v[0],particle_object.v[1], 
        head_width=0.05, head_length=0.1, color='red')
plt.show()

# Het bepalen of er een botsing plaats vindt

def collide_detection(self, other):
    dx = self.r[0] - other.r[0]
    dy = self.r[1] - other.r[1]
    rr = self.R + other.R
    return  dx**2+dy**2 < rr**2 
        
def particle_collision(p1: ParticleClass, p2: ParticleClass):
    """ past snelheden aan uitgaande van overlap """
    m1, m2 = p1.m, p2.m
    delta_r = p1.r - p2.r
    delta_v = p1.v - p2.v
    dot_product = np.dot(delta_r, delta_v)
    
    # Als deeltjes van elkaar weg bewegen dan geen botsing
    if dot_product >= 0: # '='-teken voorkomt ook problemen als delta_r == \vec{0}
        return
    
    distance_squared = np.dot(delta_r, delta_r) 
    # Botsing oplossen volgens elastische botsing in 2D
    p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
    p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_r

count = 0
def handle_collisions(particles):
    global count
#your code/answer
    """ alle onderlinge botsingen afhandelen voor deeltjes in lijst """
    num_particles = len(particles)
    for i in range(num_particles):
        for j in range(i+1, num_particles):
            if collide_detection(particles[i], particles[j]):
                count+=1
                particle_collision(particles[i], particles[j])
#your code/answer

In onderstaande code geven we de code voor de simulatie en volgen we de positie van het zware deeltje.

# tracken van het zware deeltje
track_x = []
track_y = []


#your code/answer
track2_x = []
track2_y = []
#your code/answer

#your code/answer
for i in range(400):
#your code/answer
    
    for p in particles:
        p.update_position(0)    # Update positie        
        p.boxcollision()         # Wandbotsing werkt per deeltje
        
    handle_collisions(particles)

#your code/answer
    
    track_x.append(particles[N-1].r[0])
    track_y.append(particles[N-1].r[1])
    track2_x.append(particles[0].r[0])
    track2_y.append(particles[0].r[1])
    
#your code/answer

plt.figure()
plt.plot(track_x,track_y,'r')
plt.plot(track2_x, track2_y, 'b')
plt.xlabel('beweging x richting')
plt.ylabel('beweging y richting')
#your code/answer
plt.show()

We zouden gevoel willen krijgen voor het aantal botsingen dat per tijdseenheid plaatsvindt. Elke keer dat er een botsing plaatsvindt, zou de counter met 1 omhoog moeten gaan. Idealiter wordt het aantal botsingen opgeslagen in een array zodat je het aantal botsingen als functie van de tijd kunt weergeven.

#your code/answer
def take_time_step(particles):
    for p in particles:
        p.update_position(0)
        
    handle_collisions(particles)
botsingen = []
botsingen1 = []

for i in range(400):
    count = 0
    for p in particles:
        p.update_position(0)    # Update positie        
        p.boxcollision()         # Wandbotsing werkt per deeltje
        
    handle_collisions(particles)
    botsingen.append(count)
    botsingen1.append(sum(botsingen))


plt.figure()
plt.plot(botsingen1)
plt.xlabel('tijd')
plt.ylabel('aantal botsingen')
plt.show()
# deze deed het tot gister nog wel maar geeft nu hele rare getallen.....ik ben verdrietig......

In zulke fysica modellen is de afgelegde weg (afstand tussen begin en eindpunt) van belang. Deze afgelegde weg zegt iets over de snelheid van difussie. Idealiter bekijken we een histogram. Maar voor een histogram hebben we veel deeltjes nodig.

import time
newparticles = []

#1, 2
for i in range (360):
    vx = np.random.uniform(-v_0, v_0)
    vy = np.random.choice([-1, 1])*np.sqrt(v_0**2-vx**2)        
    pos = Box_length_0*np.random.uniform(-1,1,2)
    newparticles.append(ParticleClass(m=1.0, v=[vx, vy], r = pos, R=.5,c='blue'))

zwaardeeltje = ParticleClass(m=10.0, v=[vx, vy], r = pos, R = .5, c='red')
newparticles.append(zwaardeeltje)

for particle, particle_object in enumerate(newparticles):
    plt.plot(particle_object.r[0],particle_object.r[1],color=particle_object.c,marker='.')
    plt.arrow(particle_object.r[0],particle_object.r[1], 
        particle_object.v[0],particle_object.v[1], 
        head_width=0.05, head_length=0.1, color='green')
plt.show()

#3
# alle simulatie dingen

startpositie = np.array([])
for p in newparticles:
    startpositie = np.append(startpositie, p.r)
eindpositie = np.array([])

for p in newparticles:
    take_time_step(newparticles)
    eindpositie = np.append(eindpositie, p.r)
weg = np.abs(eindpositie - startpositie)
print('De afgelegde weg van het zware deeltje is', weg[-1])
plt.figure()
plt.xlabel('Afgelegde weg')
plt.ylabel('Counts')
plt.hist(weg, bins=40)
plt.axvline(weg[-1], linestyle='--', c='r')

De afgelegde weg van het zware deeltje is 11.9607838140045

En nu we toch bezig zijn met twee verschillende deeltjes....

We kunnen twee “groepen” van deeltjes aanmaken, elk met een andere massa. Als we dan de zwaartekracht aan zetten, dan zouden we verwachten dat de lichtere deeltjes boven komen “drijven”.

from matplotlib.animation import FuncAnimation

newboxsize = 20
newboxlength = 10
v0 = 1

particlesA = []
particlesB = []

for i in range(200):
    vx = np.random.uniform(-v0, v0)
    vy = np.random.choice([-1, 1])*np.sqrt(v0**2-vx**2)        
    pos = newboxlength*np.random.uniform(-1,1,2)
    particlesA.append(ParticleClass(m=1.0, v = [vx, vy], r = pos, R = .5, c='green'))


for i in range(200):
    vx = np.random.uniform(-v0, v0)
    vy = np.random.choice([-1, 1])*np.sqrt(v0**2-vx**2)        
    pos = newboxlength*np.random.uniform(-1,1,2)
    particlesB.append(ParticleClass(m=20.0, v = [vx, vy], r = pos, R = .5, c='blue'))



fig, ax = plt.subplots()
ax.set_xlim(-10, 10)
ax.set_ylim(-10, 10)
ax.set_aspect('equal')
ax.set_title("Particle Animation")
ax.set_xlabel("x")
ax.set_ylabel("y")

# Create the particle as a red dot
dot, = ax.plot([], [], 'ro', markersize=10);
dotA, = ax.plot([], [], 'go', ms=10)
dotB, = ax.plot([], [], 'bo', ms=10)



# Initialization function for animation
def init():
    dot.set_data([], [])
    return dot,

def take_time_step(particles):
    for p in particlesA:
        p.update_position(-20)
    for p in particlesB:
        p.update_position(-20)

def update(frame):
    take_time_step(particlesA)
    take_time_step(particlesB)
    dotA.set_data([p.r[0] for p in particlesA], [p.r[1] for p in particlesA])
    dotB.set_data([p.r[0] for p in particlesB], [p.r[1] for p in particlesB])
    return dot,

# Create animation
ani = FuncAnimation(fig, update, frames=range(200), init_func=init, blit=True, interval=50)

# For Jupyter notebook:
from IPython.display import HTML
HTML(ani.to_jshtml())

#het werkt niet

Loading...