Discrete Element Model (DEM) User Guide
The Discrete Element Model simulates rock as an assembly of bonded spherical particles. Where a DFN generates a fracture network from statistics you supply, a DEM predicts one: you build a particle pack, cement it into rock, load it, and the fractures emerge from where the cement actually fails.
The DEM is under active development. It is usable end to end — a model survives a save, tests run and their results can be exported — but parts of it are incomplete. Please read Current limitations before planning work around this feature.
The Discrete Element Model branch only appears in the Interpretation tree when your licence includes the DFN feature. Without it the branch is not created at all, so there is nothing to right-click. If you cannot find it, check your licence before looking for the tree node.
Overview
A DEM in VRGS lets you:
- Build a packed assembly of spherical particles at a target porosity, positioned where the rock is
- Give the assembly rock properties — density, stiffness, friction, restitution, rolling resistance
- Cement it, installing bonds that carry force and moment between touching particles
- Run the simulation on the CPU or the GPU and watch it in the 3D view while it steps
- Load it in a virtual press: unconfined compression or tension, or a triaxial test at a confining stress of your choosing
- Read the stress–strain response live on the DEM Test chart, and export it as CSV
- Track damage as bonds break, permanently and irreversibly
- Export the resulting broken-bond population as a DFN with real dip and azimuth
The distinction that matters: bonds persist and break. The population of broken bonds is the fracture, and that is what the DFN export reads.
Creating a DE volume
Right-click Discrete Element Model in the Interpretation tree and choose New. VRGS creates a volume and immediately fills it with a random, non-overlapping particle pack using the current settings.
Select the volume to edit its properties in the Properties panel. Make Active marks it as the model the ribbon and the test chart act on when several exist.
Putting the model where the rock is
A model box has an Origin in project coordinates and a size. Fit To Active Outcrop sizes and positions the box over the active outcrop model so a DEM sits on the rock it represents rather than at the coordinate origin.
| Property | Meaning |
|---|---|
| Origin X / Y / Z | Minimum corner of the model box, in project coordinates |
| X / Y / Z dimension | Size of the model box |
| Show Domain Box | Draw the walls the solver enforces as a wireframe. A particle outside this box has escaped its container |
| Periodic X / Y / Z | Join the two faces on that axis: a particle leaving through one enters through the other, and pairs and bonds across the seam interact as if it were not there. For a long strike-slip or rift model, so the ends impose nothing. The walls on a periodic axis are not enforced |
| 2D Section (XY) | Model a two-dimensional section: discs in the XY plane, with the Z dimension of the box as the slab thickness. Particles are placed at mid-plane and move only in the plane, a disc weighs what a disc through the slab weighs, cement bonds have a rectangular section, and the Z walls are not enforced. Point gravity along -Y for a section. It reaches a 10 km section at a resolution 3D could not afford and compares directly with the 2D literature, but it cannot answer the 3D questions (along-strike linkage, relay ramps), so treat it as a complement to a 3D run |
Building the particle pack
| Property | Meaning |
|---|---|
| Initial Elements | How many particles to place when the volume is created or reset |
| Minimum / Maximum Size | Particle radius range. Particles are drawn uniformly between these |
| Target Porosity | Void fraction to aim for when adding elements |
A wide size range costs nothing in the neighbour search: both solvers keep one grid per size class (a class per doubling of radius, up to four) and search each with the reach it needs, so a model can be fine through the zone you care about and coarse elsewhere. Fault-zone width scales with particle size, which is why a graded pack is how a fault zone is resolved inside a model that could not afford to be uniformly fine.
Particles are placed randomly and rejected if they would overlap, so the achievable packing density is limited. Two commands help you densify:
- Add Elements — keeps inserting particles until the target porosity is reached or placement saturates.
- Grow Elements — expands each particle into the gap to its nearest neighbour. This raises density without adding particles, and is usually the more effective of the two.
Current Porosity is shown read-only and updates as you go. The domain ceiling follows a settling pack down, so porosity is measured against the volume the pack actually occupies.
Material properties
These describe the particles themselves, and apply whether or not the model is cemented.
| Property | Meaning |
|---|---|
| Density | kg/m³. Particle masses follow it |
| Young's Modulus | Elastic modulus, Pa. Sets contact stiffness and the stable timestep |
| Poisson's Ratio | Used by the timestep stability criterion |
| Contact Stiffness | Read-only. Derived from the modulus and particle size |
| Restitution | Bounciness on collision, 0–1 |
| Friction | Coulomb friction coefficient |
| Rolling Friction | Rolling resistance coefficient, 0 = spheres roll freely |
Contact stiffness is derived, not set directly. A stiffness quoted in N/m is only valid at one particle size — weight grows as the cube of radius while a linear spring does not grow at all — so a single fixed value would let large particles sink through the floor and through each other while small ones behaved correctly. Specifying Young's modulus instead keeps the model correct at any particle size, and it is the property you are more likely to know for a given rock.
Rolling friction gives smooth spheres the interlocking that real grain shapes have. Without it an assembly's friction angle stays near 20° whatever Friction is set to, and a broken shear band has almost no residual strength. Values of 0.1–0.3 are usual for sands and crushed rock. It is off by default. Both solvers carry it; on the GPU it costs a second history buffer the size of the slip history, allocated only while it is set.
Stiffer contacts are shorter, so they need a smaller timestep and the simulation advances more slowly in real time. Granite's ~65 GPa is physically honest but expensive. Reducing the modulus is standard practice in discrete element modelling when you care about bulk behaviour rather than contact detail — it exaggerates how much particles deform at their contacts, but leaves packing and flow broadly intact, and it can speed a model up substantially.
The contact model cannot deliver very inelastic collisions. Requesting 0.2 gives about 0.32 in practice, and the response floors out near 0.21 however low you set it. Values of 0.5 and above behave as expected; treat anything lower as "as inelastic as this model gets". To settle a pack, use Settling Damping rather than a low restitution.
Running the simulation
| Command | Effect |
|---|---|
| Start/Stop | Runs the simulation on a background thread. The 3D view updates live |
| Step Forward | Advances a single display step, for inspecting behaviour closely |
| Reset | Rebuilds the particle pack from scratch |
The view animates while the model runs. Progress appears in the Messages panel: elapsed time, RMS particle displacement, kinetic energy, contact count, and whether the run is on the CPU or the GPU.
The simulation runs as fast as the physics allows and redraws about thirty times a second, so how much simulated time passes per second of real time depends on model size and timestep rather than on any setting. Steps per Batch is the upper bound on how many physics steps run between checks for the next redraw; the default of 100 suits most models.
Settling Damping (default 0.7) is Cundall local damping: each particle is damped by a fraction of its own unbalanced force, so a dropped pack settles instead of ringing. It vanishes at equilibrium, so it does not change where the pack ends up, only how long the journey takes. It is a numerical device rather than a material property; set it to 0 to measure real dynamics such as a drop test.
CPU or GPU
Tick GPU Solver to run the physics on the graphics card. It takes effect when the simulation is next started and is roughly ten times faster than the CPU on large models. Loading tests run on either solver: on the GPU every wall's reaction is summed on the device each step and handed to the same servo and the same stress–strain bookkeeping the CPU uses, so a test gives the same curve whichever path ran it. VRGS falls back to the CPU, with a message, when there is no compute device or no Vulkan 3D view is open. The status line reports which path is actually running, so a fallback is never mistaken for a slow GPU.
On the CPU the contact and bond loops run across all cores once a model is large enough to repay it, around twenty thousand particles; smaller models use a cheaper single-threaded formulation.
Timestep
The solver chooses its own timestep from the material and particle size, using the Rayleigh criterion — the standard stability bound for DEM. Two read-only values are shown:
- Critical Timestep — the stability limit for the current model
- Timestep (s) — what the solver is actually using
Timestep Safety (default 0.1) is the fraction of the critical timestep to use. It controls
accuracy as well as stability: a collision lasts roughly 1.1 / safety timesteps, so 0.1 resolves
one in about eleven steps. Halving it doubles the run time and roughly doubles the resolution. If
you set a timestep larger than the critical value the solver overrides it rather than integrating
something that would blow up.
Run the model for a while before cementing. A freshly generated pack is loose and randomly placed; letting it settle under gravity produces a denser, more realistic contact network for the cement to bond.
Cementing the model
Cement Model installs bonds between every pair of particles that are touching or nearly so. This is what turns a loose granular pack into rock. Until you do it the assembly is cohesionless: it can rearrange, but it cannot fracture.
| Property | Meaning |
|---|---|
| Bond Radius Multiplier | Bond cross-section radius, as a fraction of the smaller particle |
| Bond Modulus | Elastic modulus of the cement, Pa. Sets the bond stiffness for the particle sizes in the model |
| Bond Stiffness Ratio | Bond shear stiffness over bond normal stiffness |
| Bond Stiffness | Read-only. The stiffness per unit bond area the modulus gives, for the largest particle |
| Tensile Strength | Bond fails in tension above this stress, Pa |
| Cohesion | Shear strength at zero normal stress, Pa |
| Friction Angle | Slope of the Mohr–Coulomb failure envelope, degrees |
| Strength Scatter | Coefficient of variation of bond strength, 0 = every bond identical |
Like the contact stiffness, the bond stiffness is derived from a modulus, so the cement is as stiff as it should be at any particle size. The bulk Young's modulus is the usual choice for the bond modulus. Projects saved before this change carry the explicit stiffness they were tuned on, and keep it.
The bonds carry both force and moment. The moment matters: without it particles rotate too freely and the model comes out far weaker in compression relative to tension than real rock.
Strength Scatter gives each bond its own strength when the model is cemented, drawn around the value you set. Real rock is heterogeneous, and that heterogeneity is what turns failure from one clean plane into diffuse microcracking that then localises — which is also what the DFN export is tuned to read. Values of 0.1–0.3 are typical. Cementing is reproducible: the same model and the same settings draw the same scatter.
Re-running Cement Model heals the model — it replaces the bond set rather than adding to it.
Bonds act in parallel with particle contacts, so their stiffnesses add and the stable timestep gets shorter. The solver accounts for this automatically, but a heavily cemented model will step more slowly.
Loading the model
Cementing gives the model strength; loading is what makes it fail. Set Mode in the Loading group and the test begins immediately; the DEM Test chart opens to follow it. Setting Mode back to Off stops the test and puts the containing box back.
| Mode | What it does |
|---|---|
| Compression | One face becomes a platen and advances at the strain rate while the opposite face holds. The four side walls are removed, so this is an unconfined compression test |
| Tension | The platen retreats instead. Only meaningful once cemented |
| Triaxial | Every face is first brought to the Confining Stress and the assembly allowed to settle there. Then the side walls keep holding that stress while the platen advances — a triaxial test |
| Isotropic | Every face is brought to the confining stress, and the test stops once the assembly has settled |
| Property | Meaning |
|---|---|
| Axis | Which axis the platen travels along |
| Confining Stress | Triaxial and Isotropic: the stress the servo-controlled faces hold, Pa |
| Strain Rate | Fraction of the starting length per second of model time |
| Max Strain | The shearing stage stops once this much strain has accumulated |
| Servo Alpha | How aggressively a face chases its target stress, 0–1. 0.5 is usual; lower is slower and more stable |
| Sample Interval | Physics steps between rows of the test history |
| Suggested Strain Rate | Read-only. The rate that keeps the test quasi-static; see below |
| Stage | Read-only. Which stage of the test is running: consolidate, shear, and so on |
Read-only measurements update as the test runs: Axial Strain, Volumetric Strain, Axial Stress, Lateral Stress, Deviator Stress (axial minus lateral), and Peak Stress. In an unconfined test on a cemented model the peak is its unconfined compressive strength.
Starting a compression test switches off the four side walls, leaving only the two platens. A specimen under axial load expands sideways, and that expansion is what puts the transverse bonds into tension — tension being how compression actually breaks rock. Hold the sides rigid at zero strain and no tension ever develops, while the axial bonds only get stronger under load, because their shear capacity rises with the stress pressing them together. A triaxial test holds the sides at a stress, not a position, which is a different thing: the specimen is still free to expand against it.
Running a triaxial test
- Cement the model.
- Set Confining Stress — a few percent of the strength you expect is a good first choice.
- Set Mode to Triaxial. The consolidate stage brings every face to the confining stress; watch Lateral Stress approach it and Stage change to shear.
- Watch Deviator Stress climb on the chart, then fall as the specimen fails. Peak Stress is the strength at that confinement.
- Repeat at three or four confining stresses. The peaks plotted against confinement are the failure envelope; a straight line through them gives cohesion and friction angle.
How to run a Discrete Element Model triaxial test walks through a complete test on a laboratory-sized specimen, from an empty box to a failure envelope.
How fast to load
A strain rate has to be slow enough that the specimen stays in equilibrium as it is squeezed, or what you measure is a dynamic impact response rather than a strength. Two read-only diagnostics say whether it is:
- Unbalanced Force Ratio — mean net force on a particle over the mean force a contact or bond carries. Below about 0.01 the assembly is in equilibrium.
- Inertial Number — strain rate × grain size × √(density / mean stress). Below about 10⁻³ the response is rate-independent, which is what quasi-static means.
Suggested Strain Rate is the rate that gives an inertial number of 10⁻⁴ at the confining stress, or at the current mean stress when no confinement is set. Slower than that is not more accurate, only more expensive: for millimetre grains the previous default rate was some twenty times slower than a test needs, which turned a quarter-hour test into hours.
The DEM Test chart
The DEM Test tab, docked beside the Properties panel, plots the test as it runs: axial, lateral and deviator stress in MPa against axial strain in percent, and volumetric strain in percent on the right-hand axis. Falling volumetric strain after the peak is dilatancy — the specimen bulking as it breaks.
The chart opens when a test is started and refreshes twice a second. It follows the model whose test you started, or the active model; Show Test Chart on a model's Simulation menu points it at that model. DEM Test Tab on the View ribbon shows or hides it.
Keeping the results
Every test records a row every Sample Interval steps: time, strains on all axes, the stress on each face, the averaged stress tensor, bonds broken in tension and in shear, and the two diagnostics. Export Test History on the model's Simulation menu writes it as CSV, with a header naming the columns, for a spreadsheet or a script. The history is also written beside the model's state file whenever the project is saved.
Reading damage
Three read-only values track failure:
- Bonds — how many were installed
- Broken Bonds — how many have failed. This population is the fracture
- Damage — broken bonds as a fraction of installed, 0–1
Breakage is permanent. A bond that fails never re-forms, even if the particles are pushed back together, because accumulated damage is the whole point.
Exporting fractures to a DFN
Export Fractures to DFN converts broken bonds into a fracture network and adds it to the Interpretation tree.
The export groups broken bonds into spatially connected clusters, fits a plane through each, and emits a fracture disc with the resulting dip and azimuth. Clusters that are not planar enough are rejected — rock accumulates diffuse microcracking as well as discrete fractures, and emitting a disc for every blob of damage would fill the DFN with features that do not exist.
Tensile and shear failures go into separate sets, so joints and shear bands stay distinguishable.
The Messages panel reports what happened, including what was rejected and why. If nothing is exported it tells you whether there were no broken bonds at all, or whether clusters were found but discarded as too small or insufficiently planar.
Once exported, the network is an ordinary DFN — see the DFN User Guide for connectivity, permeability and stress analysis.
Saving and reopening
A DE volume is saved with the project. Its parameters go into the project database; the particles, bonds and each contact's slip history go into a state file beside the project, and the test history beside that as CSV. A reopened model steps exactly as the saved one would have.
Cut sections
Add Section XZ / YZ / XY creates a cut plane through the volume. Sections currently appear in the tree but are not drawn in the 3D view; see the limitations below.
Current limitations
The DEM is experimental. These are known and specific.
Specimens are boxes. A triaxial cell here is a box with six independently controlled faces — a true-triaxial cell — rather than a cylinder with a flexible membrane. Rigid faces are the usual choice for calibration and strength envelopes.
Seeding from an existing DFN is not exposed. The engine can pre-break bonds that cross fractures from a mapped network, so a model can start with inherited discontinuities rather than as intact rock. This works but has no menu command yet.
Sections do not render. Add Section creates the tree entry and computes the intersection, but nothing is drawn.
Convert to Voxel Model produces an empty model. The command creates a correctly sized voxel grid and adds it to the project, but no voxels are filled.
Particles cast no shadows. The whole assembly is drawn in a single pass as camera-facing sprites shaded to look like spheres, which is what lets large models draw at all. Sprites do not take part in the shadow passes, so particles neither cast shadows on each other nor on the rest of the scene.
Very close zoom flattens the largest particles. Sprite size is capped by the graphics hardware. Zoom in far enough that a single particle fills a large part of the window and it stops growing with the others. This only appears at magnifications well past where a particle assembly is normally viewed.
A complete workflow
- Create a DE volume, and Fit To Active Outcrop if it should sit on a model.
- Run the simulation and let the pack settle under gravity, then stop it.
- Grow Elements if you want a denser pack.
- Cement Model to bond the settled assembly into rock.
- Set Loading Mode to Triaxial with a Confining Stress, or Compression for an unconfined test, and Run again. Watch the DEM Test chart.
- Export Test History to keep the curve, and Export Fractures to DFN once bonds have broken.
That gives you a strength, a stress–strain curve and a fracture network predicted from simulated rock failure rather than generated from statistics — which is the reason to use a DEM at all.
What the DEM is for
The DEM is most useful for exploring how a cemented granular assembly fails: how packing, particle size distribution and cement strength interact, what stress the assembly carries at a given confinement, and what the resulting damage looks like as a fracture network.
Treat the absolute numbers with care. Peak stress depends on your particle size, packing and cement parameters, and a discrete element model reproduces the character of rock failure far more reliably than it reproduces a specific rock's strength without calibration. The way to calibrate is the way a laboratory does: run unconfined and triaxial tests, compare the envelope with the rock you are modelling, and adjust the cement until they agree. Comparing models against each other is sounder than reading any single peak stress as a laboratory value.