Event Driven Hard Sphere Collision

Test with Granular Dynamics Simulation

  • No fluid flow
  • Force free, e.g, no gravity.
  • Elastic bouncing walls: $e=1, \mu=0$ on all 6 bounding walls.

Efficiency Improvement in Collision Detection

Test case

  • L=100 cubic box with elastic walls.
  • 1000 spheres with dp=1
  • location: kalmar:/mnt/UofA/sekidata/test/granularFlow/GranularGas
  • With simulation of t=10, there are 182 collisions happened.

Brute Force Double Loop

Each double loop takes about 0.08 seconds to finish. Total walltime 0.08*181=14.6 seconds.

Linked Cell Algorithm

Each double loop takes about 0.003 seconds to finish. Total walltime 0.003*181=0.5 seconds.

Cell Size time per collision search (seconds) Identical to Brute Force result?
200 0.1 Yes
20 0.02 Yes
10 0.0055 Yes
5 0.003 NO
2 0.0025 NO

Further Improvement

It is possible to further improve the algorithm but it does not change the fact we have improved from $O(N^2)$ of Brute Force to $O(N)$.
Examples are the use of tournament tree data structure. Refer to the reference for details.

Reference: Chapter 3 of Computational Granular Dynamics: Models and Algorithms by Thorsten Pöschel and Thomas Schwager

Molecular Mono Gas at Constant State

With inter-particle $e=1, \mu=0$, the kinetic energy and granular temperature are kept as constants.
Note the units. For example, for Helium

  • density:$g/dm^3=kg/m^3$
  • length: $A= 10^{-10}m$
  • time: $10^{-13} s$
  • velocity: $km/s=10^{3}m/s$
  • engery: $10^{-24}J$
  • mass: $10^{-30} kg$
  • pressure: $MPa=10^6 Pa$

With this, the Boltzmann constant is scaled by 1e-24 to have the computed temperature with a unit of Kelvin. In the code, $k_B=13.806504$.

Typical time step is O(10), ie, 1ps.

Use the following Matlab code to generate the 3d.particles input file for Helium.

function nobelGasParInput3D(filename)
R=1.3;%uni A=1e-10 m
V=4/3*pi*R^3;%unit 1e-30 m^3
m=4*1.660538782e-27;%4 Dalton,unit kg
rhop=m/(V*1e-30)%unit kg/m^3, NOT the gas density
T=300;%temperature, unit K
k=1.3806504e-23;%Boltzman constant, unit J/K
stdV=sqrt(k*T/m)/1e3%unit km/s
p3d=load(filename);
nSp=length(p3d(:,1));
fp=fopen('3d.particles','w');
fprintf(fp,'#nSp=%d\n',nSp);
fprintf(fp,'#rho_f=1, t=0\n');
fprintf(fp,'#1-id 2-radius 3-rho_p 4-fixed 5-x 6-y 7-z 8-u 9-v 10-z 11-w1 12-w2 13-w3 14-theta 15-phi\n');
rhoRand=round(rand(nSp,1));
for i=1:nSp
        fprintf(fp, '%d %f %f %d %f %f %f %f %f %f %f %f %f %f %f\n' ...
                ,i-1,R,rhop,0,p3d(i,1),p3d(i,2),p3d(i,3),stdV*randn(3,1),0,0,0,0,0);
end
fclose(fp);

For example, in the 1003 box T=300K.

  • 1000 Helium spheres
    • $\rho=6.64216 kg/m^3$
    • Computed T=305K
    • $P=\rho R T=4.2MPa$
    • Average pressure computed in 1ns: $P=4.5 \times 10^6 Pa$
  • 1000 Neon spheres
    • $\rho=33.2108 kg/m^3$
    • Computed T=305.546 K
    • $P=\rho R T=4.2MPa$
    • Average pressure computed in 1ns: $P=4.5 \times MPa$
  • 4000 Helium spheres
    • $\rho=26.5686 kg/m^3$
    • Computed T=294.81 K
    • $P=\rho R T=16.281MPa$
    • Average pressure computed in 1ns: $P=MPa$
  • 100 Helium spheres
    • $\rho= 0.664216 kg/m^3$
    • Computed T=257.241 K
    • $P=\rho R T=0.35516 MPa$
    • Average pressure computed in 1ns: $P= 0.36598 MPa$

Reference

Homogeneous Cooling

With $e<1, \mu=0$, the kinetic energy and granular temperature are decaying with time. Use $e=0.9, \mu=0$

Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License