2D shearing simulations and particle density

Asked by Kristine Jarsve on 2017-05-16


I'm trying to run some 2D shearing experiments and I've run into some problems with the particle density I don't really understand.

My primary units are:

units of length, [L] = 1.0e-3 m;
units of time, [T] = 1 s;
units of force, [F] = 1 N.

This gives the conversion

units of velocity, [V] = [L]/[T] = 1.0e-3 m / 1 s = 1.0e-3 m/s = 1 mm/s
units of stress, [s] = [F]/[L]^2 = 1 N / (1.0e-3 m)^2 = 1.0e6 N/m^2 = 1.0e6 Pa = 1 MPa
units of mass, [M] = [F] [T]^2 / [L] = 1 N * (1 s)^2 / (1.0e-3 m) = 1.0e3 N.s^2/m = 1000 kg

The model has a size of 1000x1000 model units with a predefined fault at y = 500 model units. The particle radii are in the range [1,5] model units.

I've created 3 different samples where 2 corresponds to sandstones and 1 to a granite (all samples have the same geometry) by running some uniaxial compression tests. I had some problems during this stage as well, but resolved those problems by setting the particle density = 1 using

mySim.setParticleDensity (
  tag = 0,
  mask = -1,
  Density = 1

During my shearing simulations I apply a normal stress on both walls of the sample, compressing it until it reach the desired normal stress. After it has resettled I start shearing it while still applying the normal stress. The shearing slowly ramps ut before hitting its full velocity. There should be no gravity in my simulations.

The problem:

The shearing simulations works fine for the sandstone samples, but not for the granite. When I try to load the granite sample the bonds quickly disintegrates. When I tried to set the particle density = 2 however, it worked.

Is there a reason for this? How exactly does the particle density affect all of this? How will a particle density = 2 compared to 1 affect my results? (As mentioned above, there's no gravity (or at least it shouldn't matter).)

I've read (on this site) that the default particle density is 1 model unit. When I use setParticleDensity and set the density to one, does that even change anything?

The scripts I've been using can be found here: https://www.dropbox.com/sh/2nrxhhqgkzn2bwa/AAAefHIL62KuZofGC-FwOhfqa?dl=0

The particle density has been set to 2 in the Granity.py script. In the Sandstone.py script it is 1.

Thank you so much for your help! As you can see I'm rather confused by this!


Question information

English Edit question
ESyS-Particle Edit question
No assignee Edit question
Solved by:
Last query:
Last reply:
SteffenAbe (s-abe) said : #1

Hi Kristine,

from your description my first guess (without looking at the actual script) would be a timestep issue. Reducing particle density increases the wave speeds in the model, and therefore reduces the maximum stable time step. OTOH triggering this effect by changing the density by a factor of 2 sounds unlikely.
Anyway, can you try re-running the density=1 model again with the time step size reduced by a factor of 2?


Qi Shao (uqqshao) said : #2

Hi Kristine,

As we have discussed in a previous thread (https://answers.launchpad.net/esys-particle/+question/552493), try to use the following values in your simulations and see how it goes:


mySim.setParticleDensity (
  tag = 0,
  mask = -1,
  Density = 2500.0/1.0e12*1000


Hi guys, thanks for your answers.

Steffen: It seemed to work (at least for the one simulation that's completed). Any idea why/how (and why wasn't this a problem in the other runs)?

Also, setting the density to 1 model unit means that my particles have a density of 10^12 kg/m^3. I am aware that this is reasonably high, but how much will my results be influenced by that?

Qi: I probably should have tried that immediately, but then one of my supervisors didn't agree.
If I'm to use that time step for my shearing simulations, it would cause the number of time steps to be 1992081873, meaning the simulations would take between 3458.47 and 4611.3 days with my current set up. That's simply not an option.


Density should be 10^9 (rho = M/L^2 = 10^3/((10^-3)^2)).


SteffenAbe (s-abe) said : #5

Hi Kristine,

the reason why the runs don't work if the time step is set too large in relation to density, stiffness and particle size is that the time integration scheme used in the code becomes unstable when the timestep is larger than the time needed by an elastic(1) wave to travel between the two closest particles (in terms of center distance) multiplied by some safety factor (Courant-condition). The safety factor is usually around 0.2.

(1) _normally_ in a linear elastic material the compressional (P-) wave is fastest, so the relevant wave speed is on the order of sqrt(E*\rho) - give or take some 20% depending on the Poisson's ratio of the material. However, because a DEM with rotational particles is not a linear elastic material in the strict sense but rather a micropolar material (i.e. you have local rotational degrees of freedom) you also have rotational waves which, under certain circumstances, can be faster than P-waves. I don't have an equation handy to calculate their speed, but I've encountered cases where I had to reduce the time step significantly below the Courant condition to get a simulation to run stable.

As for the number of time steps needed for the shear model: Welcome to the club! This unfortunate problem affects most shear simulations I know. You've got two ways to get abound this restriction:
(A) shear your model faster than the real process you're modelling (cf. Mair & Abe, EPSL 2008). This should be fine as long as you stay well within the quasi-static regime (i.e. kinetic energy << elastic energy)
(B) reduce the wave speed by increasing the particle density to unnaturally high values - the so called "density scaling" (don't have a reference handy). Same caveat as above.


Hi Steffen,

In regards to the "density scaling", I thought that was what I had done? By setting the particle density = 1 model unit (which should be 10^9 kg/m^2 with my unit conversion), the particle density is unnaturally high. Extremely high I would say.

The granite sample does have a significantly larger E than the sandstone sample (17500 MPa and 6600 MPa respectively). So by your equation of the wave speed, the granite sample will have wave speeds x2 the wave speed of the sandstone. Is that why it worked reducing the time step with a factor 2 for the granite shearing simulations?


Best SteffenAbe (s-abe) said : #7

Hi Kristine,

what you have done with the high density is indeed "density scaling". I was pointing out the two options available.

W.r.t. the granite vs. sandstone situation it looks like your timestep was just on the edge of stability, so the relatively small difference in the P-wave speed made the difference between stable (sandstone) and unstable (granite). And then reducing the timestep for the granite model shifted it back into the stable region.

Btw, isn't E=6.6GPa rather low for sandstone? Literature (Zoback, "Reservoir Geomechanics") suggests 40-60GPa for most intact sandstones.


Hi Steffen,

Thank you for explaining! Things are a bit more clear now.

E = 6.6 GPa is indeed rather low for a sandstone. I forgot to say it was the Young's modulus of the bonds that was set to 6.6 GPa. The calibration (uniaxial compression test) then gives me a bulk E = 19 GPa. (Structural Geology, Haakon Fossen) states that E of a sandstone is between 10-20 GPa so I aimed for a bulk E in that range. However, the fact that Fossen's suggestion is so much lower than Zoback is a bit worrying (and odd).


Thanks SteffenAbe, that solved my question.