CODEWITHSUNDEEP
julialinear-algebracollisionparticles
Struggling_Student
Struggling_Student

Reputation: 464

simulating the collision of particles in Julia

I would like to simulate the collision of particles inside a box.
To be more specific I want to create a function (lets call it collision!), that updates the particles velocities after each interaction, like shown in the image.

enter image description here

I defined the particles (with radius equal 1) as followed:

mutable struct Particle
    pos :: Vector{Float64}
    vel :: Vector{Float64}
end

p = Particle( rand(2) , rand(2) )
# example for the position
p.pos
> 2-element Vector{Float64}:
  0.49339012018408135
  0.11441734325871078

And for the collision

function collision!(p1::Particle, p2::Particle)
     
    # ... #
    
    return nothing
end

The main idea is that when two particles collide, they "exchange" their velocity vector that is parallel to the particles centers (vector n hat).
In order to do that, one would need to transform the velocity vectors to the orthonormal basis of the collision normal (n hat). enter image description here

Then exchange the parallel component and rotate it in the original basis back.

enter image description here

I think I got the math right but I am not sure how to implement it in the code

Upvotes: 2

Views: 274

Answers (1)

cbk
cbk

Reputation: 4370

With the caveat that I have not checked the math at all, one implementation for the 2d case you provide might be along the lines of:

struct Particle
    pos :: Vector{Float64}
    vel :: Vector{Float64}
end

p1 = Particle( rand(2) , rand(2) )
p2 = Particle( rand(2) , rand(2) )

function collision!(p1::Particle, p2::Particle)
    # Find collision vector
    n = p1.pos - p2.pos
    # Normalize it, since you want an orthonormal basis
    n ./= sqrt(n[1]^2 + n[2]^2)
    # Construct M
    M = [n[1] n[2]; -n[2] n[1]]
    # Find transformed velocity vectors
    v1ₙ = M*p1.vel
    v2ₙ = M*p2.vel
    # Swap first component (or should it be second? Depends on how M was constructed)
    v1ₙ[1], v2ₙ[1] = v2ₙ[1], v1ₙ[1]
    # Calculate and store new velocity vectors
    p1.vel .= M'*v1ₙ
    p2.vel .= M'*v2ₙ
    return nothing
end

A few points:

  1. You don't need a mutable struct; just a plain struct will work fine since the Vector itself is mutable
  2. This implementation has a lot of excess allocations that you could avoid if you could work either in-place or perhaps more feasibly on the stack (for example, using StaticArrays of some sort instead of base Arrays as the basis for your position and velocity vectors). In-place actually might not be too hard either if you just make another struct (say "CollisionEvent") which holds preallocated buffers for M, n, v1n and v2n, and pass that to the collision! function as well.
  3. While I have not dived in to see, one might be able to find useful reference implementations for this type of collision in a molecular dynamics package like https://github.com/JuliaMolSim/Molly.jl

Upvotes: 2

Related Questions