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.

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
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$

page revision: 36, last edited: 29 Mar 2009 21:54