Discrete element simulation of railway ballast: modelling cell pressure effects in triaxial tests

The paper investigates reproducing the effects of confining pressure on the behaviour of scaled railway ballast in triaxial tests in discrete element models (DEM). Previous DEM work, using a standard Hertzian elastic contact law with an elastic–perfectly plastic tangential slip model, has been unable to replicate the behaviour observed in laboratory tests across a range of confining pressures without altering both the material stiffness and the inter-particle friction. A new contact law modelling damage at the contacts between particles is introduced. Particle contact is via spherically-capped conical asperities, which reduce in height if over-stressed. This introduces plasticity to the behaviour normal to the contact surface. In addition, the inter-particle friction angle is varied as a function of normalized contact normal force. At relatively low normal forces the friction angle must be increased for peak mobilized friction angles to match the laboratory data, an effect that is attributed to interlocking at the scale of surface roughness. Simulation results show close agreement with laboratory data.


Introduction
Triaxial test results on railway ballast and other granular materials (e.g., [3,9,13,17,20,25]) have shown that, when other factors are held constant, the peak mobilized strength and volumetric strain are strongly influenced by the con-B John Harkness J.Harkness@soton.ac.uk 1 University of Southampton, University Rd., Southampton SO17 1BJ, UK fining cell pressure applied. In particular, at low confining pressure the peak mobilised shear strength is significantly greater. This has particular relevance in the context of railway ballast. Railway ballast is placed as a near surface layer and operates entirely within a low horizontal confining stress regime, and is required to perform acceptably over a large range of vertical to horizontal stress ratios.
Prior discrete element modelling (DEM) studies have modelled the effects of confining pressure on the behaviour of granular matter in different ways. Grain crushing has been modelled explicitly by various means (see, for example [4]). Other studies have modelled pressure effects without modelling crushing (see [19,28] for sands, [15,22] for direct shear tests of railway ballast and [26] for a qualitative study) by using constant-stiffness linear (rather than Hertzian) contact laws, which allow greater particle overlap at higher confining pressures and hence a softer behaviour. In this paper, it is demonstrated that for triaxial tests on a scaled granite ballast material, although no particle breakage was observed in the laboratory, some contact damage must be taking place with increasing confining pressure. This is consistent with anecdotal experimental evidence of ballast attrition observed in the laboratory. Furthermore, it is shown that the observed behaviour, including the correct loading/unloading response, cannot be modelled using a standard elastic contact law and a new model is proposed to address this need.
Ahmed et al. [2] presented discrete element simulations of railway ballast that showed close agreement with monotonic laboratory triaxial tests carried out by Aingaran [3], on onethird scale ballast. The experimental test rig is shown in Fig. 1. The numerical model used a Hertzian contact law with elastic behaviour normal to the contact plane and elastic/perfectly plastic behaviour in the tangential direction. The ballast was tested at confining cell pressures of 15 and 30 kPa; example results for a cell pressure of 15 kPa are Fig. 1 Photo of one-third scale ballast specimen of height = 300 mm and diameter = 150 mm in a GDS Instruments [10] triaxial apparatus. The triaxial cell was mounted within a reaction frame and bolted to a pedestal seated by self-weight onto a bottom-up mechanical drive unit used to apply axial strain or stress  Fig. 2. However, if we now extend and apply the same DEM model to test an identical numerical specimen at 200 kPa, a significant difference between laboratory and DEM results can be seen, as shown in Fig. 3. In comparison with results at 15 kPa, the real ballast exhibits much more volumetric compression and a slower rise to a lower peak strength. In contrast, the simulated response at 200 kPa is  generally similar to that at 15 kPa, except for a slight increase in volumetric compression and, perhaps surprisingly, in peak strength. These results are summarized in Table 1. According to established soil mechanics principles, for tests on soil samples with the same initial void ratio, increased confining pressure will lead to a reduction in both the peak strength (unless the peak strength coincides with the critical state) and the specific volume at the critical state. There are two main possible causes of the difference in stiffness between the laboratory and simulated tests: 1. The model is too stiff or is otherwise incorrectly calibrated. Contacts between real ballast particles are not Hertzian, as the particles have rough, non-spherical surfaces, so there is some scope for uncertainty in the choice of stiffness magnitude. However, the 10 GPa shear modulus used for the Hertzian contact stiffness in the simulation is already considerably less than that of granite, which is in the region of 17-29 GPa [12] (based on a Poisson's ratio of 0.2). 2. The ballast in the laboratory is experiencing some sort of damage that increases with cell pressure. Visual and sieve analysis before and after triaxial testing revealed no discernible particle breakage and no measurable change in particle size distribution. This was also observed in tests on full scale ballast [1]. Therefore, if damage is For the one-third scale particles, D 50 = 13.5 mm occurring to the material used in these tests, it must take the form of very small-scale crushing or abrasion of contacting asperities, perhaps reducing the effective surface roughness.
Other factors, such as approximations in the modelling of the triaxial test itself are discussed by Ahmed et al. [2].
First, the possibility is explored that variation of the model parameters could improve the fit to the laboratory data over the range of cell pressures 15-200 kPa. The main model parameters that can be adjusted are: 1. Particle shape and size distribution (PSD). 2. Initial void ratio. 3. Particle contact stiffness. 4. Inter-particle friction angle.
The particles were chosen so that the distribution of particle forms and PSD (shown in Table 2) matched the real ballast [2,18], and the void ratio at the start of the test was the same as the void ratio measured in the laboratory. Furthermore, approximately the same number of particles (2800) were used in the real and simulated tests. The D 50 for the one third-scale ballast was 13.5 mm.
The same initial specimen of particles was used for all of the numerical studies presented in this paper (see Fig. 4). The simulated particle shapes are approximations, which will potentially give rise to some errors both in the simulated behaviour and in the measurement of void ratios in comparison with the real ballast. However, for the purposes of this study, the particle size distribution and initial void ratio are considered to be fixed. This leaves only the inter-particle friction angle and stiffness as variables.
The effects of varying the inter-particle friction angle, for a constant shear modulus of 10 GPa and a constant confining pressure of 200 kPa, are shown in Fig. 5. These results were obtained using the discrete element model and the specimen of potential particles [11] presented in [2]. Reducing the interparticle friction angle has little effect on the initial response, which remains stiffer than the laboratory test. Also, the critical state strength is not significantly altered, with all three simulations being slightly weaker than the laboratory test.  However, the peak strength increases with increasing interparticle friction angle and is accompanied by an increase in the rate of dilation. Thus it is clear that a change of interparticle friction angle alone cannot correct the fit to the laboratory data at the higher confining pressure of 200 kPa (it would also spoil the fit at lower confining pressures). Figure 6 shows the effects of varying the Hertzian shear modulus, G, for a constant inter-particle friction angle of 40 • . With a value of G = 1 GPa, the initial stiffness of the response is much improved. However, the peak strength, rather than reducing as desired, has increased along with the dilation rate. At G = 100 MPa the peak mobilized friction angle is reduced, but the initial response is too soft. Once again, it seems that the critical state response is not significantly affected.
By altering both the Hertzian shear modulus (G = 0.5 GPa) and the interparticle friction angle (30 • ), a reasonable fit can be obtained at a confining pressure of 200 kPa (Fig. 7). The response at 15 kPa, using the same low Hertzian   shear modulus (G = 0.5 GPa), but retaining the higher interparticle friction angle of 40 • , was evaluated and is also shown. At 15 kPa, the reduction in stiffness has affected the peak mobilized friction a little, but the volumetric behaviour and the strength at 200 kPa is quite well matched. Although this choice of parameters yields quite a close match for the monotonic tests, universally changing the interparticle friction angle to match the data at different confining pressures is not a workable solution for more general loading cases in which the confining pressure is not held constant. Furthermore, although there is some room for manoeuvre in the choice of stiffness (ballast particles are not spherical and particle contacts are therefore not perfectly Hertzian), the very low value used (about one fiftieth of the value for granite) is difficult to justify. This is borne out in a comparison of cyclic loading tests shown in Fig. 8, for the laboratory and simulation data. These results are summarized in Table 3. This shows that a much stiffer model (G of at least 10 GPa) is required to match the per-cycle deflection observed in the laboratory scaled triaxial test.
One explanation for the observed behaviour is that plastic damage occurs at contacts as they are loaded. Cavarretta et al. [5] showed that for coarse grained sand the interparticle contact stiffness included an irrecoverable proportion of deflection beyond which the deflection vs. stress plot showed Hertzian contact behaviour. This would make the contacts appear initially soft in monotonic loading (as the asperities are crushed), but stiffer as loading progressed and stiffen cyclic loading behaviour if the model was subsequently unloaded and reloaded.

Conical damage model
To address the inconsistencies observed between the measured and simulated specimen stiffness at different confining pressures, a new contact damage model is proposed. The basis of the model is an idealized asperity represented by a cone (Fig. 9). The sides of the cone make an angle α with the vertical and the cone has a rounded apex with radius of curvature, R. This asperity contacts an elastic half plane via a Hertzian contact with the rounded tip. A minimum radius of curvature R min is specified, as shown in Fig. 9a. The maximum stress, σ 0 , for a Hertzian contact (which occurs at the centre of the contact) can be expressed as a function of the radius of contact, R, the normal force, P, and the overlap of the spherical cap with the elastic half plane, δ, as: Alternatively, σ 0 can be expressed as: where: and where G 1 , G 2 and ν 1 , ν 2 are, respectively, the shear moduli and Poisson's ratios of the two contacting particles. If the maximum stress exceeds the maximum compressive strength of the material, σ cmax , the contact can no longer support the load, P max . In this case, the radius of curvature is recalculated such that σ 0 is equal to σ cmax . Rearranging Eq. 2 and setting σ 0 = σ cmax , the new radius R is given by: The tip of the asperity then vanishes and is replaced with a new, more rounded, cap, as shown in Fig. 9b. The rate of recession of the asperity is governed by the cone angle α together with material strength. If the tangent to the cone is parallel to the surface of the cap at the point where the cap and the cone meet, the new offset from the apex of the cone to the highest point on the cap, δ c , may be calculated from geometry as: In implementing the model, damage to the contact surface is calculated incrementally and changes at each new simulation step are assumed to be small. If σ cmax is exceeded at a contact, then the new radius of contact, R, and offset, δ c , are calculated. The new normal force, P, is calculated on the basis of the new overlap between the two particles, δ, (which will have reduced as the surface of the asperity has receded) and the new contact stiffness, k n , which has changed due to the change of contact radius: where: The new shear stiffness, k s is given by: and the new shear force is calculated incrementally as usual. This model is relatively simple and can potentially satisfy the simultaneous requirements of an apparently low initial stiffness in monotonic loading at high confining pressures (through crushing of the asperity), while retaining high elastic stiffness for the case of cyclic loading, which is of crucial importance in modelling railway ballast. Finding the values of minimum contact radius, R min , and cone angle, α, which are assumed to be the same for all contacts, is addressed in Sect. 2.1.

Calibration of the conical damage model
The parameter values for R min and α were refined using the following iterative procedure: 1. Select the elastic shear modulus, G. (Poisson's ratio, ν, was set to 0.2 for all simulations).

Fig. 10
Percentage of contacts that do not crush as a function of contact radius, R, for shear stiffness, for G = 5 GPa, and a confining pressure of 15 kPa. As a reference, the D 50 of the specimen is 13.5 mm 2. Select R min : it was assumed that very little crushing occurs at a confining pressure of 15 kPa, thus the behaviour is dominated by the influence of R min on the stiffness of the contacts.
(a) Using the contact force data at peak strength for the best fit simulation at a cell pressure of 15 kPa (Fig. 2), from Eq. 4 calculate the value of the contact radius required to support the load, based on the elastic modulus chosen in (1). The proportion of non-crushing contacts as a function of contact radius, for a shear modulus of 5 GPa, is shown in Fig. 10. (b) Choose a value of R min so that most of the contacts will not crush.
3. Select α: with all other parameters held constant, the cone angle, α, determines the rate of plastic deformation normal to the contact plane in response to a normal force increment beyond the elastic limit. An initial estimate for alpha can be obtained by studying the simulation data for the best fit at a cell pressure of 200 kPa with purely elastic behaviour normal to the contact plane (Fig. 7). In this case, the fit to the laboratory data was obtained by using a low inter-particle stiffness, resulting in relatively large particle overlaps. By reapportioning these particle overlaps between plastic and elastic deformations, a more realistic material stiffness can be used, and a value of α may be determined as follows: (a) At the peak of mobilized friction, examine the contact forces, P, and particle overlaps, δ. (b) Calculate the radius of curvature, R, required to support the force at each of the (non-zero force) contacts, using Eq. 4. (c) From Eq. 6, calculate the elastic deformation, δ e , that would result from the application of this force, using the radius, R, and the shear modulus, G. (d) Calculate the plastic displacement, δ p , that would give rise to the same overall contact displacement, δ, as δ p = δ − δ e . This minimum value of δ p is limited to zero for the purpose of the next calculation. (e) Use geometry to calculate the angle α for each contact such that the new contact radius, R, is obtained for a plastic reduction in height of the asperity, δ p , as: A value of alpha can be calculated for each contact. Next, select a single value that satisfies the largest number of contacts. For example, in the case of G = 5 GPa and R min = 4 mm, a histogram of values of α for all particle contacts is shown in Fig. 11, which shows a peak at α ≈ 78 • . 4. Run a cyclic loading test and examine the resulting resilient modulus. Observation of the results for different parameter values shows that: (a) The relationship between shear modulus and resilient modulus is rather non linear. (b) The resilient modulus is not very sensitive to α or the interparticle friction angle, φ. (c) The plastic axial strain is affected by α and φ.
Adjust the shear modulus, by iteration, to obtain the correct residual modulus and then repeat steps (1)-(4) to update R min , α and G.
Finally, these parameters can be fine-tuned through simulation of cyclic and monotonic triaxial tests. For σ cmax = 200 MPa, G = 5 GPa, R min = 4 mm and α = 78 • , a close match to the residual modulus was obtained and the results of the monotonic triaxial test simulations are shown in Fig. 12. At 15 kPa, the results are relatively unaffected by the damage At 200 kPa, the initial stiffness of the response matches the data very closely, but the peak strength is much higher than for the real ballast. This is not all that surprising, as the model was conceived to address the discrepancy in monotonic and cyclic stiffness noted in Sect. 1. A variable friction model is introduced in the next section to address the reduction in peak strength observed at higher cell pressures.

Variable friction model
In Sect. 1 it was shown that at a confining pressure of 200 kPa, a good fit to the modelling data could be obtained through a reduction of both the particle contact stiffness (to an unrealistically low value) and the inter-particle friction angle. The stiffness anomaly was addressed through the introduction of a damage model in Sect. 2.1, but the peak strength remained too high. Reduction of the inter-particle friction angle would improve the response at higher confining pressures, but this is not a viable approach for more general loading conditions. This section describes the rationale for and implementation of a per-contact variable friction model to address the reduction in peak effective strength with confining pressure.
As an alternative to setting the inter-particle friction angle globally as a function of the confining pressure, determination on a per-contact basis could provide a more generally applicable model. Leaving aside for the moment the justification for such a measure, the interparticle friction angle would need to be varied with some aspect of the current state of the contact, for example the normal force. Figure 13 shows the distributions of normal force at peak mobilized effective  Owing to the greater number of contacts per particle at 200 kPa than at 15 kPa (about 60 % more) and to the slightly lower mean pressure at peak strength, the ratio of mean normal forces between the two simulations is about 7. A relationship between interparticle friction angle and contact force would need to devised in such a way that most of the contacts at 15 kPa have a value of around 40 • , reducing to about 30 • for the 200 kPa case.
A possible justification for a reduction in strength at higher confining pressures is that asperities can provide some degree of interlocking between particles at low confining pressures, but tend to break at higher pressures. Evidence for this type of polishing is presented by Indraratna and Salim [14] for ballast subjected to cyclic loading. This process is imagined in the spirit of the well-known saw-blade analogy for dilation in Fig. 14, focusing on a change in apparent friction angle rather than dilation. In Fig. 14a, there is very little wear and the contact normal is at an angle β to the vertical. The apparent friction angle is then φ s +β, where φ s is the friction angle of the smooth material. As the asperities become more worn down, as represented by an increasing radius of curvature in Fig. 14b-d, β tends to zero as the contact normal nears vertical. Thus the apparent friction angle approaches the smooth limit, φ s . The concept proposed by Pöschel and Herrmann [24] in the context of earthquake modelling, in which friction reduces with increasing sliding velocity, offers a possible precedent.

Variable friction model: a first attempt
Consider a pair of contacting particles with a single breakable asperity as shown in Fig. 15a. In this idealized case, before sliding can occur the magnitude of the shear force must exceed the sum of the asperity strength and the frictional limit, which is proportional to the magnitude of the normal force. This limiting shear force is indicated by the solid line in Fig. 15b. If this limit is exceeded, the asperity breaks and the limiting shear force would then be purely frictional, as indicated by the dotted line in Fig. 15b. Trial simulations with this model demonstrated that a very modest asperity strength of about 1 Newton combined with a residual friction angle of about 28 • provided a close fit to the data at confining pressures of both 15 and 200 kPa. At 200 kPa, the asperity strength has very little effect. However, between these values, it was found that the triaxial strength fell away too quickly with increased confining pressure. To address this, a power law relationship is proposed for the shear force limit of the form: where Q lim is the shear force limit, φ residual is the residual friction angle and a and b are constants. The coefficient of friction, μ, can then be written as: Fig. 15 a Two particles in contact; the lower particle has a nominally flat upper surface but includes an asperity that must break before the particles can slide. Arrows indicate the relative forces on the particles. b Shear force limit as a function of normal force magnitude for single asperity case The results of monotonic triaxial test simulations using the conical damage model in combination with the powerlaw friction model are shown in Figs. 16, 17, 18 and 19. The parameters for the simulation are summarized in Table 4 and were obtained as follows: 1. The residual friction angle, φ residual was set to 28 • . The coefficients α and R min were set to the values obtained in Sect. 2.1. 2. The magnitude of the normal force above which the friction coefficient is constant was set to 80 N (see Fig. 16). This results in most of the contacts at a confining pressure of 15 kPa having an increased friction coefficient, and many of the higher-load carrying contacts at 200 kPa the residual value (see Fig. 13). 3. A trial value for b was chosen. The parameter a was calculated from Eq. 11 with P = 80 N (from step 2) and  In all cases, the agreement in the mobilized friction response is excellent and in the volumetric response at least acceptable. Apart from the improved fit to the data, one piece of circumstantial evidence of increased effective interparticle friction at low normal forces is that the mobilized friction response is more similar in nature to the laboratory data with smaller fluctuations (for example, compare Fig. 17 with the responses at a confining pressure of 15 kPa in Figs. 3 or 12). This suggests that the specimen is stabilized (exhibiting less erratic stick-slip behaviour) by lightly-loaded contacts which have increased tangential strength due to surface interlocking between particles.

A damage-dependent friction model
The friction model proposed in Sect. 3.1 relates the interparticle friction angle to normal force at the contact. Although the simulation results offer a good fit to the data, there are two deficiencies with the model: 1. The normal force should be expressed in a non-dimensional form so that the model can be applied at different scales. 2. There is no link between the state of damage and the inter-particle friction angle.
As a first step in addressing these issues, consider the model for the inception of sliding of a spherical contact on a plane, proposed by Kogut and Etsion [16]. In this model, rather than assuming a friction angle at which sliding will occur, the area in contact is assumed to be fused and sliding only takes place if the material itself fails. This leads to an increase in the apparent friction coefficient at lower normal force magnitudes. Based on a fit to results from a finite element model of a sphere contacting a plane, Kogut and Etsion propose a relationship for the friction coefficient, μ, as a function of the normal force, P, normalized by the critical force, P c , i.e. the value of normal force at which plastic deformation will occur for a given sphere radius and set of elastic properties: valid in the range 0 ≤ P ≤ P c , where: where R is the radius of the sphere, K is a hardness coefficient related to the Poisson's ratio of the material by K = 0.454 + 0.41ν and H = 2.8Y , according to [27], where Y is the yield strength of the material. Note that Eq. 13 is very similar to Eq. 2, which can be rearranged to give a slightly different expression for the critical force, P c , as: The difference is due to the fact that the maximum shear stress occurs below the surface for a spherical indenter acting on a plane [6]. Also, note that the form of Eq. 12 is similar to the power law relationship in 3.1 (Eq. 11). The main differences are the normalization of the normal force in Eq. 12, and the lack of a lower limit on the friction coefficient, μ.
Kogut and Etsion were studying the behaviour of metals and their approach may not be directly applicable to rock contacts. However, several authors have reported increased friction with rough surfaces (see for example [7,8,23]). In the case of railway ballast, it is possible to imagine that surface roughness, on a smaller scale than the asperities considered in this paper, could provide a resistance to lateral movements through interlocking. At low normal forces, the interlocking effect would dominate, becoming less significant with increasing normal force. However, the lateral strength of the interlocked region depends on the contact area, which in the case of nominally Hertzian contacts is a function of the effective radius of contact. For constant surface roughness and normal load, an increase in the nominal radius of contact will increase the area of contact and could generate a greater degree of interlocking. This provides a possible link between the conical damage model (see Sect. 2), which alters the effective contact radius, and the frictional behaviour.
Building on Kogut and Etsion, the conical damage model can be combined with a friction law of the form given by Eq. 12, by using the radius of contact from the conical damage model (Eq. 4) to calculate the critical force, P c . The proposed model for μ is given by: where β > 0, γ < 0 and P c is defined by Eq. 14, using the value for R calculated in Eq. 4. It is important to note that the normal force, P can never exceed the critical force P c , as this would result in damage to the cone, an increase in the contact radius and, consequently, in the critical force, P c . As a result, the friction coefficient, μ, has a lower limit given by μ = β.

Model calibration
The parameters β and γ were derived from the values of residual friction, φ residual , and the exponent, 1/b − 1 (see Eq. 11) determined in Sect. 3.1. Thus, β = tan(28 • ) and γ = 1/1.2 − 1 = −0.16. The cone angle, α, and shear modulus, G were initially set to the values obtained in Sect. 2.1. By studying the distribution of contact normal force magnitudes shown in Fig. 13, and making the assumption that most of the contacts do not crush at a confining pressure of 15 kPa, the load at which previously unloaded contacts begin to crush, P 0 c , was initialized to 80 N (from step 2 in Sect. 3.1). The minimum radius, R min was then calculated by rearranging Eq. 14 to give: where E * is given by Eq. 3. Note that the rate of plastic deformation normal to the contact plane for a monotonically increasing load, P max , is determined by both the cone angle, α, and the stiffness, E * . Thus, to maintain similar plastic behaviour when changing the stiffness, α must also be changed. This can be seen by substituting for the cap radius, R, from Eq. 4 into 5 to give a relationship between the normal force and the plastic settlement, δ c , as: is a constant.
Cyclic and monotonic loading simulations were then carried out. The stiffness, G, was then adjusted to obtain a good fit to the resilient modulus measured from the laboratory cyclic loading data. At each stiffness, R min was recalculated using Eq. 16 and α was adjusted to maintain a constant value of A in Eq. 18. Finally, small trial-and-error adjustments were made to R min and α to obtain the correct monotonic loading response across the range of confining pressures.
The best fit parameters are presented in Table 5 Fig. 23; these show excellent agreement for both the resilient modulus and the initial axial strain. After several cycles the rate of axial strain in the simulation is larger than in the laboratory test; this needs further investigation.

Summary and conclusions
A discrete element contact law has been proposed to model the behaviour observed in laboratory triaxial tests of scaled railway ballast at a range of confining pressures from 15 to 200 kPa. It was shown not to be possible to match the laboratory behaviour using the standard Hertzian elastic contact law with an elastic-perfectly plastic tangential slip model, presented in [2]. This suggests that some damage must be occurring at particle contacts. A new contact model was proposed, based on the supposition that particles contact at asperities. The asperity contact was modelled by a spherically-capped cone. In an elastic regime, the contact is Hertzian with a radius defined by the cap radius. If the yield stress is exceeded the cap is crushed and forms a lower cap with a larger radius, able to support the load. This model was shown to provide good agreement with the initial loading behaviour in monotonic laboratory tests for the different confining pressures tested and also satisfies the requirement for a stiffer unloading response. However, the damage model did not significantly alter the peak strength, which reduces with confining pressure in the laboratory tests. Therefore a supplementary model was proposed that varies the interparticle friction coefficient as a function of a normalized load, which is the contact normal force divided by the critical force, P c (the load at which plastic failure will start). For contact normal loads less than P c , the friction is increased; an effect attributed to interlocking at the surface roughness scale. For virgin monotonic loading of a contact, the friction will decrease until the load is equal to P c and the friction coefficient reaches a constant minimum value. A continued increase in loading will result in plastic damage and a consequent increase in P c to the current load value. Subsequent unloading results in an increase in friction coefficient.
Simulation results show excellent agreement with the laboratory data for all monotonic triaxial tests. Good agreement was obtained with cyclic loading data over the first few cycles, in terms of both the resilient response and the axial strain. In later cycles, the simulation displays larger plastic axial strain than the laboratory results and this remains a subject for further investigation.