Cylindrical triaxial test with flexible membrane


The task of this study was to investigate the role of support friction in triaxial tests using the discrete element method. Spherical grains and elastic membrane around the cylindrical specimen were used. Parametric study was performed varying specimen size, aspect ratio and support friction (normal or reduced values). The conclusion was that the role of support friction is in the same order of magnitude as scatter due to random factors in the simulation (such as particle arrangement), in contrast to laboratory experiments. A likely cause for this outcome is are spherical grains which do not have any rolling friction, unlike irregular shapes of real soils.


The simulation presented here was done using the woo.pre.cylTriax.CylTriaxTest preprocessor. The results can be obtained by downloading


woo-batch --job-threads=3 cyl-triax.batch.xls cyl-triax.preprocessor

(adjust --job-threads according to how many cores you want to allocate for each simulation; or don’t specify for using one core per simulation) and waiting a few hours for simulations to finish will result in reports for individual simulations being created, and also the cyl-triax.batch.hdf5 file with aggregate results. Afterwards, run:


to produce graphs shown below.


This project was supported by a commercial contract. We honor the wish to not disclose the identity of the contractor.


The simulation runs in coordinate system where \(z\) is the axial direction, \(x\) and \(y\) are radial; average quantities, such as radial stress, are noted with the \(rr\) index, as in \(\sigma_{rr}\equiv(\sigma_{xx}+\sigma_{yy})/2\).

The simulation runs in three distinct stages:

  1. Compaction. Non-overlapping spherical particles are generated in cylinder-shaped volume, of which radii satisfy the PSD (particle size distribution) function. Cylindrical volume, with rigid boundaries, is compressed to reach isotropic stress value \(\sigma_{\rm iso}\); friction parameters are set to zero to ensure high compacity of the resulting packing, as rearrangements are not constrained by tangential forces. Stresses are evaluated from global stress tensor \(\sigma\) computed from contact forces. The controlling logic is independent for each axis and uses damped feedback loop [Mod13] to reach desired stress state \(\sigma_{rr}=\sigma_{zz}=\sigma_{\rm iso}\).

  2. Membrane stabilization. Rigid cylindrical boundary is made compliant (it is transformed into membrane), and each membrane element is loaded with surface load equal to \(\sigma_{\rm iso}\). Membrane surface load is adjusted so that \(\sigma_{rr}\) reaches \(\sigma_{\rm iso}\), and \(\eps_{zz}\) deformation is allowed to maintain \(\sigma_{zz}=\sigma_{\rm iso}\). Friction angle of particles/membrane/supports is incremented during this phase linearly to reach real levels used in the triaxial phase.

  3. Triaxial phase. The initial triaxial configuration is used as reference configuration with zero strain tensor. The value of \(\dot\eps_{zz}\) is increased linearly to reach prescribed value, and axial compression continues until prescribed axial strain \(\eps_{zz}^{\rm max}\) is attained. During this phase, volume enclosed by the membrane+supports triangulation is monitored. Besides axial strain \(\eps_{zz}=\log\frac{l_z}{l_{z0}}\), volumetric strain \(\eps_v=\log\frac{V}{V_0}\) and stress tensor \(\sigma_{ii}\), we can compute mean and deviatoric stresses

    \begin{align} p&=\frac{1}{3}\tr\tens{\sigma}, & q&=\sigma_{zz}-\sigma_{rr} \end{align}

    and radial and deviatoric strain

    \begin{align} \eps_{rr}&=\frac{1}{2}(\eps_{v}-\eps_{zz}), & \eps_d=\eps_{zz}-\frac{2\eps_{rr}+\eps_{zz}}{3}. \end{align}

    These are used to plot the \(q/p(\eps_d)\) and \(\eps_v(\eps_d)\) functions.

Input parameters

