# 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 100^{3} 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**

- Relevant relations: http://en.wikipedia.org/wiki/Kinetic_theory_of_gases
- Noble Gas Properties: http://en.wikipedia.org/wiki/Noble_gas

## Homogeneous Cooling

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