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
- Start from a template for a sandstone core or a loose sand
- Confine the pack to an outcrop mesh rather than a box
- Give the assembly rock properties — density, stiffness, friction, restitution, rolling resistance — with different materials by layer if you need them
- 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, and watch the fractures form
- Follow how far the rock has moved from marker horizons, and put the results on a voxel grid
- 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: 1,000 particles of radius 0.3–0.9 in a 100 × 100 × 50 box at the coordinate origin. Those are starting values for an outcrop-scale model. Change them in the Properties panel, then choose Reset on the model's Simulation menu to rebuild the pack.
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, and shows the Discrete Element Model ribbon tab; choose it again to clear it. The commands on a model's right-click menu act on the models selected in the tree, whether or not they are active. The tree page lists every command.
Import on the same menu reads a model from a text file (*.txt): an element count and a second
number, then for each element its X, Z and radius and four whole numbers, which become the Lyr,
TimeBrks, Bonds and Count attributes. All the elements are placed at a Y of 0.
Starting from a template
New from Template on the same menu builds a model set up for a test rather than an outcrop.
| Template | Box and grains | Test settings |
|---|---|---|
| Sandstone Core - UCS | A 5-inch core, 12.5 × 12.5 × 25 cm once consolidated, of grains 1.5–3 mm in radius. Use Gravity off | Axis Z, Confining Stress (Pa) 1 MPa for the consolidation, Strain Rate (1/s) 0.5, Max Strain 0.02 |
| Sandstone Core - Triaxial 10 MPa | The same core | Axis Z, Confining Stress (Pa) 10 MPa, Strain Rate (1/s) 0.5, Max Strain 0.05 |
| Loose Sand - Settling Box | A 0.5 × 0.5 × 1 m box of grains 5–10 mm in radius. Use Gravity on, Young's Modulus (Pa) 100 MPa, Rolling Friction 0.2, Target Porosity 0.4 | None: start the simulation and the grains settle |
All three tick GPU Solver and leave Mode at Off, because choosing a mode is what starts a test. A core's section is square because the walls are a box; a strength is force over face area, so the shape does not change it. A core's elements go into a larger, loose box sized so that consolidation brings it down to the specimen. The Messages panel reports what was built and what to do next.
The core's grains and cement are aimed at wet Darley Dale sandstone (UCS 42 MPa, E 12 GPa): Young's Modulus (Pa) and Bond Modulus (Pa) 12 GPa, Rolling Friction 0.1, Bond Radius Multiplier 1, Tensile Strength (Pa) and Cohesion (Pa) 40 MPa, Strength Scatter 0.2. They are starting values for your own calibration, not a calibrated rock. To run a core:
- Set Mode to Isotropic and choose Start/Stop. The walls close in to the confining stress and the stage ends once the pack settles. This is the slow step.
- Stop the simulation, set Mode to Off, and choose Cement Model on the Elements menu.
- Set Mode to Compression for the UCS template or Triaxial for the triaxial one, and choose Start/Stop.
Scale Tensile Strength (Pa) and Cohesion (Pa) together until the peak stress matches your rock, and both moduli until the slope does. A bonded sphere model's tensile strength comes out at about UCS/3 to UCS/8, where sandstone's is nearer UCS/14, so calibrate to UCS and modulus first. For other points on the strength envelope, change Confining Stress (Pa) before step 1: 0, 5, 10, 20 and 40 MPa is the usual series.
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 mesh or point cloud, with a margin of 5% of its span on each side, so a DEM sits on the rock it represents rather than at the coordinate origin. The outcrop has to be active and visible: make it active on the Data tree first.
Fit To Active Outcrop empties the model. Particles placed for the old box would sit outside the new one, so the particles, bonds and any test go, and Reset fills the new box.
| Property | Meaning |
|---|---|
| Origin X, Origin Y, Origin Z | Minimum corner of the model box, in project coordinates |
| X dimension, Y dimension, 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. On by default |
| Periodic X, Periodic Y, Periodic 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. Off by default, and rarely wanted on Z, which gravity acts along |
| 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. Use Gravity pulls along -Y rather than -Z. 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. Off by default |
The origin and dimensions are in the Model Dimensions group; the other three are in Display Options.
On a 2D section, set Axis to X or Y before a test: the Z faces are not enforced, so a platen on Z presses on nothing. Triaxial then holds the one in-plane lateral pair at the confining stress, which makes it a biaxial test. Grow Elements still stops at the Z faces, so a disc cannot grow past the slab thickness.
Clipping the pack to the rock
Fit To Active Outcrop sizes a box, and random placement fills all of it, so a cliff-face model carries particles in the air in front of the face and every porosity is measured against a volume that is mostly not rock. Clip Pack To Outcrop confines the pack to the rock of the active mesh. Point clouds have no inside and are not used.
- A watertight mesh is used as a solid: the rock is inside it.
- An open mesh bounds the rock on one side, so VRGS asks which. It suggests the side behind the mesh's faces, because an outcrop surface faces the air it was photographed from. Yes keeps that side, No keeps the other and Cancel does nothing. If the triangles do not agree on a facing direction, the prompt says the suggestion is a guess.
Every element outside the rock is removed, and if any are removed all the bonds go with them, so clip before you cement. The Messages panel reports the elements removed and kept, the rock volume as a share of the box, and how much of the box the mesh covers. An open mesh says nothing beyond its own edge, so elements past it are kept: read the coverage figure before trusting the shape.
From then on Add Elements and Reset place particles only in the rock, and Current Porosity is measured against the rock's volume. Since Fit To Active Outcrop empties the model, the usual order is fit, clip, then Reset.
The bounds are stored by mesh name and rebuilt from the mesh when the project opens. If the mesh has gone, the model is left unbounded and a warning names it. Running the command again with another mesh active replaces the bounds.
Building the particle pack
| Property | Meaning |
|---|---|
| Initial Elements | How many particles to place when the volume is created or reset. Default 1000 |
| Minimum Size, Maximum Size | Particle radius range. Particles are drawn uniformly between these. Default 0.3 and 0.9 |
| Target Porosity | Void fraction to aim for when adding elements. Default 0.3 |
| Total Elements | Read-only. How many particles the model holds |
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. It costs elsewhere, though: above about 5:1 the smallest grains set the timestep while mostly carrying no load, and the largest need several times the contact history. 1.5:1 to 3:1 is usual, about 2:1 for bonded rock, and the Messages panel says so when a model goes past 5:1.
Particles are placed randomly and rejected if they would overlap, so the achievable packing density is limited. A request for more than half the box by volume is cut to what fits, with the same size range, and the Messages panel gives the numbers. Three commands help you densify:
- Add Elements — keeps inserting particles until the target porosity is reached or placement saturates, up to 20,000 particles each time you choose it.
- Prepare Pack — fills the box to the target porosity by depositing particles in lifts and letting each one settle. See below.
- Grow Elements — expands each particle into the gap to its nearest neighbour. This raises density without adding particles, but it widens the size range, because isolated grains end far larger than crowded ones: a pack placed at 5:1 can come out at 15:1. It only closes the gap to the nearest neighbour, so a loose cloud stays loose. If the Messages panel says it converged above the target, run the simulation to settle the pack, then grow again.
Current Porosity is shown read-only, and in the status line while the model runs. The domain ceiling follows a settling pack down, so porosity is measured against the volume the pack actually occupies.
Prepare Pack
Prepare Pack on the model's Elements menu fills the box to Target Porosity with a pack settled under gravity. Here porosity is measured against the whole box you set, so the solid volume it aims for is (1 − porosity) × box volume.
It needs a stopped model with Use Gravity on and Mode at Off, one that has not been cemented or clipped to an outcrop, and a gravity axis that is not periodic. It works on a 2D section.
- Any elements already in the model are deposited straight down onto the floor, keeping their radii and material sets. The clock, the contact history and any marker horizons are reset. An empty model starts with its first lift, so Initial Elements does not need setting.
- VRGS starts the simulation by itself and lets the pack relax.
- When kinetic energy has stayed low, it adds another lift above the pack, up to 5% of the box volume at a time, and relaxes again. The status line shows the fill as a percentage of the target and how close the motion is to the limits.
- When the target volume is reached and the pack has come to rest, the simulation stops and the Messages panel reports the element count and the porosity.
It stops with a message rather than claiming success if there is no headroom for another lift, if placement is saturated above the pack, or if the pack does not settle. Gravity alone cannot reach every porosity: a low target may need wall compaction or a larger model. If it does not settle, check Settling Damping. Changing the box, Target Porosity or the particle sizes during a fill also stops it.
Material properties
These describe the particles themselves, and apply whether or not the model is cemented.
| Property | Meaning |
|---|---|
| Density (kg/m3) | Particle masses follow it. Default 2650 |
| Young's Modulus (Pa) | Elastic modulus. Sets contact stiffness and the stable timestep. Default 65 GPa |
| Poisson's Ratio | Used by the timestep stability criterion. Default 0.25 |
| Contact Stiffness (N/m) | Read-only. Derived from the modulus and particle size, for the largest particle |
| Restitution | Bounciness on collision, 0–1. Default 0.5 |
| Friction | Coulomb friction coefficient. Default 0.5 |
| Rolling Friction | Rolling resistance coefficient, 0 = spheres roll freely. Default 0 |
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 the model by one physics step, for inspecting behaviour closely |
| Reset | Rebuilds the particle pack from scratch, from Initial Elements. Bonds, the test history, marker horizons and the elements' material-set assignments go with the old particles |
Stop the simulation before using the model's other commands: most are disabled or refused while it runs.
The view animates while the model runs. Progress appears in the Messages panel: elapsed time, element count, RMS particle displacement, kinetic energy, porosity, contact count, and whether the run is multi-threaded, single-threaded or on the GPU.
The simulation runs as fast as the physics allows and redraws up to 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.
Use Gravity (on by default) lets gravity act on the model: down Z in a box, down Y in a 2D section. Turn it off for a specimen whose self-weight is only noise against the test stress.
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, when no Vulkan 3D view is open, or when the model has more than eight material sets. 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 (s) — 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.
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, with surfaces within 5% of the smaller radius. 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. The Messages panel reports how many bonds went in; none means the pack is not touching.
| Property | Meaning |
|---|---|
| Bond Radius Multiplier | Bond cross-section radius, as a fraction of the smaller particle. Default 0.5 |
| Bond Modulus (Pa) | Elastic modulus of the cement. Sets the bond stiffness for the particle sizes in the model. Default 65 GPa |
| Bond Stiffness Ratio | Bond shear stiffness over bond normal stiffness. Default 0.4 |
| Bond Stiffness (Pa/m) | Read-only. The stiffness per unit bond area the modulus gives, for the largest particle |
| Tensile Strength (Pa) | Bond fails in tension above this stress. Default 1 MPa |
| Cohesion (Pa) | Shear strength at zero normal stress. Default 2 MPa |
| Friction Angle (deg) | Slope of the Mohr–Coulomb failure envelope. Default 30 |
| Strength Scatter | Coefficient of variation of bond strength, 0 = every bond identical. Default 0 |
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. A model saved with an explicit stiffness rather than a modulus keeps that stiffness.
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, so it also discards anything Seed From DFN did.
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.
Material sets
A model can hold several material sets: a mechanical stratigraphy, such as strong beds over weak ones or a cover over a rigid basement. Every model starts with one set, Set 0, named Base, whose properties are the Material and Cement groups. The commands are on the model's Material Sets menu, and the sets are listed in the Material Sets property group.
- Add Set adds a set as a copy of the base material as it stands, so it is a working layer at once. Under its row (Set 1, Set 2 and so on), edit Name, Density (kg/m3), Young's Modulus (Pa), Friction, Tensile Strength (Pa), Cohesion (Pa) and Bond Modulus (Pa). Elements shows how many elements are in the set.
- Assign Material Set... asks which set, then puts elements in it. It offers a depth band, as the lowest and highest Z in project coordinates. When the pack has been clipped to an open mesh it first offers one side of that surface instead, above or below it, and leaves alone any element the surface does not cover.
- Remove Set... asks which set to remove, the last by default. Its elements return to the base material and the sets above it move down one. The base material cannot be removed.
Colour the elements by the Material Set attribute to check the layers.
Where two sets meet, a contact takes the harmonic mean of their moduli and the lower friction, restitution and rolling friction. A bond takes the weaker set's cement, judged by tensile strength, on the harmonic mean of the bond moduli. So a strong bed over a weak one fails in the weak one, and an interface is no stronger than its weaker side.
Assign the sets before you choose Cement Model: the bond radius and the strength scatter are fixed when the bonds are installed. An assignment belongs to the element, not to its position, so it moves with the element. Reset builds new elements, all in the base material, so assign again afterwards; the sets themselves stay. Add Elements and new lifts from Prepare Pack also go into the base material.
The CPU solver takes up to 256 sets and the GPU solver up to eight. Past eight, a run falls back to the CPU with a message.
Loading the model
Cementing gives the model strength; loading is what makes it fail. With the simulation stopped, set Mode in the Loading group, then choose Start/Stop. Setting Mode lays the test out from the other values in the group and opens the DEM Test chart to follow it, so set those first: a value changed afterwards does not apply until the next test. Changing Mode while the simulation runs is refused, with a message.
Setting Mode back to Off stops the test and restores the containing walls at the positions the test left them, and the model's origin and dimensions take those positions. That is what lets an isotropic consolidation carry over into the test that follows it. Set Mode to Off before Reset, then set it again on the fresh pack.
| 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 (Pa) 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. Default Z |
| Confining Stress (Pa) | Triaxial and Isotropic: the stress the servo-controlled faces hold. Default 0 |
| Strain Rate (1/s) | Fraction of the starting length per second of model time. Default 0.05 |
| Max Strain | The shearing stage stops once this much strain has accumulated. Default 0.15 |
| Servo Alpha | How aggressively a face chases its target stress, 0–1. 0.5 is usual; lower is slower and more stable |
| Sample Interval (steps) | Physics steps between rows of the test history. Default 100 |
| Suggested Strain Rate (1/s) | Read-only. The rate that keeps the test quasi-static; see below |
| Stage | Read-only. Which stage of the test is running: consolidate, shear, isotropic compression, uniaxial compression or uniaxial tension |
The read-only measurements are Axial Strain, Volumetric Strain, Axial Stress (Pa), Lateral Stress (Pa), Deviator Stress (Pa) (axial minus lateral), and Peak Stress (Pa). In an unconfined test on a cemented model the peak is its unconfined compressive strength. The Properties panel reads these rows when you click the model in the tree, so click it again for current values; the chart follows the test live.
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 (Pa) — a few percent of the strength you expect is a good first choice.
- Set Mode to Triaxial and choose Start/Stop. The consolidate stage brings every face to the confining stress; watch the lateral stress approach it on the chart before the test moves on to shear.
- Watch the deviator stress climb on the chart, then fall as the specimen fails. Peak Stress (Pa) 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 (1/s) 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; it reads 0 when there is no stress to judge by. Slower than that is not more accurate, only more expensive: for millimetre grains at a megapascal, the default rate of 0.05 per second is more than ten times slower than a test needs.
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's title gives the model, the mode, the number of samples and the peak stress so far.
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, under Project Explorer in the Windows panel of the Home ribbon tab, shows or hides it.
Keeping the results
Every test records a row at each Sample Interval (steps): time, step number, stage, strains on all axes, the axial, lateral, deviator and mean stress, the stress on each face, the averaged stress tensor, bonds broken in tension and in shear, the two diagnostics, and kinetic energy. 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.
Watching fractures form
These are in the Display Options group. Particles, bonds and the domain box are drawn only in a Vulkan 3D view.
| Property | Meaning |
|---|---|
| Show Elements | Draw the particles. On by default |
| Show Bonds | Draw each intact bond as a line between the centres of the two elements it joins. Off by default |
| Show Broken Bonds | Draw each broken bond as a small disc across the bond where it broke: red for tensile failure, blue for shear. Off by default |
| Colour Attribute | What colours the elements: Sphere (a fixed colour per element, the default), Velocity, Density, Mass, Volume or Radius. An attribute ticked under the model's Attributes row takes precedence |
Turn Show Elements off and Show Broken Bonds on to watch fractures form while a test runs.
Element attributes
Each model has an Attributes row beneath it. Right-click it and choose Add Attribute to give every element a value computed from the model as it stands.
| Attribute | Value for each element |
|---|---|
| Radius, Volume, Mass | Its radius, volume and mass. In a 2D section the volume is that of a disc through the slab |
| Speed | Its current speed |
| Contacts | How many elements touch it: the coordination number |
| Intact Bonds, Broken Bonds | How many of its bonds are intact, and how many have broken |
| Bond Damage | The fraction of its bonds that have broken |
| Tensile Breaks, Shear Breaks | How many of its bonds broke in tension, and in shear |
| Material Set | The number of its material set, to colour the layers |
| Marker Band | The marker horizon it started in, coloured in discrete steps |
| Displacement | How far it has moved from its marked position |
| Vertical Displacement | Its vertical movement from its marked position, positive up |
All adds every attribute, leaving out the last three if no marker horizons are set.
Adding one attribute colours the elements by it. All leaves the colouring alone; tick a layer under Attributes to colour by it. A layer describes the moment it was added and does not follow a run, so stop the simulation and add it again to refresh it: that recomputes the same layer rather than stacking a copy. Attribute layers are not saved with the project.
Marker horizons and displacement
A fault model is read from where the rock started and how far it moved. Set Marker Horizons... records every element's current position as its reference, and bands the model into layers counted up from the floor of the box as it stands. It asks for the layer thickness, by default a tenth of the model's height: enough bands to read an offset, each many elements thick. The elements are coloured by Marker Band at once. The bands move with the elements, so a throw, a fold or an offset shows in the view as the test runs.
Marker Band, Displacement, Vertical Displacement and the displacement fields of Convert to Voxel Model... then read against this reference, and it is saved with the project. Without it these attributes read 0, with a warning. On a periodic axis, displacement is folded to the nearest image, so a particle that wrapped through a face has moved a little rather than a whole period.
Set the markers on the finished pack. Adding elements invalidates them, and Reset and Prepare Pack clear them. In a 2D section every element shares one Z, so Marker Band gives a single band and Vertical Displacement reads 0: use Displacement there.
Exporting fractures to a DFN
Export Fractures to DFN converts broken bonds into a fracture network, named after the model, 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 with the resulting dip and azimuth, shaped as an ellipse along the damage it fits. Clusters that are not planar enough are rejected — rock accumulates diffuse microcracking as well as discrete fractures, and emitting a fracture for every blob of damage would fill the DFN with features that do not exist.
Tensile and shear failures go into separate sets, DEM tensile (joints) and DEM shear, so joints and shear bands stay distinguishable.
The settings are in the Fracture Export property group.
| Property | Default | Meaning |
|---|---|---|
| Cluster Distance Ratio | 1.5 | Broken bonds within this multiple of the mean element diameter belong to the same fracture surface |
| Min Cluster Size | 6 | Smallest group of broken bonds that may become a fracture. Three is the least a plane can be fitted through; more suppresses isolated grain-boundary failures |
| Max Condition Score | 0.05 | How planar a group must be: near 0 is a sheet, near 1 a diffuse cloud. Groups above this are not exported |
| Polygon Segments | 12 | Edges of the polygon each fracture is drawn with |
| Clip To Domain | On | Clip exported fractures to the model box, except along a periodic axis, which has no face |
| Include Inherited | Off | Also export the fractures Seed From DFN put into the model, in a set of their own, DEM inherited (seeded). Off because they are the input, not a prediction |
| Elliptical Fractures | On | Export each fracture as the ellipse of its broken bonds' spread in the plane, long axis along its elongation. Off gives the disc that circumscribes them, which can give an elongated shear band up to twice its area. The radius is that of the disc of equal area, so set statistics stay comparable with a stochastic network's |
| Tensile Set Threshold | 0.5 | Fraction of a fracture's bonds that must have failed in tension for it to go into the tensile set. 0.5 classes by majority; 1 demands every bond, so one shear-failed bond at a joint's tip reclassifies the joint |
| Min Break Step, Max Break Step | 0 | Export only the bonds that broke within this window of solver steps, inclusive; 0 is no limit on that side |
To follow growth, export once per window. The step column of the test history
gives the solver step of each point on the curve. Seeded bonds have no step and are not windowed.
Each export also writes a table of its fractures as CSV beside the model's state file, in the
project's DEMODELS folder, replacing the last one. One row per fracture gives the identifiers that
match it to the DFN, the centroid, dip and dip azimuth, the radius, the ellipse's semi-axes with the
trend and plunge of its long axis, the area, the planarity score and residual, the half width of the
damage zone, the tensile fraction, the bond count, and the first and last break step. The damage
width and break steps have nowhere to go in the DFN itself.
In a 2D section each cluster is read as a trace. The fracture stands through the slab as a rectangle with the trace's true strike, and Max Condition Score judges how straight the trace is. On a periodic axis, a fracture across the seam comes out whole.
The Messages panel reports what happened, including what was rejected and why, and names the table. If nothing is exported it tells you whether there were no broken bonds at all, whether they all fell outside the break-step window or were all seeded, 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.
Seeding a model from a DFN
Seed From DFN takes its fractures from the active fracture network. Right-click the DFN in the Interpretation tree and choose Make Active first; it shows in bold. Only one DFN is active at a time, and with none active the command stops with a message saying so.
Seed From DFN on the model's Elements menu is the reverse of the export: it starts a model with the discontinuities mapped from the outcrop rather than as an intact block. Every intact bond that crosses a fracture of the active network is broken or weakened. Choose Cement Model first, because there is otherwise nothing to break, and check the model sits where the network is. Choosing Cement Model again discards the seeding.
| Property | Default | Meaning |
|---|---|---|
| Break Outright | On | Break every bond a fracture crosses. Off keeps the bond and scales its strength by Weakening Factor, for a healed or partly cemented discontinuity that can still fail under load |
| Weakening Factor | 0.1 | Strength multiplier for crossed bonds when Break Outright is off, 0–1. 0.1 leaves a tenth of the cement's tensile strength and cohesion |
These are in the Fracture Seeding property group. Bonds broken this way count towards Broken Bonds and Damage, and draw as tensile. They are left out of the next export unless Include Inherited is on, so that a seeded model does not hand its input back as a prediction.
Converting to a voxel model
Convert to Voxel Model... puts the model on a regular grid that the voxel tools can slice and compare with mapped surfaces. It asks for a cell size, by default twice the mean element diameter, so a cell holds a handful of elements. The grid covers the domain box and is coarsened if it would pass 50 million cells; a 2D section gets one layer of cells through the slab.
The voxel model is named after the DE model and added to the Data tree. See CT Image Stack and the Voxel Volumes guide for what to do with it.
| Attribute | Value for each cell |
|---|---|
| Solid Fraction | The share of the cell's volume inside an element, from each element's exact overlap with the cell |
| Broken Bond Density (per m^3) | Bonds that broke in the cell, per unit volume, seeded ones included |
| Displacement (m) | The mean distance moved by the elements centred in the cell. Needs marker horizons |
| Vertical Displacement (m) | The mean vertical movement of the elements centred in the cell. Needs marker horizons |
A layer that would be zero everywhere is left out: no broken bonds, no marker horizons, or nothing has moved. The Messages panel lists the layers written and those left out, with the reason, and any elements or bonds outside the box that were not counted. An empty cell is null, and so transparent. Local stress is not on the grid, because the solver keeps no per-contact forces to build it from.
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, and a test in progress picks up where it was when you next start the simulation. Material sets, marker horizons and the bounds from Clip Pack To Outcrop come back too. Attribute layers and cut sections do not.
Cut sections
XZ Plane, YZ Plane and XY Plane on the Add Section menu create a cut plane through the centre of the volume. Sections 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.
Sections do not render. Add Section computes the intersection, but nothing is drawn, a tree row appears only when the tree is next rebuilt, and the section is not saved with the project.
One mesh bounds a pack, and it cannot be removed. Only one mesh can be active at a time, so Clip Pack To Outcrop bounds a pack by one solid or one surface, and running it again replaces that one rather than adding to it. No command removes the bounds, and Prepare Pack refuses a bounded model, so start a new model if you need a plain box again.
Driven particles have no command. The solver can move a particle at a prescribed velocity, ignoring the forces on it, which is the basis of a rough floor, an indenter or a rigid basement block. Nothing in the application sets one.
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 with New, and Fit To Active Outcrop if it should sit on a model. Use Clip Pack To Outcrop to confine the pack to the rock, then Reset. For a laboratory specimen, use New from Template instead.
- Start the simulation and let the pack settle under gravity, then stop it, or use Prepare Pack to fill a box to Target Porosity.
- Grow Elements if you want a denser pack.
- For a layered model, Add Set and Assign Material Set....
- Cement Model to bond the settled assembly into rock.
- Set Marker Horizons... if you want to read displacement.
- Set Loading Mode to Triaxial with a Confining Stress (Pa), or Compression for an unconfined test, and choose Start/Stop. Watch the DEM Test chart, and turn on Show Broken Bonds to see the fractures form.
- Export Test History... to keep the curve, Export Fractures to DFN once bonds have broken, and Convert to Voxel Model... to put the result on a grid.
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.