Particle parameters were using usual values for stiffness and the \(k_T/k_N\) ratio, with internal friction angle \(\tan\phi=0.4\). This value was valid for inter-particle contacts, and for particle-membrane contacts. Supports had the same parameters as particles (stiffness), but their friction was reduced in to 4° in some simulations. Dimensions of the cylinder had varying aspect ratios; number of particles was around 5000 for all simulations, with the exception of the 1:1A one, which had approximately 14.000 particles. Membrane parameters were taken roughly from [ML81].

Density of particles was scaled in order to establish quasi-static conditions with reasonable amount of computation [TA98].

All the parameters are shown in the following; quantities marked with “†” are scaled for the purposes of the simulation.

Membrane parameters

Young’s modulus



density †





friction angle \(\phi\)







Particle parameters

Young’s modulus



density †





friction angle \(\phi\)







(see PSD)

Support parameters

Young’s modulus





friction angle \(\phi\)

a, A








Global parameters

hydrostatic stress \(\sigma_{\rm iso}\)



maximum axial strain



cylinder height×diameter
















Particle size

Fig. 30 Prescribed and real PSD for the 2:1 simulation


Each simulation resulted in \(q/p(\eps_d)\) and \(\eps_v(\eps_d)\) curves. They compare simulations with every feature varied:

Number of particles

The simulation was done with either the usual number (a/b) or an increased number (A) of particles:


Fig. 31 Influence of the number of particles on the result – 4.000 particles vs. 12.000 particles (larger specimen).

Support friction

\(\phi\) was either 4° or 22° − those simulations were respectively denoted with a/A or B in their simulation title.


Fig. 32 Cylinder 1:1, with normal and reduced support friction.


Fig. 33 Cylinder 1.5:1, with normal and reduced support friction.


Fig. 34 Cylinder 2:1, with normal and reduced support friction.


Fig. 35 Cylinder 3:1, with normal and reduced support friction.

Specimen aspect ratio

The height/diameter was 1:1, 1.5:1, 2:1 and 3:1.


Fig. 36 All simulations with reduced friction − influence of specimen shape.


Fig. 37 All simulations with normal friction − influence of specimen shape.


The influence of reduced support friction seems to be negligible – it is of the same order of magnitude as scatter of simulation results shown in tab. Influence of the number of particles.

The most obvious explanation for this discrepancy between our simulations and experiments is the shape of particles.

  • Spherical particles as a part of force chains have, by virtue of rotational symmetry, the property of not transforming shear force to normal and vice versa. In another words, shear force leads to rotation, and rotation only causes shear force (supposing small displacements).

    As there is a dissipation mechanics (Coulomb friction) in the shear direction, the amount of stored elastic energy, originating from shear force anywhere in the system, is globally limited.

  • On the other hand, with non-spherical particles, rotation around the centroid may lead to increment of normal force with some other particle, as the contact direction does not intersect the centroid. It follows that shear loading of particles at the support may be transmitted much further into the specimen, thus having substantial influence on the global behavior.

A possible remedy for future simulations is to

  • introduce non-physical rolling stiffness at contacts, accounting for non-spherical grain effects and micro-asperities of grain surfaces, i.e. macro-level and micro-level interlocking;

  • use non-spherical particle shapes, such as clumped spheres.



Fig. 38 Membrane mesh during simulation (front view).


Fig. 39 Membrane mesh during simulation (front view); triangles show the undeformed (reference) configuration.


Fig. 40 Cut through the specimen mid-plane.


Fig. 41 Cut through the specimen mid-plane (horizontal displacements scaled 3×).


Fig. 42 Force chains transmitting normal loads (displacements scaled).


Fig. 43 Node traces (top view); membrane nodes on the perimeter (displacements scaled).


Fig. 44 Node traces (lateral view) clipped in the vertical mid-plane, without membrane (radial displacements scaled, axial displacements not shown).


Fig. 45 Relative rotation of particles; specimen clipped in the vertical mid-plane.


  • Global view of the simulation (displacements scaled):

  • Force chain evolution during the compaction phase (displacements scaled):


Report issues or inclarities to github.