# Granular pressure (Matteo)

###### tags: `M2`
- LAMMPS based simulations (including friction) of spherical beads on a vibrating plate

## To Do

To help with the interpretqtion of the experiments, here are a few plots of the Fourier transform of the force felt by the shaker by a single particle for different height values (here r=2.5 mm, amp = 0.03) h=1.1 (partially synchronized):
h=1.2 (partially synchronized):
h=1.3 (synchronized) :
h=1.4 (synchronized) :
h=1.6 (fully asynchronous) : 

So, I cheked the time dependence of the number of active particles in the case of a synchronous transition (r=2.5 mm , h = 1.55, amp = 0.03), and there seems to be a meta-stable state, phi=0.175:
phi=0.15 :

However for the continuous case, the curve seems to be as an exponential near the critical value (r=2.5 mm , h = 1.2, amp = 0.03), without any meta-stable state, phi=0.3 :
phi=0.325 (for this case, maybe I should check by making the simulation longer):
phi=0.35 :

> [name=gfoffi] Can you put the definition of the synchronization parameter you are using?

So, here the synchronization parameter is given by : $s=\frac{v_{\text{cm}}}{\frac{1}{N}\sum{v_i}}$, where $v_{\text{cm}}$ is the variance in position of the center of mass, and the $v_i$ are the individual variances of the position of the different particles.
So in a synchronous region, where essentially all the $v_i$ are equal, we must have $s=1$, and if they are completely asynchronous, we expect $s\approx\frac{1}{\sqrt{N}}$.
Then there are the intermediate regions, where the degree of synchronicity seems to take intermediate values, for which we have seen some plots during the last brainstorming session, and which indicates intermediate synchronicity.

Upgrade : new and more accurate figure for the case r=2.5 mm : 

And for r=1.25 mm :

It seems that overall, the effect of decreasing the radius, effectively increases the size of all the regions, allowing us to see a bit better the valleys inside the synchronous regions, we saw before.

So here are the first plane figures for a radius of 2.5 mm, with the degree of synchronization defined during our last talk, the data isn't great because I didn't do very long simulations, but it still gives a first idea of what to expect (I tried some interpolation on the function imshow, to smooth the edges of the regions, but the result is never really great) :
And if anyone fancies 3D diagram :

We can clearly see two well-defined regions of synchronicity and a synchronicity, with at least one well-defined valley in the synchronicity region, with maybe another one (or other ones), though these valleys are thinner when we increase the amplitude.

There is also some funny behavior for the lowest amplitude, but at this amplitude, the particles probably barely "jumps". Indeed when I look at the ovito simulations, only a part of the particles continues to jump during the whole simulation, it probably depends on the initial z coordinate, thus for the lowest amplitude, the behavior comes mostly from the random variables I use to create the initial conditions. Though this behavior also depends on the height.

The area explored before, was for amp = $0.03\sigma$, for which we see that the aforementionned valley is at its widest and deepest, and the synchronocity region is at its smallest, which explains the very diverse phase diagrams we observed before.

Here are a few spectrum of the vertical motion of the particles for different height values, first h=1.1 :
For h=1.2 :
h=1.3 :
h=1.4 :
h=1.5 :
h=1.6 :
h=1.7 :
h=1.8 :

Notice that we can always see a pic at $\omega=53 \text{Hz}$, and some of its multiples. What seeems to be the defining feature of the synchronous behavior, is that there seems to be less noise.

Here is a more complete phase diagram : 

Phase diagram for only synchronous parameters : 

And for non-synchronous parameters : 

(For the last figures, amp = 0.00014827 m, which is about $0.03 \sigma$)

For the phase diagram at fixed amplitude, there is a slight problem in that the transition between the curves is not very nice : 
(Notice the very late discontinuous transition at height 1.3 the green points)

We can notice on the ovito animations, that in the height where there is a discontinuity, the particles tended to synchronize their vertical movements, around the same time the system actually dies.

