Modelling granular soil behaviour using a physics engine

The discrete element method (DEM) has proved increasingly popular as a research tool in the field of geomechanics since its introduction to the community by Cundall & Strack (1979). However, it has been somewhat handicapped by long simulation run times, and the hope that improvements in computer power alone would solve this problem, as expressed by Cundall (2001), may have been misplaced. In fact processor performance gains have been slowing down with every new generation, and this trend is unlikely to change in the near future (Esmaeilzadeh et al., 2012). Significant reductions in DEM run times are possible by using high-performance parallel computing platforms (O'Sullivan, 2015), but the beneficiaries, at least in the short term, are primarily likely to be researchers with access to supercomputers. To address this it seems worthwhile to look at alternatives to DEM. In particular, the physics based animation tools used in the computer games and film industries, referred to as ‘physics engines’, appear worthy of investigation. Physics engines are software libraries which employ a wide variety of simulation techniques; considering rigid body simulation, a good overview was provided by Erleben (2005), and, more recently, by Bender et al. (2014). For games the simulation needs to be real-time; thus traditionally physics engines favour speed, robustness and stability over accuracy. There is, however, a demand for more physical realism in video games and simulation methods are continually being improved. A list of commercial and open-source physics engines is provided by Bender et al. (2012). Early attempts at modelling soil have been described by Izadi & Bezuijen (2015) using the open-source Bullet physics engine (Coumanns, 2012) and by Pytlos et al. (2015) using the open-source Box2D physics engine (Catto, 2011), both with promising results. The aim of this paper is to provide a more thorough verification of the physical correctness of the simulations that can be obtained using Box2D in the context of soil modelling; this will help to establish whether Box2D is a viable alternative to the two-dimensional DEM codes currently used in the geomechanics field.

Box2D is a two-dimensional physics engine which simulates the dynamic interaction between discrete bodies. The continuous motion of bodies is discretised in the time domain and the simulation progressed using a time-stepping scheme. Each time step can be viewed as a sub-problem, where the task is first to calculate the rate of change of movements, and then to update the variables describing the state of each body. Objects are idealised as rigid bodies and free body motion is governed by the Newton–Euler equations.

Simulation of an assembly of granular particles requires the equations of motion to be augmented in order to prevent bodies from inter-penetrating and to model friction; this is achieved by means of the contact model. Whereas a traditional DEM code based on the distinct element method (Cundall & Strack, 1979) uses a penalty based contact model, in Box2D a constraint based contact model is used, which is more akin to the ‘contact dynamics’ approach considered by Jean (1999) and Radjai & Richefeu (2009). This can be considered to be the main difference between Box2D and traditional DEM approaches.

Consider two bodies in contact as shown in Fig. 1(b). The state of each body i, which has mass m_i, is defined by the position of the centre of mass, x_i, translational velocity, v_i, rotation, Ω_i and rotational velocity, ω_i. The contact is defined in terms of the contact normal unit vector, $n^$⁠, contact tangential unit vector, $t^$⁠, coefficient of friction, μ, and, for both bodies, the position of the contact point, p_i, and the moment of inertia about the contact point I_i. The non-penetration constraint can be formulated as follows

If the bodies are in contact (C_n ≤ 0) the tangential constraint, which is used in conjunction with the Coulomb friction law to simulate friction between the bodies, can be defined as

In the distinct element method violation of the constraints, illustrated in Fig. 1(c), is allowed and can be considered to be an inherent feature of the contact model. In a given time step the violations are calculated for the current state of the system. Force–displacement laws are then used to add forces to the system in order to minimise the violations. The penalty method is easy to implement and theoretically very fast because of the simple computations involved. However, in practice, stability of the simulation necessitates the use of a very small time step size, so that run times can be very long. A second disadvantage of the penalty method is that the contact parameters do not represent real physical properties of the particles, and are instead ‘tuned’ in order to achieve the desired macro-level behaviour (O'Sullivan. 2011).

In Box2D, a constraint-based approach is used. For a pair of bodies in contact (C_n ≤ 0), the constraints are formulated at the velocity level

where $p˙i=vi+ωi(pi−xi)⊥$⁠. Assuming that, to a given tolerance, points p₁ and p₂ are coincident, that is, $(p1−p2)=0$⁠, then these constraints can be simplified to

In a given time step the solver first computes tentative, or ‘pre-impulse’, velocities $vi−$ and $ωi−$ for each body, assuming free body motion. For a pair of bodies in contact (C_n ≤ 0) these velocities are then tested against the tangential constraint (⁠$C˙t=0$⁠). The solver computes an impulse, $Jt=jtt^$⁠, which instantly changes the relative velocity $p˙1−−p˙2−⋅t^$ in order to prevent violation of the constraint. The magnitude of the impulse is calculated from

where r_i = p_i − x_i. The magnitude of the impulse is limited by Coulomb's friction law

where j_n is the magnitude of the impulse in the normal collision direction. If j_t calculated in equation (7) is outside the friction limit, its value is reduced and the tangential constraint will not be satisfied (the bodies will slide). The post-impulse velocities for each body, $vi+$ and $ωi+$⁠, are then calculated from

A similar process is carried out for the non-penetration constraint; that is, if $C˙n<0$ then the magnitude of the corrective impulse is calculated from

where e is the coefficient of restitution. The new velocities $vi++$ and $ωi++$ are calculated from

(if $C˙n≥0$ then j_n = 0).

The solver treats constraints sequentially and when multiple bodies are in contact, several iterations are required in order to converge to an accurate global solution. Once the correct velocities v_i(t + Δt) and ω_i(t + Δt) are established, the integrator is restarted and the new positions are calculated from

where x_i(t) and ω_i(t) are, respectively, the position and the rotation from the previous time step and Δt is the time step size. (Note that in Pytlos et al. (2015) it is incorrectly stated that Box2D uses an explicit Euler integrator.)

Although per time step the contact model implemented in Box2D is computationally more expensive than the penalty method, it is still relatively fast because of its sequential nature. At the same time, the time step size can be much larger than in the penalty method without the danger of numerical instability or unduly sacrificing accuracy. This means that the overall simulation time in Box2D can be potentially much lower than in a code based on the distinct element method, although direct comparison is beyond the scope of the present paper. Another advantage of the presented contact model is that, considering frictional soil, there is only one relevant contact parameter, μ, and this directly represents a physical property of the soil being modelled.

Body inter-penetrations due to numerical errors cannot be removed by the collision solver because the constraints are formulated at the velocity level. The position error correction method implemented in Box2D (Catto, 2014) is not fully based on physics and is a potential source of simulation inaccuracy. It is therefore important to keep the position errors small by reducing the time step size if necessary. Note, however, that in practice the minimum usable time step size is still much larger than that required when using the penalty based contact model.

Additional information on the simulation method employed in Box2D can be found in Catto (2006, 2009, 2014), and also directly in the freely available source code.

To verify the ability of the physics engine to accurately model an assemblage of bodies, a simulation involving biaxial compression of hexagonally packed discs was run, and then validated against the theoretical solution derived by Rowe (1962). The model, consisting of 32 discs, is shown in Fig. 2(a). The theoretical stress ratio for this problem is

where σ₁ is the major principal stress, σ₂ is the minor principal stress, ϕ_μ is the angle of particle surface friction, defined as $ϕμ=tan−1⁡(μ)$⁠, and β is the deviation of the tangent at the contact point between sliding discs from the direction of the principal stress. Following Rowe (1962), the confining stress was simulated by applying forces F₁ and F₂, of magnitudes adjusted for disc spacing, to the centre of mass of the outer boundary discs as indicated in Fig. 2(a). The test was strain controlled, with the deviatoric stress applied by the top platen moving vertically at a constant velocity while the major principal stress σ₁ was measured at the base of the sample.

The top platen and the bottom boundary were rigid and frictionless and the test was conducted under zero gravity. The results of the test on discs with ϕ_μ = 10°, the same friction as in the experimental validation tests conducted by Rowe (1962), are shown in Fig. 2(c). For convenience, the strain in this problem was defined as the ratio of the vertical displacement of the top row of discs, δ, to the initial length per row in the vertical direction L (illustrated in Fig. 2(a)). The failure mechanism is shown in Fig. 2(b). At δ/L of about 0·56, when the stress ratio dropped to below one, a slip occurred. At this point the sample was unstable under the confining stress, and the problem became increasingly dynamic. At lower strains, up to δ/L = 0·16, the accuracy of Box2D is excellent. Past that point the accuracy is slightly lower with the error on the stress ratio typically within a few percent. The peak stress ratio is 101·0% of the theoretical value (100·5% in terms of angle of friction). This accuracy can be compared with the results reported by O'Sullivan & Bray (2003), who modelled the same problem with the PFC2D and modified DDAD software codes; their simulations yielded a peak angle of friction of about 92% of the theoretical value.

A critical state type response, as defined by O'Sullivan (2015), is one of the fundamental response characteristics for granular materials. A biaxial compression of a loose and a dense system were simulated in order to verify whether Box2D could successfully capture this type of response.

Although simple and relatively computationally inexpensive to model, an assembly of discs is problematic as a representation of a real soil because it leads to excessive rolling of individual particles. In DEM more realistic shapes are usually achieved by combining multiple discs into clusters. Box2D provides an efficient algorithm for modelling convex polygon-shaped bodies and it was decided to take advantage of this feature in the present study. The particles were modelled as randomly shaped convex dodecagons. Fig. 3 shows an example of the particles used in the simulations. The length to width ratio of particles was set to 1·0 in order to limit the effect of the initial fabric on the results. Particles were of uniform size d_b, where d_b is defined as a diameter of a circle bounding the full extent of a given particle.

The accuracy of the simulation can be controlled by adjusting the time step size, Δt, and the maximum number of velocity iterations per time step available to the constraint solver, N_i. Additionally, for given values of Δt and N_i, accuracy is affected by the ratio of the force experienced by a given particle to its mass. A high ratio will lead to a high tentative velocity, and consequently to large velocity errors that then have to be corrected by the constraint solver (as target velocities will be close to zero in a quasi-static analysis). Therefore, in the biaxial compression test described in this paper the values of confining pressure and density could be selected with a view to maximising accuracy for a given run time, rather than to replicate laboratory test settings (i.e. as this test simply involves an assemblage of rigid particles under zero gravity, these parameters will not affect the physics of this quasi-static simulation).

The simulation input parameters are given in Table 1. Particles were created simultaneously at a random position and at a random orientation within a two-dimensional zone of approximately 50 × 166 d_b; the particles were not allowed to overlap. The sides of the zone were bounded by temporary, rigid, frictionless walls. Once created, the particles were allowed to fall under gravity. After all the particles had come to rest, the top of the sample was levelled, gravity was gradually reduced to zero and the temporary side walls were removed. The final sample size was 50 × 100 d_b. The initial sample density was controlled by the particle friction coefficient (μ_g) used at the sample generation stage.

The exact shape and drop position for each particle were determined by the random number generator; in order to check repeatability both test setups were run three times with different seeds. The initial void ratios for all specimens are listed in Table 2.

The confining pressure in the horizontal direction was simulated by applying forces to the centre of mass of the particles at the perimeter of the specimen. Boundary particles were determined automatically in each time step by casting multiple horizontal rays along the height of the specimen. The force applied to a boundary particle was proportional to the number of rays hitting the particle in a given time step. The confining pressure in the vertical direction was applied by a ‘servo-controlled’ rigid top cap. The confining pressure was applied incrementally until the sample was in equilibrium at 1 kN/m. The particle coefficient of friction was then gradually changed to the value used for shearing (μ_s).

The compression was strain controlled with deviatoric stress applied by the top cap moving vertically at a constant velocity. Fig. 4(a) shows the mobilisation of the angle of friction with the axial strain for both the loose and dense states (three simulations per state). The corresponding evolution of volumetric strain and global void ratio are shown in Figs 4(b) and 4(c), respectively. Figs 5(a) and 5(b) show particle arrangement and accumulated rotation at 15% axial strain of the seed 1 dense and the seed 1 loose samples, respectively.

Overall the results show reasonably good repeatability between simulations. Qualitatively all the dense and the loose samples display behaviour typical of that obtained in laboratory experiments. The loose samples contract and the dense samples dilate upon shearing, consistent with critical state behaviour. However, the global average void ratios plotted in Fig. 4(c) do not converge as expected. Figs 5(a) and 5(b) show that in the loose system, shear deformations are spread across most of the sample, whereas in the dense system shearing is localised in distinctive bands, with large parts of the sample undisturbed by shearing, and having local void ratios largely unchanged from the initial values. In Figs 5(a) and 5(b) a volume element (RVE) considered to be representative of the part of the sample that underwent shearing is indicated by a red circle. Each RVE contained approximately 10% of the sample volume and was positioned in the zone with the highest density of particles having high accumulated rotation. The void ratios in the RVEs of all samples at 15% axial strain are given in Table 2. The results show that parts of the samples that underwent shearing converged, or are very close to convergence, to a consistent critical state void ratio.

The mobilised angle of friction also converges to the critical state value ϕ_crit, of about 21·5°, for both the loose and the dense states. The dense samples exhibit peak strength behaviour, with a peak angle of friction of about 29°; the peak strength corresponds to the maximum rate of dilation. Considering that ϕ_crit of sands observed in laboratory experiments is in the range of 32–37° (Bolton, 1986), the simulated shear strength appears low. However, this can be attributed to the relatively simple soil model used in this study. It is widely accepted that particle angularity and eccentricity have significant influence on the soil shear strength (Cho et al., 2006) and that the particle size distribution is also relevant (Morgan, 1999). This means that it is possible to achieve higher shear strength in Box2D by using model properties which are more representative of that of a real soil. For example, Pytlos et al. (2015) observed a higher value of ϕ_crit (29° as opposed to 21·5° found here) when using particles of a different geometry and non-uniform size.

A physics engine is a software library used in the film and computer games industries to realistically animate a physical system. Box2D, a widely used open-source rigid body physics engine, has been shown to be capable of accurately simulating disc interaction dynamics. The accuracy achieved was superior to that achieved using the PFC2D and DDAD DEM codes quoted in the literature. It has been demonstrated that Box2D can successfully capture the critical state type response of granular media with particles modelled as randomly shaped polygons. This work suggests that Box2D is a viable alternative to the current two-dimensional DEM tools for modelling particulate media, providing a simple and intuitive physics model.

The simulation method presented here has two potential advantages over a traditional DEM code based on the distinct element method that make it worthy of further investigation. First, physics engines used in computer games are designed to enable fast real-time simulations and this has the potential to translate into corresponding speed benefits for large particulate systems. Second, unlike traditional DEM approaches, the contact model does not appear to require extensive tuning and therefore the soil macro-scale behaviour would only be controlled by particle shape, size distribution and coefficient of friction, all of which can be a direct representation of the physical properties of the soil being modelled.

The authors would like to thank the reviewers for their constructive comments on an earlier version of this manuscript. This work was supported by the Engineering and Physical Sciences Research Council, under grant number EP/l014489/1.

d_b

diameter of a circle bounding the full extent of a given particle

e

coefficient of restitution

I

moment of inertia about contact point

J_t

tangential corrective impulse

j_n

magnitude of normal corrective impulse

j_t

magnitude of tangential corrective impulse

L

initial length per row of discs in vertical direction

m

mass

$n^$

contact normal unit vector

p

position of contact point

$p˙+$

post-tangential impulse translational velocity of contact point

$p˙−$

pre-impulse translational velocity of contact point

r

vector from centre of mass of body to contact point

t

time

$t^$

contact tangential unit vector

v

translational velocity of body

v⁺⁺

post-normal impulse translational velocity of body

v⁺

post-tangential impulse translational velocity of body

v⁻

pre-impulse translational velocity of body

x

position of centre of mass

β

deviation of tangent at contact point between sliding discs from direction of principal stress

Δt

time step size

δ

vertical displacement of top row of discs

μ

coefficient of friction

μ_g

particle coefficient of friction during sample generation stage

μ_s

particle coefficient of friction during shearing stage

σ₁

major principal stress

σ₂

minor principal stress

ϕ_μ

angle of particle surface friction

ϕ_crit

critical state friction angle

Ω

rotation of body

ω

rotational velocity of body

ω⁺⁺

post-normal impulse rotational velocity of body

ω⁺

post-tangential impulse rotational velocity of body

ω⁻

pre-impulse rotational velocity of body