Brownian motion - Thermoproject IP2

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):
        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

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')) 

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()

We gaan nu de functies van de simulatie weer aanroepen:

# 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

collisions = 0
collision_history = []

def handle_collisions(particles):
#your code/answer
    """ alle onderlinge botsingen afhandelen voor deeltjes in lijst """
    global collisions
    collisions = 0
    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]):
                particle_collision(particles[i], particles[j])
                collisions +=1
#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 = []

#tracken licht deeltje
track_xL = []
track_yL = []

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

#licht deeltje track list appenden
    track_xL.append(particles[N-2].r[0])
    track_yL.append(particles[N-2].r[1])
    
#zwaar deeltje track list appenden    
    track_x.append(particles[N-1].r[0])
    track_y.append(particles[N-1].r[1])
    
    collision_history.append(collisions)


plt.figure(figsize=(12,12))
plt.xlabel("x")
plt.ylabel("y")
plt.plot(track_x,track_y,'r',label="zwaar deeltje")
plt.plot(track_xL,track_yL,'b-',label="licht deeltje")
plt.legend()
plt.title("Afgelegde pad licht & zwaar deeltje")
plt.show()
print(collisions)

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
plt.figure(figsize=(10,6))
plt.plot(collision_history)
plt.xlabel("Tijdstap")
plt.ylabel("Totaal aantal botsingen")
plt.title("Totaal aantal botsingen tijdens simulatie")
plt.show()

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.

# Particle class

class ParticleClass:
    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 = R
        self.c = c

    def update_position(self):
        self.r += self.v * dt

    def boxcollision(self):
        if abs(self.r[0]) + self.R > Box_length:
            self.v[0] = -self.v[0]
            self.r[0] = np.sign(self.r[0]) * (Box_length - self.R)

        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 0.5 * self.m * np.dot(self.v, self.v)



# Collision functions

def collide_detection(p1, p2):
    dx = p1.r[0] - p2.r[0]
    dy = p1.r[1] - p2.r[1]
    rr = p1.R + p2.R
    return dx**2 + dy**2 < rr**2


def particle_collision(p1, p2):
    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)
    if dot_product >= 0:
        return

    distance_squared = np.dot(delta_r, delta_r)
    p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
    p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_r


def handle_collisions(particles):
    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]):
                particle_collision(particles[i], particles[j])



# Simulation setup


N = 361
N_light = N - 1
v_0 = 1
dt = 0.04

# Box schaalt met aantal deeltjes
density = 1e-2                     # aanpasbare "dichtheid" voor variabele box maat
Box_size_0 = np.sqrt(N / density)
Box_length_0 = Box_size_0 / 2
Box_length = Box_length_0

particles = []

# lichte deeltjes
for i in range(N_light):
    vx = np.random.uniform(-v_0, v_0)
    vy = np.random.choice([-1, 1]) * np.sqrt(max(v_0**2 - vx**2, 0))
    pos = Box_length * np.random.uniform(-1, 1, 2)
    particles.append(ParticleClass(m=1.0, v=[vx, vy], r=pos, R=0.5, c='blue'))

# zwaar deeltje
particles.append(ParticleClass(m=25.0, v=[0, 0], r=[0, 0], R=1.0, c='red'))


# Tracking afgelegde weg

distance_travelled = np.zeros(N)
previous_positions = np.array([p.r.copy() for p in particles])

steps = 500


# Simulatie loop

for step in range(steps):

    for idx, p in enumerate(particles):

        p.update_position()
        p.boxcollision()

        dr = np.linalg.norm(p.r - previous_positions[idx])
        distance_travelled[idx] += dr

        previous_positions[idx] = p.r.copy()

    handle_collisions(particles)



# Plot histogram

plt.figure(figsize=(12, 12))
plt.hist(distance_travelled[0:], bins=40, label='Normale deeltjes')
plt.hist(distance_travelled[-1], bins=7, label='Zwaar deeltje', color='red')


plt.xlabel("Afgelegde weg")
plt.ylabel("Aantal deeltjes")
plt.title("Histogram van afgelegde weg van alle deeltjes")
plt.legend()
plt.show()

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”.


# Particle class

class ParticleClass:
    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 = R
        self.c = c

    def update_position(self):
        global g_vector
        # positie-update met zwaartekracht
        self.v += g_vector * dt
        self.r += self.v * dt

    def boxcollision(self):
        # botsing met de verticale muren
        if abs(self.r[0]) + self.R > Box_width:
            self.v[0] = -self.v[0]
            self.r[0] = np.sign(self.r[0]) * (Box_width - self.R)

        # botsing met de horizontale muren
        if abs(self.r[1]) + self.R > Box_height:
            self.v[1] = -self.v[1]
            self.r[1] = np.sign(self.r[1]) * (Box_height - self.R)


def collide_detection(p1, p2):
    dx = p1.r[0] - p2.r[0]
    dy = p1.r[1] - p2.r[1]
    rr = p1.R + p2.R
    return dx**2 + dy**2 < rr**2


def particle_collision(p1, p2):
    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)
    if dot_product >= 0:
        return

    dist_sq = np.dot(delta_r, delta_r)
    p1.v -= 2 * m2 / (m1 + m2) * dot_product / dist_sq * delta_r
    p2.v += 2 * m1 / (m1 + m2) * dot_product / dist_sq * delta_r


def handle_collisions(particles):
    Np = len(particles)
    for i in range(Np):
        for j in range(i + 1, Np):
            if collide_detection(particles[i], particles[j]):
                particle_collision(particles[i], particles[j])



# Simulation parameters


# totaal aantal deeltjes *verdubbeld*
N = 722                       
N_half = N // 2               

v_0 = 1
dt = 0.02

Box_width = 40
Box_height = 80

# sterke artificiële zwaartekracht
g = -5.0                      
g_vector = np.array([0, g])

particles = []


# Twee groepen maken


# Groep lichte deeltjes
for i in range(N_half):
    vx = np.random.uniform(-v_0, v_0)
    vy = np.random.uniform(-v_0, v_0)
    pos = np.array([np.random.uniform(-Box_width, Box_width),
                    np.random.uniform(-Box_height, Box_height)])
    particles.append(ParticleClass(m=1.0, v=[vx, vy], r=pos, R=0.4, c='blue'))

# Groep zware deeltjes
for i in range(N_half):
    vx = np.random.uniform(-v_0, v_0)
    vy = np.random.uniform(-v_0, v_0)
    pos = np.array([np.random.uniform(-Box_width, Box_width),
                    np.random.uniform(-Box_height, Box_height)])
    particles.append(ParticleClass(m=5.0, v=[vx, vy], r=pos, R=0.4, c='red'))



# Simulatie loop

steps = 400

for step in range(steps):

    for p in particles:
        p.update_position()
        p.boxcollision()

    handle_collisions(particles)


# Plot eindconfiguratie

plt.figure(figsize=(8,16))
plt.title("Verwachte scheiding: lichte deeltjes 'drijven' omhoog")
plt.xlabel("x")
plt.ylabel("y")
plt.xlim(-Box_width, Box_width)
plt.ylim(-Box_height, Box_height)

for p in particles:
    plt.plot(p.r[0], p.r[1], '.', color=p.c, markersize=4)

plt.grid(True)
plt.show()