*(AP: this vertical synchronization is not an artifact! I found this paper where they mention it http://dx.doi.org/10.1103/PhysRevLett.106.088001 they also have videos and it is exactly the same. My guess is that you can have the transition both with sync and non sync. In the first case it is discontinuous in the second case it is continuous. We can talk about that next week.)*

After trying with a randomized initial velocity along the vertical axis (as before it was not the case, which might have accounted for the weird behaviors), I still find the same kind of plots : 
Also, on the ovito simulations, the particles synchronize their vertical movements, just as before.

Here is a phase diagram for h=1.8, but different values of amplitudes (which are, this time around, expressed in terms of diameters of the balls) : 
We see that the transition seems to get steeper as we increase the amplitude, I guess it is possible that there is a transition from 2nd to 1st transition, as Raphael suggested.
The curve for the smallest amplitude is a bit weird, but I think that it is due to the fact that for this value, the particles are free.

- Particles in a shaker, aligned -> we see how far they go depending on their initial velocity

First plot, for amplitude $A=218.05 \mu\text{m}$ : 

The height here is expressed as the factor of balls diameter.
For amp = 0.00014827 m : 

And amp = 0.00030527 m : 
This curve seems a bit off, I'll have to check if I didn't make a mistake somewhere.

Here is for different amplitude as given by amp = [0.00014827 0.00021805 0.00030527]m and height =1.2 : 
I am not sure about this last figure, the next two seem better.

For h=1.8 : 
For h=2.4 : 

Upon further study of these small simulations, I have found some trends which might me be interesting. I plotted the velocity of the quickest ball as a function of time, if I put the y-axis in a logarithmic scale, I obtain, for h=1.2 : 
They look like rather clean lines, suggesting an exponential decay of the velocity as a function of time (also I switched the labels of the x- and y-axis again...)

If I do the same for h=1.8, the lines are less clean : 
Actually, for the case of low amplitude, the ball doesn't always hit the ceiling and only lose velocity stepwise.

Then for h=2.4, the lines become really disgusting, for the same reason as above : 

If now, for each of these cases, I fit an exponential to the velocity as a function of time for each particle and plot this decay rate as function of the initial velocity, I obtain for h=1.2 : 
For h=1.8 : 
For h=2.4 : 

Thus, we see that the decay rates seem almost constant for each of these cases. For now, I have not seen any clear relations between these different decay rates, maybe I should do more experiments, if we want to look further into it. What is however noteworthy is that the dependence of the decay rate with the height seems to flip between h=1.2 and h=1.8.

For completeness, here are the decay rates for amp = 0.00014827 m : 

For amp = 0.00021805 m : 

For amp = 0.00030527 m : 

The last figure is the most weird, in that the decay rate doesn't depend monotonically with the height.

- Full simulation for experimental values (400 particles)

Bad news, in my simulation, for $\phi=0.05$, the particles are still active : 

And here is the phase diagram : 
With the timescale to stationary state : 

- Full simulation for a lot of other values (4900 particles)

First results of the experiments for different values of height but same amplitude = 0.00014827 m : 

Then for amp = 0.00021805 : 

For amp = 0.00030527 : 

If this time, we fix the height at h=1.2, but vary the amplitude, we have : 

And h=1.8 : 

And h=2.4 : 

From these few diagrams, we can see that the small range of parameters have a huge impact on the phase transition (this is really not an understatement), just to be sure I checked some extreme cases, just to be sure that I didn't mess up the plots. Overall, from the animations, my feeling is that the transition seems to really depend a lot on the how much a particle can travel without hitting another one.

Some radial Structure factor plots for different packing fractions, I don't really trust them, but in case it seems good, I am still putting them (I am the most suspicious about the $\phi=0.275$). The main result is that there is absolutely nothing special at the critical value.
For $\phi=0.1$ : 

For $\phi=0.25$ : 

For $\phi=0.275$ : 

For $\phi=0.3$ : 

For $\phi=0.35$ : 

For $\phi=0.45$ : 

## Pressure

If we want more details on the behavior of the pressure, we can look at the density and the temperature near the wall. For that, I compute these densities for particles near the wall, i.e. the ones whose distance to the walls are less than $1.5r_{high}$. We can then plot these quantities as a function of the packing fraction. Here is the plot for the density : 
The density near the wall seems to grow linearly with the packing fraction.

Then we can plot, the temperature near the walls : 
As seen before, the temperature grows more quickly, which we assumed to be the main factor behind the rapid growth of the pressure as a function of the packing fraction.

In order to have an idea of how the computed quantities relate to the pressure, we can plot along with the pressure, the product of the temperature and the density near the walls. This produce this figure :

The agreement seems between the two curves seems to agree with our naive hope to make a link between the measured quantities, as one would expect from an ideal gas. However, if we directly plot the ratios, we can see that the agreement is not very good, and shows us that the system is quite far from an ideal gas, as one would expect in this case a constant function : 

To complete our study of pressure, we look at the study of pressure in two very similar systems : in the first, the walls wiggle with the shaker, and in the second, we remove the roof of the shaker.

If we superpose the pressure measured in the case of wiggling walls and fixed walls : 

As we can see, there does not seem to be any difference between these two cases.

If we look at the plot of the density, of the temperature and of the ratio of their product with the pressure, it reinforces the idea that the wiggling of the walls doesn't affect the system a lot : 

However, if we compare the scenario with and without a roof (both situations with fixed walls), there seems to be quite a difference : 

Note : I didn't put the situation where the packing fraction is equal to 0.8, because I quickly lost a lot of balls (I should probably change the initial conditions to prevent that).

If we first look at the density near the wall, there doesn't seem to be much differences between the two scenarios, except maybe at low packing fractions where the particles are much more eager to stay near the walls, and the linear trend of the desnity near the wall is much more obvious :

However, the biggest differences are seen when we compare the temperatures near the walls : 

When there is no roof, the temperature doesn't increase as much, as it is already quite important at the beginning.

As a bonus figure, here is the comparison of the ratios $\frac{E_{\text{kin}}n}{p}$ between the two aforementionned cases :

## Active Particles Transition

As requested, here is the plot of the proportion of active particles between the cases $N=4900$ and $N=400$ : 

The critical point seems to not have changed much as previously stated. The non-zero behavior of some points below the transition point stems probably from the fact that the simulations for $N=400$ were done on a much smaller time interval.

## Meeting 10/05/22

- For Pressure :

Try and do same experiments, if the walls are wiggling with the shaker and see if the pressure and kinetic energy change because of that.

Plot the density and kinetic energy near the wall as a function of the packing fraction -> look at these quantities for approximately a radius away from the wall.

- For active particles :

Plot the proportion of active particles for the two sets of experiment (different number of particles) on top of each other, using same velocity threshold just to better see how similar the two sets of experiments really are.

Plot radial structure factors -> to have better curves do the experiments multiple times

## Pressure (2nd and 3rd week)

We perform different experiments, with different packing fraction, and different wall viscosity (we change just the normal component). We always keep the same random force on the particles with $\sigma=0.001$. We then measure the pressure on the wall and we have : 

It seems that the wall viscosity doesn't have much an influence on the overall pressure, except maybe at low packing fraction (I should plot the errorbar, but I am having troubles computing the standard deviation for now...). We can however note that at high packing fraction, the points converge even more, probably from the fact that there are more collisions, making the energy loss from the wall collisions negligible.

For a broader range of values, we obtain : 
The importance of the viscosity is now more apparent, though again for low packing fraction and only for the most important viscosity.

To have a better idea of the effect of viscosity, I can plot the coefficient of restitution for a collision normal to the wall : 

and for a broader range of parameters : 

As expected the highest coefficient of restitution appear for the lowest value of viscosity.

If we measure the pressure for different wall stiffness, we have : 

Then for the 1 particle case : 

If I take away the gravity and the plates, I instead obtain this figure :

Which exhibits the complete opposite behavior as the previous figure !

- Density and kinetic energy near the wall.

As a follow-up on the discussion, I plotted the denisity of particles as a function of their distance to the nearest wall.

Here are the figures, first for $\phi=0.05$ : 

for $\phi=0.1$ : 

for $\phi=0.2$ : 

for $\phi=0.4$ : 

for $\phi=0.6$ : 

for $\phi=0.8$ : 

We observe some discrepancies when the radius is large, but I think it is mainly linked to the fact that for large radii, the uncertainty is bigger ( as we are then in the middle of the box, and the nearest all can be any of the 4 ones). This is not so important as we are mainly interested at the curve for small distances of the wall.

The first observation we can make is that the density near the wall increases when we increase the paclking fraction, as expected. However the increase is definitely not exponential, and therefore doesn't, explain the exponential increase of the pressure.

A second observation we can make - though it hass no impact on the pressure, I think - is the pics appearing at high packing fraction, which signals us the ordered structures appearing then (which is also a goodway for me to know, that the curves I display seem alright).

Then, I can do the same for translational kinetic energy in the plane.

We obtain the figures for $\phi=0.05$ : 

for $\phi=0.1$ : 

for $\phi=0.2$ : 

for $\phi=0.4$ : 

for $\phi=0.6$ : 

for $\phi=0.8$ : 

In these figures, lies probably our best cahnce to explain the behavior of the pressure, as we notice that the near the wall, the temperature increases drastically, as we increase the packing fraction.

- Updated plot for the pressure

After making 5 different simulations for each packing fraction, I have the following figure for the pression :

It seems to confirm the trend exponential trend seen before. Also, we can notice that at low packing fraction, the standard deviation is higher, which may provide an explanation to the fact that in the figures above, the pressure varied more in the low packing fraction cases than in the higher ones, even though the figure above, the conditions of the experiments varied for each measure of the pressure.

- Moving wall with two different sides. Start with the same density initially on both sides.

I made these experiments at one point, and kind of forgot about them...

We make the same kind of experiments as above, for the measure of the pressure at the walls, except that we put a wall in the middle of the box, and measure the difference of forces there. By first putting the same density on either side of the wall, we obtain (labeling issue here the x-axis is where the middle wall is situated, 0.5 being the middle of the box): 
Again, I need to compute the actual errorbar, but from what we see, it seems that the net force difference at the wall is null.

If however, we put the same numbers of particles for both regions of the box, we obtain (again I need to put sone error bars and the x-labeling is wrong): 

In this case, we note that that for each value of packing fraction, the force difference is indeed dependent on the placement of the wall. We also note that the higher packing fraction, the higher the magnitude of the force difference, as one would expect.

## Extract the fraction of *active* particles (2nd and 3rd week)

For the computation of the "active" particles, we are presented the issue of defining what is an actual "active" particle in our case.
After examining multiple histograms of the norm of the velocity, there doesn't seem to be a region of the packing fraction where the velocity norms would split. Thus as a hand-wavy argument, I just take a threshold value for the velocity, above which the particle is said to be active. I then plot the proportion of active particles as a function of the packing fraction. In our case, I found the threshold value of $v=0.005 m.s^{-1}$ seemed to give a rather good distinction between "dying" systems and systems that are more "alive", thus I obtained : 
The figure looks like a what we would expect, with a nice phase transition. If we change the threshold value, we obtain the same kind of figure, except that the proportion of active particles will grow differently. We also observe that the "critical" packing fraction will change according to the threshold.

> [name=gfoffi] Nice! I think the phenomenon is there and I am sure that we can came with a better definition for the fraction of mobile particles. I am really wondering what is happening for the case in which we have a small fraction of active particles, i.e, $10\%$. Do we have a population of active and quiescent one? Do they exchange status?
It would be very interesting to see a video of the steady state as well as the way to the absorbing state. I think you can do this easily with ovito. Do they exchange status?
It would be very interesting to see a video of the steady state as well as the way to the absorbing state. I think you can do this easily with ovito. It would be also interesting to see one of the histograms at the same point. Perhaps we need to use a criterion based on the MSD?

Here is the histogram for packing fraction $\phi=0.325$, for which we have about $10\%$ of active particles : 
It seems that the particles tend to exchange status, which is why, the PDF is rather smooth, as ex-active particles tend to decrease their energy smoothly.

> [name=gfoffi]
Interesting... this would mean that we have a very different behaviour from Pine and co. where a particle when dead is dead forever! So it would be interesting to test this hp. Perhaps you could plot the kinetic energies of all the particles as function of time. This might be a messy plot, but will show if particles are dying or reactivating. Alternatively, you could plot an histogram of the stdv per particle. Perhaps we can see something ? Bottom line, it would be nice to confirm this observation.

Here is a plot of the kinetic energy of 10 particles of 400 as a function of time : 

We can see that there is not one particle that stays alive forever, ordead for that matter.

> [name=gfoffi]Now I am curious, can you make a video of this trajectory with Ovito?

> [name=MatteoS] I have one, it's just I don't know how to put it on HackMD
> First try of a video $\phi=0.325$ : 

Here is a plot of the proportion of active particles as a function of time for several values of packing fraction along the "critical" intervalle : 

Note : I chose the threshold value as there was often a small bump in the curve slightly above this value which is not present in this histogram, but it is for $\phi=0.3$ : 
or $\phi=0.35$ : 
or $\phi=0.45$ : 
or $\phi=0.5$ : 

But I don't know if these bumps are meaningful or just a coincidence. Perhaps a criterion based on the MSD, would be more suited here.

Reproducing the experiments with more particles (4900), using the same method as above, we obtain a similar figure : 
In particular, the critical packing fraction doesn't, seem to have changed a lot. It might have decreased a little bit, as now for $\phi=0.275$, there are some active particles, though the time to reach a stationary state is quite long.

Indeed, as expected the timescales near the critical point tend to infinity. If we fit an exponential decrease to the proportion of active particles as a function of time, we obtain this figure of timescales per packing fraction : 

For now, I haven't really succedded at using MSD to determine the number of active particles.

One possibility would be to try and plot the MSD of every particles on a certain interval of time in the stationary state. This approach shows a nice figure where we observe indeed a phase transition, and as the packing fraction increases the activity of the system decreases due to the inherent structure, but it doesn't really give us any proportion of active particles : 

If I just put some threshold on the distance travelled by the particles, to determine the particles, I haven't obtained any "nice" figure, in the sense that I generally obtain a step function or a constant function. I don't put any figure now, because maybe I just have to tweak the value of the time interval or the threshold value to have a better figure.

## First Week

For now I have conducted three types of numerical experiment (using LAMMPS). The first ones were with one size of steel particle confined in a plexiglass shaker (also with plexiglass walls). The second ones were done without the top cover plate of the shaker, and finally the last ones were done with two types of particles in the shaker.

For each of the series of experiments, I gave to each of the particles a random initial velocity in the x-y plane. Furthermore, each of them had a random force in the x-y plane applied to them (which could in theory emulate roughness). This random force was given as a random Gaussian variable which changes discontinuously for each period of the shaker and whose standard deviation is a one of the two variables of the experiment $\sigma$, the second variable being the packing fraction $\phi$.

For the two sets of experiment, I simulated 400 particles, then for the last ones, I simulated 400 big particles and 400 small ones.

Then for each of the types of the experiments, we wait until the simulated state reach a stationary state, and we can then measure meaningful values.

The first of these meaningful values is the translational kinetic energy, that we plot as a function of the packing fraction.

For the first experiments, we obtain :  
The first figure shows only the total translational kinetic energy in the x-y plane, and the second one includes the vertical contribution.

We notice that in the first figure, for low packing fraction and for low $\sigma$, the kinetic energy is almost null, however as the packing fraction becomes more important, the value of the kinetic energy rises and reaches values of more important $\sigma$. This seems to suggest that for such packing fraction the collisions are more important than the random kicks.

The second figure doesn't tell us much more, we can only notice that the largest part of the translational kinetic energy comes from the veritcal velocities, and that for high packing fraction, the total kinetic energy decreases as one would expect from granular materials (in this case, the collisions which are more numerous dissipate more energy).

For the second series of experiment, we obtain a similar figure (we show only the energy from the x-y plane velocity) : 

However, that the main difference is that for low kicks, the kinetic energy is not zero. This can be probably explained that without the top cover plate, a lot less energy is being dissipated and thus the system can retain some energy in the stationary state.

Finally, for the case of two different particles, we may look at different plots, according to the energy of which particles we are looking at. However as we are mainly interested in the kinetic energy per the packing fraction, I am only putting, for now, the total kinetic energy of both the particles, as whether we are looking at one particle or the other, the overall behavior doesn't change much : 

We can notice that, like our previous experiment, for low packing fraction and for low kicks, the kinetic energy in the x-y plane is not zero. We can try and explain this, by saying that small particles essentially doesn't feel the top cover plate, thus they are in a situation similar as above, allowing our system to maintain some kinetic energy in the x-y plane. Otherwise the behavior of the rest of the figure is pretty much identical to the ones above.

We can also plot the ratios of the kinetic energy for the different situation(here we plot the quotient of the kinetic energy of the bigger particles over the smaller ones) : 

- Perform simulations of a single particle

I have performed a simulation of just one particle with or without boundary conditions. The particles quickly gains a rotation which makes it cancel its tangential velocity at each collision with the shaker, so that it doesn't lose any energy, so its rotational velocity is given by : $v=\omega R$, with $R$ the radius of the particle. Indeed for the periodic boundary condition case, we find at the end : $v_x\simeq-5.2 \text{ mm.s}^{-1}, v_y \simeq 5.8 \text{ mm.s}^{-1}$ and $\omega_x\simeq-2.8 \text{ rad.s}^{-1}, \omega_y\simeq-2.5 \text{ rad.s}^{-1}$, thus, as $R\simeq2.08\text{ mm}$, we have both $v=\omega R$ and $\overrightarrow{v}\bot\overrightarrow{\omega}$ (actually the relation are verified at all the decimals given by the simulation).

> [name=gfoffi] not clear perhaspt you meant $v=\omega R$ for the tangetial velocity? In principle we should be able to show, that the tangential velocity is equal to the translation velocity as suggested by Andrea. It wouldbe really interesting to have the same graph for the case with ceiling just as a curiosity.

To makes this more apparent here are the graphs of the kinetic energy in the plane for the experiment without a top plate and with periodic boundary conditions : 

With walls, we have : 
Here, we see a few bumps in the kinetic energy corresponding to the collisions with the walls.

Then, the same graphs for the rotational kinetic energy : 

When there is a ceiling, we obtain : 
and : 

We notice, in particular, that at first, the particle picks up some rotation, but then it almost loses it entirely, when it hits the ceiling (notice that the time axis has been zoomed in when compared with the other figures, in order to show more details).

The next experimental values, we can take a look at, are the probability distribution function (PDF) of the velocities. We can notice that the PDF of hte velocities look a lot like Gaussians, whose standard deviation we can measure using scipy. As an example, here is the PDF of one of the first experiments with $\sigma=0.01$ and $\phi=0.2$, plotted with a fitted Gaussian and a "reference" Gaussian whose standard deviation is given by the strength of the kick $\sigma$ : 

To have a general feeling of the difference between all of the experiments, we can plot the effective $\sigma$ computed with scipy as a function of the packing fraction.

We have thus for our first set of experiments : 
We find much like our look of the kinetic energy, that as we increase the pakcing fraction, the overall behavior of the system for different low-valued $\sigma$ tend to converge.

For our second set of experiment : 
As before, we find that the behaviors of low kicks experiments converge much sooner than our first batch of experiment (we however note that the function curve_fit form the module scipy tends to change wildly its value according to its "first guess", even though this is not shown in the error bars computed by this same function).

We have the same observation for our last set of experiment : 

Finally we can look at the Mean Square Deviation (MSD) as a function time for different experiments. As we would expect, we can see a ballistic phase and a diffusion phase for most experiments. For example, this curve from the first batch with $\sigma=0.1$ and $\phi=0.1$, (actually I inadvertently plotted the squqre root of the MSD): 

However, for this first batch experiments, there were also some stranger curves. One example of these strange behaviors are for high packing fraction and/or kick (whose explanation might lie in the fact that the particles are confined by a crystalline structure for high packing fractions, and by the box for high values of $\sigma$). For example, here we have $\sigma=0$ and $\phi=0.8$ : 

Another type of weird behaviors happen for low packing fraction and no kick, here we have $\phi=0.1$ (I have no idea what happens in this one) : 

## Future questions
- Theory for underdamped behavior?
- Can we observe osmotic-like pressure in mixtures with semi-permeable membrane?
- Can we define a bulk pressure? (Also check literature)
- BONUS: Do we see random organization similar to what discussed [here](https://www.dropbox.com/sh/440k3glywzm3e9x/AACfGf9MqgmDdestTGn1D8AKa?dl=0)

- Analytical treatment? Overdamped is manageable in the dilute limit.
- Test whether there is a simulation regime where the behavior is captured by an overdamped description.
- As training, try to rederive Nature results.
- Density profiles near wall.

## References