# Granular Project V2.0
###### tags: `Granular Simulations` `SoftQC` `self-assembly`
---
### Owners (the only ones with the permission to edit the main text)
Etienne, Andrea Jr., Andrea Sr., Frank, Giuseppe
---
## Background
Previous simulations showed that we can probably self-assemble binary structures in 2D vibro-fluidised granular systems.
We now fix the material parameters such that:
|$\gamma_n^{aa}$|$1.08\times10^7\,s^{-1}m^{-1}$|
|:---:|:---:|
|$\gamma_t^{aa}$|$3.26\times10^7\,s^{-1}m^{-1}$|
|$\gamma_n^{ab}$|$3.2\times10^6\,s^{-1}m^{-1}$|
|$\gamma_t^{ab}$|$2.6\times10^4\,s^{-1}m^{-1}$|
### Goals
The two main objective of this page is to:
* understand better the key parameters for self-assemby
* simulate the systems we will likely use in experiments.
### To Do (In green what has been finished)
:::success
* Effect of the roof
* Effect of real boundaries on larger ($N\sim10^4$ particles) systems
:::
* Explore other state point (QC?)
* QC DEM: Explore the possibility of longer run for better QCs. Smaller systems?
* QC DEM: Melting simulatios. Start from Random tiling and see if they melt
* Explore in MD other phases: do T1 and maybe H3 form?
* S1: Check Parameters (realistic parameters, no soften)
## Bouncing dynamics
### Roof effect : negligible
Without the roof, we observe a similar qualitative behaviour than in the simulations with the roof.
The red lines on the next plots depict the limits of the accessible region for the particles in the $z$ direction. The red shading represents the amplitude of the vibration.
For large particles:
[30/10/2020]
For small particles:
> [name=gfoffi] Just check also that the crystall forms in the same way, but it should be fine.
> [name=Etienne] Yes, simulations are running to check this
> [Side note] Large particles initial $z$ is set to $3r_L+0.1r_S$. There could be a problem if small particles were big enough to cause this value to grow beyond the maximum allowed $z$ for the large spheres. Luckily, it turns ou that this would require $q=23$, which is impossible.
Even though the roof plate seems to have a negligible effect on the dynamics, it could help LAMMPS not loosing small particles that, sometimes, get an unexpected kick and fly out of the box.
The roof is super important for smaller particles. See *Why do we loose particles* section.
## Attempt at QC composition
With the parameters that self-assembled S1 ($A=1.3e-5m, fr=200Hz, r_S = 1mm$ and $q=0.457$), I tried running some simulations, close to the composition at which we expect the QC : $N_L=1000, N_S=317$.
In thermal systems, decreasing the amount of small particles tends to lower the crystallisation pressure (ie. packing fraction). Therefore, I tried simulations at $\phi=0.855, 0.850, 0.845$ and $0.840$.
Little mouvement is observed overall, so the packing fraction might still be too high. In all cases, large particles tend to crystallise into a hexagonal crystal and expel small one to grain boundaries and small clusters.
At $\phi=0.840$:
Initial configuration
Final configuration after $5.74e9$ timesteps
[30/10/2020]
At smaller packing fractions, large particles crystallise into a hexagonal crystal and small ones remain in a fluid phase.
In all simulations, a few small particles are lost, probably due to numerical errors (see *Why do we lose particles* subsection)
## Beads size
We plan on buying metal bead of diameters 1.19mm, 2.381mm and 2.5mm. With those, we can achieve size ratios of 0.476 and 0.49979 that are both in the random tiling region.
Preliminary simulations were run with small particles of 2mm in diameter. The scaling of the vibration amplitude (possibly frequency) with particles sizes is not obvious : if one thinks with respect to the imposed displacement, a linear scaling could be considered. But one could also think in terms of force, hence mass, which scales as the size cubed.
Without further insight, it seems that our best option is to carry out a scan of those parameters, and see what works.
The simulations below are performed without a roof. When the amplitude of the vibration is too large, small particles are lost.
### $d_L=2.381mm, d_S=1.19$
Here, $q=0.49979$.
There is a slight increase in $q_4$ for $A=0.9e-5 \; or\; 1.1e-5$, but it is very hard to relate it to a S1 growth in the snapshots.
$q\sim0.5$ is the largest $q$ one can have in the random tiling region, with spheres lying on a plane. The small particles exactly fit into the vacancy at the center of the square and hence, it takes a very long time to enter or leave a square, slowing dynamics down.
For $A=1.1e-5\,m$, 2 small particles are lost in the course of the simulation. For $A=1.3e-5\,m$, all small particles are lost eventually.
### $d_L=2.5mm, d_S=1.19$
Here, $q=0.476$.
A clear increase in $q_4$ is observed for all simulations, more pronounced as the vibration amplitude increases. Chunks of S1 are clearly visible for $A=0.9e-5\,m$ and $A=1.1e-5\,m$.
2 small particles are lost for $A=1.1e-5\,m$. The corresponding boop plot stops at the time of the first particle loss, but it is clear from the snapshot of the last timestep that $q_4$ keeps on increasing.
At $A=1.3e-5\,m$, all small particles are lost.
The dynamics here is much faster, owing to the larger size difference between the two types of spheres that eases the formation of squares (small spheres can rattle in the central volume).
It seems that the "good" amplitude is around $1e-5\,m$. This confirms that the scaling of A with the size of the particles is not simply linear, and calls for a proper scan of amplitude and frequency to find a sweet spot where S1 forms quickly and reliably.
### Why are we loosing particles ?
[30/10/2020]
LAMMPS considers that a particle is lost if it exits the simulation box. Two reasons at least could lead to small particles loss:
* Numerical instability : in the input script definition, the simulation timestep scales linearly with the size of the particles, hence it is smaller for smaller particles as expected. But the real scaling might not be linear, and the timestep we use might simply be too large, leading to numerical error accumulation until a particle skyrockets.
* Physical rare events : the simulation box is 3 large diameters thick, and the particles lies on a vibrating plate located at one large diameter height. This leaves 2 large diameters in the vertical direction before a particle is lost. I could be possible that some small particles just sometimes fly high enough to leave the box.
To try distinguish between these two possibilities, I ran a new simulation in the worst case scenario $A=1.3e-5\,m$ with a shorter timestep.
The simulation with shorter timestep loses particles as well, but slower than the first one. Particles loss rate increases as time goes, which could be due to numerical errors accumulation.
In addition, the standard deviation of $v_z$ for the small particles increases as simulation goes, which is consistent with momentum accumumlation due to numerical errors.
I also ran a simulation with the roof plate added back at a height of $6/5 d_L$ above the bottom plate, and with the fastest timestep. Then, no particles are lost and a S1 crystal forms.
With the roof, there is still an increase in the standard deviation of $v_z$ for the small particles, but shallower. Hence, the roof helps keeping the particles inside the simulation box.
The results above suggest that small particles accumulate vertical momentum due to numerical errors and gradually escape the box. The roof slows down this accumulation to some extent.
To confirm that numerical instability is actually at play, I started a simulation at the same state point with an even smaller timestep (dtfactor = 0.2).
[02/11/2020]
The simulations with dtfactor = 0.2 also looses all of its particles, at a slower rate than when dtfactor = 1.0, but faster than dtfactor = 0.5.
Is that really numerical instability ? I doubting now.
[04/11/2020]
I ran a simulation resolving short time dynamics around the time of first particle loss.
In the usual regime, the small particles are kind of trapped under the layer of large ones. It appears that the particle that is about to get lost escapes from its hole and starts bouncing on top of the large particles layer.
From there, it is much closer to the box boundary, and a sufficient kick will push it out of the system. Also, since the large particles are also bouncing, they act as a driving plate with even larger oscillation amplitude than the actual driving plate. This reminds a bit of the *Galilean cannon*.
{%youtube l6EFJJKmvlc%}
I suspect (but did not check) that the same phenomenon occurs for all lost particles.
The trajectory of the particle looks sensible.
Moreover, the observed drift in the standard deviation of $v_z$ could be due to the gradual change in packing fraction as particles escape. The lower the packing fraction, the less collisions to transfer vertical momentum to the horizontal plane and hence the more energy in the vertical direction. Similarly, the small drift observed even with the roof is present could be due to the self-assembly of the S1 crystal leading to a more efficient packing, leaving more space for each particle to rattle.
I now believe that we are not facing numerical instability, and that the roof plate prevents small particles to get on top of the large particles layer.
## Refixed parameters
With the corrected formulas for kts, at the point where Andrea did the first promising simulations ($q=0.855, x_S=0.5, A=1.3e-5\,m, fr=200\,Hz$), S1 still forms so the small change in the parameters indeed does not change the physics dramatically.
Crystallisation speed is comparable to that observed with the other parameters at the same point.
|Roof | No roof|
|:---:|:---:|
|||
|||
## Beads size, second take
[05/01/2021]
A more thorough scan of the parameter space is conducted, with the refixed parameters and keeping the roof plate in place. Beads diameters are set to 1.19 mm and 2.5 mm, corresponding to the real ons. The frequency and amplitude of the vibrating plate are varied.
Simulations last $10^{10}$ timesteps, corresponding to a real time of ~55.8h.
For a rougher overview, the next figure displays the value of $q_4$ at the last time of the simulation as a colormap.
As ususal, for small amplitudes and frequencies, nothing moves.
Then, for larger amplitudes and frequencies, S1 self-assembly takes place. Finally, when to much energy is pumped into the system, small particles pile up and a hexagonal lattice of large particles forms.
[12/01/2021]
In an attempt to reduce the number of parameters, I computed the reduced acceleration $\Gamma = (2 \pi f)^2 A / g$ for each simulation and plotted the last value of $q_4$ against it. The colors of the dots correspond to the frequency from 100Hz (purple) to 400Hz (yellow).
A value of $\Gamma \sim 4$ seems preferable. However, the reduced accelaration does not quite accurately collapse the simulation results : some simulations at $\Gamma \sim 4$ have a poor $q_4$.
It turns out that the product $A \times f$ (average speed of the driving plate $v_{plate}$) yields a better collapse.
This quantity is traditionnaly considered as the *shaker's strength* defined as $(2 \pi f A)^2 / (x_S d_S + (1 - x_S) d_L) g$.
S1 self assembly takes place around ${v_{plate}} = 0.003$m/s.
At low frequency and amplitude, particles dont move except vertically, following the plate.
For the smallest frequency (100Hz) and amplitudes (7 and 8 µm) though, same movement occurs.
For an amplitude of 8 µm, a strange behaviour is observed in the early stages of the simulation.
## Strange early dynamics
[11/01/2021]
In the early stages of the simulation at 100 Hz and 8 µm amplitude, some voids appear around specific particles. These voids move around the system, merge together sometime, and eventually vanish. Then, the system follows a more sensible dynamics.
{%youtube L6s5QHqaZfE%}
The simulation was restarted, and the same phenomenon occured. It also occured for a system of 2000 particles of each type (hence with different initial conditions).
It appears that the voids contain particles bouncing offsync with the rest of the system. The next figure displays a snapshot of the simulation with particles colored according to their z-velocity. This is observed on all frames. Moreover, off-sync particles tend to have a velocity component in the z direction ~10 times larger than the others in absolute value.
These energetic particles certainly create the voids by pushing their neighbours away.
Why they end up bouncing off-sync is a mystery to me though. The initial configuration is a desordered dense packing obtained by a rapid compression in MC simulation. All particles are initialised at the same z-position. There is no reason why some would dramatically differently than others. The voids seem to appear homogeneously in the early steps.
Only the simulation with these parameters (100 Hz, 8 µm) displays this strange dynamics. At 200 Hz or 100 Hz at $A \leq 9$ µm, nothing moves. At 100 Hz and 7 µm, surprisingly, particles move, but no voids appear.
[15/01/2021]
**P.B. Umbanhowar, F. Melo and H. L. Swinney, *Localized excitations in a vertically vibrated granular layer*, Nature Letters, 1996, 10.1038/382793a0**
Authors observe stable circular excitations oscillating between crater and peak shapes (they call them *oscillons*) in vertically vibrated granular systems. They use bronze beads of 0.15-0.18 mm in diameter. Oscillons are observed for driving frequencies between 20 and 30 Hz and shaking accelerations $\Gamma$ between 2.5 and 2.5. Oscillons require a particle layer of at least 13 particles diameter depth.
Oscillons do not appear spontaneously. They are obtained from patterned states at higher $\Gamma$ and nucleate at the defect of the patterns when $\Gamma$ is reduced (hysteresis).
The oscillons interact. Two oscillons of opposed phase attract while those in phase repell. Dipole, polymer chains, triangle clusters and square lattices are observed. Adding harmonics to the driving vibration allows to stabilise a hexagonal "crystal" of oscillons.
Scaling of critical frequency with particles size $f \propto D^{-1/2}$ suggest that oscillons form in a region of constant (high) dissipation.
This does not seem to be super close to our system : oscillons require a 3D stacking of spheres to form, which we don't have.
Moreover, We observe voids formation at 100 Hz and our particles are much larger than theirs, so we would expect to find oscillons at frequencies smaller than 20 Hz according to the proposed scaling.
Very nice to see self assembly of such excitations though !
**E. Clément, *Granular packing under vibration*, EPS, 1998, 10.1007/978-94-017-2653-5_42**
Overview of unique features of vertically vibrated granular systems, including oscillons of Umbanhowar et al.
No proper references unfortunately.
## Hard walls
[13/01/2021]
I modified the LAMMPS input script to add hard walls on the sides.
Initial configurations need to be taken from a proper MC compression with hard walls, otherwise the space required around the particles chunk to avoid initial overlaps with the walls is too large to reach packing fraction at which crystallisation occurs.
Some preliminary simulations show that the setup is working though.
Therefore, I modified my MC program to support hard walls. Simulations using these initial configurations are now running.
[03/02/2021]
Granular simulations with 2000 particles of diameter 1.19mm, 2000 particles with diameter 2.5mm.
Simulations are performed with either fixed or wiggling hard walls. I am not 100% sure that the wiggling hard walls are in phase with the plate, but this would make sense.
At first sight, the impact of the dynamics of the walls on the self assembly is unobservable.
* Frank's intuition: squares are favoured by the wiggling of the central small disk. At the wall, there lack one small particle on one side, so the entropy gain by forming a square might be smaller than that obtained from forming other structures (eg. hexagonal packing).
Surprisingly, in all simulations, S1 self assembly takes place away from the walls. Self assembly seems more efficient in simulations with periodic boundary conditions although comparison is not straightforward since periodic systems contain twice as less particles.
It seems that a square box does not play the template role that one could have expected. Instead, a layer of fluid interfaces with the walls.
>[name=frank] twice as less = half as many? More particles is likely to lead to more domains, generally speaking. The walls also probably help with more domains (as they cannot cross the box boundaries now). So that might be the main factor in lower overall order.
Finally, no simulation seems to have reached equilibrium in $10^{10}$ timesteps.
### 300 Hz, 10 µm
* Wiggling hard walls, $\phi = 0.855$
{%youtube mrSJyhcSyAU%}
* Periodic boundary conditions, same system (less particles)
{%youtube kYoDoHRWSc4%}
### 350 Hz, 9 µm
* Wiggling hard walls, $\phi = 0.855$
{%youtube n9cC6qRUZcU%}
* Periodic boundary conditions, same system (less particles)
{%youtube 3bAqqw-U2go%}
### Longer simulations
[08/03/2021]
* Wiggling hard walls, $\phi = 0.850$, $fr = 300$Hz, $A = 10$µm, $2\times 10^{10}$ timesteps
* Wiggling hard walls, $\phi = 0.850$, $fr = 350$Hz, $A = 9$µm, $2\times 10^{10}$ timesteps
The S1 crystal at $fr = 350$Hz clearly nucleates a grow away from the walls. The other aligns with the wall at the bottom.
The presence of walls does not clearly promote S1 self assembly, but does not prevent it completely either. The proof of concept of self-assembly in granular media seems robust.
## Large systems
[15/12/2020]
Scaling analysis for systems of 10k particles of each type.
16 cores seems like a reasonable choice for these simulations.
[08/03/2021]
* $r_{LL} = 1.25$mm, $r_{SS} = 0.595$mm, $fr = 300$Hz, $A = 10$µm, $N_L = N_S = 10000$, $\phi = 0.850$, wiggling hard walls
{%youtube 3zgppKgsQEQ%}
* $r_{LL} = 1.25$mm, $r_{SS} = 0.595$mm, $fr = 350$Hz, $A = 9$µm, $N_L = N_S = 10000$, $\phi = 0.850$, wiggling hard walls
{%youtube 6WBO0MuTHxQ%}
In both cases, S1 nucleates all over the box, and the nuclei grow quickly until they start feeling one another (or the walls). Then those crystals need to rotate to eliminate grain boundaries, which is a much slower process (not completed in the simulated time $2.007 \times 10^{5}$s).
Each simulation took about 14 days to complete on 16 cores.
[01/04/2021]
Some results about large simulations ran at different composition.
S1 chunks form down to $x_S = 7/20$, but slower and slower. Two effects play in this slow down. First, with less small particles in the system, it is harder to form squares. Second, the crystallisation density is expected to go down as the number of small particles is decreased, and all simulations are run at the same packing fraction of $\phi = 0.850$.
Below $x_S = 7/20$, the system is actually jammed and small particles simply rattle in their cages.
I restarted some simulations with a lower packing fraction at these compositions.
* $r_{LL} = 1.25$mm, $r_{SS} = 0.595$mm, $fr = 350$Hz, $A = 9$µm, $N_L = 10000$, $N_S = 9000$, $\phi = 0.850$, wiggling hard walls, after $1.58 \times 10^{10}$ timesteps.
* $r_{LL} = 1.25$mm, $r_{SS} = 0.595$mm, $fr = 350$Hz, $A = 9$µm, $N_L = 10000$, $N_S = 8000$, $\phi = 0.850$, wiggling hard walls, after $1.64 \times 10^{10}$ timesteps.
* $r_{LL} = 1.25$mm, $r_{SS} = 0.595$mm, $fr = 350$Hz, $A = 9$µm, $N_L = 10000$, $N_S = 7000$, $\phi = 0.850$, wiggling hard walls, after $1.74 \times 10^{10}$ timesteps.
[21/05/2021]
Some news of the simulations at lower density.
* $r_{LL} = 1.25$mm, $r_{SS} = 0.595$mm, $fr = 350$Hz, $A = 9$µm, $N_L = 10000$, $N_S = 7000$, $\phi = 0.840$, wiggling hard walls, after $2.41 \times 10^{10}$ timesteps.
The configuration is not jammed, there are fast a slow moving regions. The diffraction pattern of the last frame displays 12 peaks, quite blurry (plus buggy pixels from my program).
The 12 peaks are located at a $|k| \sim 5000$ which correspond to a direct space distance of $1.26\;mm$. No order is visible at smaller $k$ so the direction of the bonds is ordering in 12 directions, but the resulting random tiling does not have a long-range 12-fold symmetry.
The next video displays the bond network between the large particles. There is a squeletton that barely moves, and highly mobile region inbetween, with a high concentration of small disks. At the boundary local realisation of what looks like a RTQC are sampled. (I can share a higher quality video if necessary, youtube compresses to a crazy degree)
{%youtube AliK8dQszvM %}
The same behaviour is observed at $x_S = 6/16, \; \phi = 0.835$.
At lower concentrations of small disks the system demixes.
| $N_L=10000\; N_L=5000\;\phi=0.835$ | $N_L=10000\; N_L=4000\;\phi=0.825$ | $N_L=10000\; N_L=3170\;\phi=0.825$ |
| -------- | -------- | -------- |
|  | | |
I have made an incredibly stupid mistake. The QC *composition* is $x_S \approx 0.317$ and in the simulation, I looked for it at a *number ratio* $0.317$ of small disks... This results in a much smaller $x_S$, which of course favours hexagonal packing.
For $10^4$ large disks, we need $4641$ small disk to reach the taget concentration. I will start simulations at this point.
The simulations with $5000$ small disks that we already have demixe at packing fractions of $0.835$ and $0.840$ though.
[06/07/2021]
The simulation with $N_L=10000$, $N_S = 7000$ ran for $4\times 10^{10} \tau$ more. The diffraction pattern and direct space snapshots are essentially unchanged.
If some annealing at play, it is remarkably slow.
I started some smaller simulations in periodic boxes with $N_L = 5000$ and $N_S = 2321$ (RTQC composition). For packing fractions between $0.830$ and $0.844$ (with step-size $0.002$), the large particles demix and form a hexagonal crystal in mobile regions of large and small disks.
Example at $\phi = 0.840$ after $2.7\times 10^{10} \tau$
As was suspected earlier, it seems that an excess of small disks is necessary to assemble the RTQC.
## Comparison with thermal
For the granular experiments, we have beads of diameters 1.19mm, 2.381mm and 2.5mm. This gives us access to size ratios $q = 1.19/2.5 = 0.476$ and $q = 1.19 / 2.381 = 0.49979$.
The non-additivity parameters that correspond to spheres on a plane for these two size ratios are respectively $\Delta \approx -0.06514$ and $\Delta \approx -0.05725$.
From the infinite pressure phase diagrams at these non-additivity values, we can expect T1, and maybe H3 to be stable as well.
|||
|:-:|:-:|
|$\Delta = - 0.06514$|$\Delta = - 0.05725$|
## Granular Quasi-Crystal: QC8 edition
Lammps script to generate initial configurations with the desired: size ratio, composition, packing fraction, total number.
## Perspectives
* **Avoiding piling up**. Regarding the scan of [05/01/2021] and the plots of [12/01/2021], for low values of A and f nothing moves while for high values of them the S1 is destroyed because small particles pile up. Crystals like H2 with q>0.5 should be immune to this problem because one can design the height of the roof in such a way to be slightly higher than the large diameter but less high than the sum of two small diameters. I then suggest to check if also a crystal of this kind can be reached dynamically in MD so that we can also try to obtain it in the granular one. Pssible choices: Delta = 0.03 q = 0.6 -> H2 and H1, Delta = 0.05 q = 0.52 -> H2 and T1.
* **Changing structure with A/f**. Are there values of q (possible > 0.5) and xs with \Delta=OkForGranular where two different crystals can coexist at infinite pressure and also in thermal MD? Because an interesting thing to do in my opinion would be perform granular simulations in that geometrical conditions and try to see if varying the driving parameters one can switch from one structure to another. This would be interesting also for what Emmanuel Trizac suggested i.e. try to observe phenomena in the granular systems that are not possible in the thermal one. This would be the case: for elastic hard spheres there are just the geometrical constraints while in the granular one we have also the driving. In the conditions for which the geometry predicts the coexistence between two structures can the driving parameters select one of them as the more stable? Is it possible to observe dynamically (i.e. varying A,f to A',f' during a simulation) such a change of structure?
This would be also a field of application for the optimization protocols. A question like "which is the way to vary A and f that minimize the time to switch from a structure to another?" is suitable for them. Idea: The numerical part of this point can be done by the student Raphel.
* **Look for quasicrystals in granular simualtions** following the ideas contained in https://arxiv.org/abs/2202.12726 (pre print of Eitenne et al)
Sign in with Wallet
Connect another wallet