The effect of consolidation on undrained behaviour of granular materials: experiment and DEM simulation

B

Skempton’s pore pressure parameter

C

number of particle contacts

C_N

coordination number

D_I

densification index

D_r

relative density

d

diameter of particle

d

incremental rate

e

void ratio

e_low

void ratio on lower consolidation line

e_max

maximum void ratio

e_min

minimum void ratio

e_up

void ratio on upper consolidation line

F

fabric tensor

F_ij

scalar value of fabric tensor

F_vM

von Mises (deviatoric) fabric

I

inertial number

J₁, J₂ and J₃

invariants of sF

K

principal stress ratio (⁠$σ3′/σ1′$⁠)

K₀

principal stress ratio ‘at rest’

M

critical state stress ratio (q/p′)

N

number of particles

N₀

number of particles having no contact

N₁

number of particles having only one contact

n^k

direction (unit vector) of the kth contact

p′

mean effective confining stress

q

deviatoric stress

sF

deviatoric fabric tensor

sF_ij

scalar value of fabric tensor

u_F

flow potential

V

volume

Z_m

mechanical coordination number

Δ

rate of change

Δu

pore water pressure

δ_ij

Kronecker delta

$ϵ˙$

shear strain rate

ϵ_q

deviatoric strain

ϵ_v

volumetric strain

ϵ₁ = ϵ₁₁

axial strain

ϵ₂ = ϵ₂₂

lateral strain

ϵ₃ = ϵ₃₃

lateral strain

η

stress ratio

θ_m

mobilised friction angle

ρ_g

particle density

σ

stress tensor

σ₁

major total stress

σ₂

minor total stress

σ₃

minor total stress

$σ1′=σ11′$

major effective stress

$σ2′=σ22′$

minor effective stress

$σ3′=σ33′$

minor effective stress

χ₁

average force imbalance divided by the average contact force

χ₂

average moment imbalance by the average contact force and particle radius

ψ

state parameter

0

start of shearing

CS

feature at critical state

IS

feature at instability

peak

feature at peak deviatoric stress

PT

feature at phase transformation state

The undrained behaviours of granular materials are frequently classified as non-flow (NF), flow (F) and limited flow (LF) in the literature (Alarcon-Guzman et al., 1988; Bobei and Lo, 2005; Mohamad and Dobry, 1986); schematic diagrams of these behaviours are shown in Figures 1(a) and 1(b). Loose specimens exhibit F behaviour where the deviatoric stress (q) quickly attains an initial peak (q_peak), represented by the yellow circles in Figures 1(a) and 1(b), and then q softens until reaching an equilibrium at the critical state (CS); where $q=(σ1′−σ3′)$⁠, and $σ1′$ and $σ3′$ are the principal effective stresses in the major and minor directions. For the F behaviour, the pore water pressure (Δu) increases throughout shearing and reaches a peak Δu, where dΔu = 0 at CS (see Figure 1(c)). When Δu becomes equal to the value of initial p′, then it is referred to as complete liquefaction; otherwise, it is simply referred to as liquefaction hereafter, where p′ is the mean effective stress and equal to $(σ1′+2σ3′)/3$⁠. Medium dense specimens exhibit LF behaviour, where q hardens after initial softening. The strain softening as a form of dq/dϵ_q < 0 in both the F and LF behaviours is the triggering of static liquefaction, where ϵ_q is the deviatoric strain. The triggering of static liquefaction can be presented by the instability stress ratio (η_IS = q/p′) in the q–p′ space, as shown in Figure 1(b). Dense specimens exhibit NF behaviour, where q hardens until reaching the CS. In the case of LF and NF behaviours, Δu initially increases (contractive) to a local peak and then decreases (dilative). The transition between contractive and dilative states – that is, the phase transformation (PT) state in LF and NF behaviours – occurs at the ‘knee’ of the effective stress path (ESP) in q–p′ space or the local peak Δu in Δu–ϵ_q space, represented by the yellow diamonds in Figure 1 (Ishihara et al., 1975). Therefore, both η_IS and PT are important characteristic features of undrained behaviour. The anchor concept of the critical state soil mechanics (CSSM) framework is that the CS data points from these tests form a line in the e–log(p′) space, which is commonly known as critical state line (CSL). The CSL is the reference line to define the state of a soil. The characteristic features – for example, η_IS and PT – are assumed to be dependent on soil states with respect to CSL. Therefore, simple functional forms are suggested for such states with respect to CSL that can be measured from a laboratory triaxial test or estimated from a field test (Chu et al., 2015; Shuttle and Jefferies, 2016) and also incorporated into constitutive modelling (Imam et al., 2002; Ishihara, 1993; Jefferies et al., 2015; Li and Dafalias, 2000).

The CSSM framework is elegant and more widely evaluated with experimental data for granular materials after isotropic consolidation than K₀ consolidation. However, there are experimental evidences that the consolidation history may affect the undrained response of sand and its corresponding CS behaviour. For instance, Finno and Rechenmacher (2003) reported that the location of CSL is affected by the consolidation history and initial state (i.e. before consolidation) of sand. The study of Fourie and Tshabalala (2005) found different CSLs for isotropically and K₀-consolidated specimens. They also found that η_IS values of K₀ consolidation were higher than those of isotropic consolidation. However, many other studies found a unique CSL for a soil regardless of consolidation condition (Kato et al., 2001). The uniqueness or non-uniqueness of CSL for different consolidations is a significant consideration for establishing the CSSM framework and its application for engineering practice. Due to the disagreement in literature, the effect of consolidation on the CS behaviour of soil requires further investigation.

Recognising the limitations of element tests under laboratory conditions, such as the sensitivity of monitoring transducers, disturbance/non-uniformity of specimens and shear band formation, the discrete-element method (DEM) has been gaining popularity for the qualitative understanding of soil behaviour under the CSSM framework. DEM does not have the same limitations as laboratory experimentation and also allows the understanding of particle interactions. There are published studies on CS behaviour of the assembly of two-dimensional circular and three-dimensional (3D) spherical particles (Huang et al., 2014; Kuhn, 2016; Sitharam and Vinod, 2009; Yan and Dong, 2011; Zhao and Guo, 2013); however, there is no detailed study on the effect of isotropic and K₀ consolidations on undrained response, its characteristic features and CS behaviour.

Therefore, the objectives of this study are to (a) evaluate the effect of isotropic and K₀ consolidations on undrained behaviour of saturated yellow silty (YS) sand and DEM simulations of 3D cubical assembly of ellipsoid particles under triaxial condition. This study did not attempt to simulate undrained behaviour of YS sand; rather, it attempted to evaluate undrained characteristic features – for example, η_IS, PT and CS – of granular materials that were observed for both laboratory triaxial test and DEM simulation. It is envisaged that such comparative evaluation between laboratory test and DEM will lead to the better understanding of the fabric-dependent behaviour of granular materials – for example, Yoshimine et al. (1998) and Li and Dafalias (2012); (b) evaluate unique/non-unique CSLs for isotropic and K₀ consolidations; (c) evaluate the correlation for the undrained characteristic features – for example, PT state and η_IS – with state indices; and (d) evaluate the correlation for micromechanics measures such as particle contact and fabric anisotropy with state indices. This study enhances the understanding of the fundamental basis of the relationship between macro- or micromechanical responses and state indices of cohesionless granular materials.

This study is composed of the triaxial testing of YS sand and the discrete element simulation of 3D ellipsoid particles. The physical properties of YS sand and the numerical properties of ellipsoid particles as well as laboratory testing and numerical simulation methodologies are explained below.

The YS sand, from South Australia, was used in this study. It is poorly graded and contained about 10% fines (poorly graded silty sand (SP-SM) according to the Unified Soil Classification System (USCS)), where fines are particles that have sizes smaller than 0·075 mm. The minimum and maximum void ratios are 0·423 (1·86 g/cm³) and 0·892 (1·40 g/cm³), respectively, according to ASTM D 4253 (ASTM, 2016) and ASTM D 698 (ASTM, 2017) methods, respectively. The clean sand is uniform and subrounded in nature. The fines were obtained by wet sieving of YS sand and classified as low plastic silt (ML) according to USCS which have a liquid limit (LL) and plasticity index (PI) of 40·0% and 9·3%, respectively.

The reason for using 3D elliptical particles is to capture an essence of angularity and reduce the computational demand as required for clustered particles. The grain size distribution curve of YS sand and 3D ellipsoids with their distribution of eccentricity – that is, height to width ratios – are shown in Figure 2. It is noted that the sizes and shapes of individual grains and specimens for laboratory testing and DEM simulation were different. Ellipsoid particles in DEM were larger and smoother than natural YS sand particles. A smaller and cubical specimen was used for DEM simulation, whereas a large cylindrical specimen was used in laboratory triaxial testing. These variations are commonly adopted for DEM simulations for simplicity and less computational demand as reported in many studies (Gu et al., 2014; Huang et al., 2014). As this study did not attempt to simulate YS sand behaviour with DEM, it rather attempted to evaluate the characteristic features that was reported for both laboratory triaxial test and DEM simulations, and the objective of this study was not significantly influenced by these variations. The DEM particles were assigned a normal contact stiffness of 10⁸ N/m, ratio of tangential to normal contact stiffness of 1, coefficient of friction at particle contacts of 0·5 and the coefficient of (rotational and translational) body damping of 0·05.

Different specimen preparation methods have been used in the literature such as air pluviation, water sedimentation and moist tamping for triaxial testing. Each method offers different advantages in resembling certain field conditions and fabrics. The moist tamping method provides two advantages for this study: (a) it produces very loose specimen suitable for liquefaction studies (Yamamuro and Lade, 1997) and (b) it reduces segregation of fines particles (Baki et al., 2012; Ladd, 1978; Rahman, 2009), which was a concern for YS sand. Therefore, a modified moist tamping method was adopted in this study (Bobei et al., 2009; Lo et al., 2010). A moist tamped specimen provides an initial soil fabric more representative of sandy deposits placed on land by end tipping, such as embankment fill or colluvium (Rahman et al., 2014a).

A 100 mm dia. and 100 mm high split mould of was placed over the enlarged bottom platen with free end, and a small vacuum of approximately −10 kPa was applied to hold the membrane inside this mould. The moist sand was tamped in ten layers to increase the homogeneity and reduce the variation in local density between neighbouring layers. The thickness of each layer was controlled using a height-controlled tamping rod. Free ends with enlarged platen reduce end restraint to a minimal value as verified in a number of studies on silt size to sand size materials (Bobei et al., 2009; Lo and Wardani, 2002; Lo et al., 2010). Figure 3(a) shows a specimen just after dismantling the split mould. Note that specimens with height-to-diameter ratio were unity, because it (a) was easier to prepare in terms of controlling uniformity of initial density and (b) minimised the likelihood and extent of specimen tilting during shearing. This approach has been well documented in literatures (Baki et al., 2014; Bobei et al., 2013; Chu et al., 1993; Rahman et al., 2008; Yamamuro and Lade, 1998; Zhang et al., 2018).

Carbon dioxide (CO₂) was pumped into the specimen through the bottom platen for at least 10 min before the dismantlement of the split mould. Then saturation was performed in two stages: vacuum flushing and back pressure saturation. In stage 1, the specimen was flushed under the vacuum of −20 kPa overnight with de-aired water, and then in stage 2, a back pressure of 400 kPa was applied maintaining a p′ of 20 kPa. The specimen was left to rest with the back pressure for approximately 5 h to help dissolve any remaining air in the water. In most cases, a Skempton’s B value of 0·98 or higher was achieved in each test to ensure adequate saturation.

The consolidation was conducted using a computer control module. During isotropic consolidation, the module increased cell pressure to the target cell pressure at specified intervals while allowing pore-water drainage to keep the pore-water pressure (back pressure) constant. While increasing isotropic cell pressure (σ₁ = σ₂ = σ₃), the q was maintained constant to the user-defined value, usually q = 1 kPa, in the bedding control. Although the bedding control was not necessary for isotropic consolidation, it was maintained for consistency with other testing – for example, K₀ consolidation and extension test. The digital pressure volume controller (DPVC) was used to measure the volume change during isotropic consolidation as well as keeping back pressure constant. Independent pressure transducers were used just outside the cell to check pore-water pressure equilibrium.

First, a constant axial strain rate was applied, producing an increase in the effective vertical stress (⁠$σ1′$⁠) up to the target value. For the applied vertical strain, zero radial strain was enforced by the DPVC in the manner described by Menzies (1988). In this module, a feedback system calculated the change of specimen volume (Δv) equal to the axial deformation increment multiplied by the current cross-sectional area of the specimen for a time step (default: 6 s). The time step was adjustable depending on the material response. The feedback system then commanded the DPVC to extract the volume of Δv to maintain zero radial strain – that is, ϵ_v = ϵ₁ as shown in the following.

where ϵ_v is the volumetric strain. The extraction of volume causes a subsequent reduction to the pore pressure u (back pressure). To avoid the potential problem that a lower pore pressure may cause partial desaturation, the program commands the cell pressure regulator to increase the cell pressure by the drop in pore pressure. As change in pore pressure, Δu = BΔ₃, where B is Skempton’s pore pressure parameter, and assuming that the specimen is nearly saturated (B ≈ 1), this action will restore the back pressure to the desired level. If B is somewhat less than 1, then there may be a lag in the restoration of the back pressure.

Once the consolidation is completed, the program has two exit modes, either ‘freeze stress’ or ‘freeze strain’. In the freeze stress mode, the control program keeps cell pressure, deviator stress and back pressure at their final values. Granting that pore pressures may still be equalising to some degree, and that the specimen may creep a small amount, the final K₀ condition may be compromised if the specimen is held in this state for a protracted length of time. On the other hand, in the freeze strain mode, the ram is stopped at the exit of K₀ consolidation – that is, axial strain is fixed, and the DPVC ceases volume change and reverts to pore pressure measurement. The cell pressure and pore pressure remain constant, but the deviator stress may keep changing slightly, depending on specimen stiffness. In this study, the freeze stress exit mode was used for all the tests.

The axial force was applied to the specimens at a constant deformation rate, and the cell pressure was kept constant at the initial value. In the undrained condition, the DPVC did not allow volume change, but the change of pore water pressure was measured.

The DEM software called Oval, developed by Kuhn (2006), was used for this study. The 3D ellipsoids assembly of cubical specimen under triaxial condition was simulated. Two methods were employed for preparing ellipsoid assemblies: staged isotropic consolidation and particle elimination (discussed later). A linear contact model and periodic boundary condition were considered in this study.

A staged isotropic consolidation approach, similar to that of Guo and Zhao (2013), was employed to achieve different densities. A set of specimens with a large number of particles (particle number N varied from 4350 to 5400) with the same grain size distribution was randomly generated in cubical 3D space. It is noted that there was no contact between particles at the beginning of this stage. The particles were then assigned with varying coefficients of friction at particle contacts (0–0·50) and 3D isotropic stress (⁠$σ1′=σ2′=σ3′$⁠) up to 20 kPa was applied to achieve a desired void ratio. A higher coefficient of friction was used to generate a loose assembly. Note that after this stage, the coefficient of friction of all simulations was 0·50.

The staged isotropic consolidation introduces different stress histories for different coefficients of friction at particle contacts during consolidation of up to 20 kPa. To avoid this, a new specimen preparation method was sought, and the particle elimination method was developed. In this method, a reference specimen was firstly achieved with ellipsoid with a very high density at p′ of 20 kPa by the staged isotropic consolidation. A MatLab code was developed to eliminate some ellipsoid particles from a dense specimen in such a way that the grain size distribution and particle orientation remained the same. This is discussed by Nguyen et al. (2017). An example of ellipsoid assembly is shown in Figure 3(b). The smaller specimen size for this study was approximately 9·33 × 9·33 × 9·33 mm³. The equivalence of undrained behaviour of specimens prepared by these two methods was observed and discussed by Nguyen et al. (2017).

A representative elementary volume (REV), which is the smallest volume that can yield a representative value of the whole, is important for numerical simulation. To achieve an REV condition in DEM simulations, micromechanical measures, such as inertial number (I), pseudostatic measures and specimen size, were examined in this study.

The quasi-static conditions of the DEM simulations were checked and maintained by analysing the inertial number (I), which is given by the following equation.

where $ϵ˙$ is shear strain rate (1/s), d is the particle diameter (m), p′ is the mean stress (kPa), and ρ_g is the particle density (kg/m³). The particle density (ρ_g) was 9·766 × 10⁸ kg/m³, which is approximately 3·69 × 10⁵ times larger than 2650 kg/m³ for real geomaterial as reported by Perez et al. (2016). To optimise the simulation time, the particle density was scaled up in the literature, and notably, Thornton (2000) used a scale-up factor of 10¹². A similar scale-up approach was also observed by Ng (2006). The strain rate of 0·00001% per time step was found as suitable for all DEM simulations (Nguyen and Rahman, 2017; Nguyen et al., 2015). The calculated I values for all simulations in this study were much smaller than the recommended value of 0·003 for quasi-static conditions (da Cruz et al., 2005).

Two variables, χ₁ and χ₂, were also checked to determine the numerical stability of the simulation. χ₁ is the average force imbalance on a particle divided by the average magnitude of a contact force. χ₂ is the average moment imbalance on a particle divided by both the average magnitude of a contact force and the average particle radius. For the applied strain rate, the highest values of χ₁ and χ₂ were below the recommended value of 0·005 (Kuhn, 2006) and considered that simulations were quasi-static.

The Oval software allows a flexible stress or strain control triaxial simulation scheme by a six-digit input variable ‘icontr’, which represents the components of three principals in 11, 22 and 33 directions, and three shear 12, 13 and 23 directions. It should be noted that there are nine directions (i.e. three principal and six shear directions). However, the components in the directions 12 and 21, 13 and 31, 23 and 32 have the same magnitude but opposite direction. The different conditions in DEM simulations were achieved by adjusting these six components. The details are presented in the following.

• Isotropic consolidation. The rates of change of principal stresses are the same – that is, $dσ11′=dσ22′=dσ33′$⁠. Note that single digit notations for principal stresses – that is, 1, 2 and 3 – were used in subscript when explaining the triaxial testing.

• K₀ consolidation. For an applied rate of major principal stress, the lateral strains are maintained constant or no lateral deformation – that is, dϵ₂₂ = dϵ₃₃ = 0.

• Undrained shearing. For an applied axial strain rate, the volume of the assembly is maintained constant by adjusting the lateral strains. There was no pore water pressure development (Δu) during DEM simulation. However, it was determined by the difference between p′ for drained and undrained paths as explained by Sitharam et al. (2008) and Nguyen et al. (2018).

DEM allows the examination of particle interactions or contacts to establish links with the macroresponse of granular materials. Therefore, the coordination number (C_N) and mechanical coordination number (Z_m), as defined by Rothenburg and Bathurst (1989) and Thornton (2000), respectively, were used in this study. They can be mathematically expressed as follows.

where C_N is the coordination number, C is the total number of contacts, N is the number of particles, N₀ is the number of particles having no contact and N₁ is the number of particles having only one contact in the specimen. Note that previous studies showed that Z_m > C_N, and their evolution paths are parallel during undrained shearing (Nguyen et al., 2017). Therefore, C_N, instead of Z_m, is used to synthesise the simulation data.

The structural anisotropy can be determined by the fabric tensor F, which was firstly introduced by Satake (1982) and expressed in the following equation.

where n^k is the direction (unit vector) of the kth contact and N is the number of contacts in the specimen. The bold character (e.g. F) indicates the tensor, whereas the italic character (e.g. F_ij) indicates the scalar value of the tensor. F₁₁, F₂₂ and F₃₃ are the fabric in the principal directions (i = j), and F₁₂, F₁₃, F₂₁, F₂₃, F₃₁ and F₃₂ are the fabric in the shear directions (i ≠ j). As mentioned earlier, the scalar values of shear components in opposite direction are equal – that is, F₁₂ = F₂₁, F₁₃ = F₃₁ and F₂₃ = F₃₂. The fabric tensor can be formulated based on the theory of Cauchy stress tensor (σ), and the deviatoric fabric tensor (sF) can be expressed as

where

where the notation ‘tr(…)’ indicates the trace of a tensor and δ_ij is Kronecker delta. sF has a set of invariants – for example, J₁, J₂ and J₃ – and its characteristic equation can be written as

where λ is the constant of proportionality. The invariants of sF can be defined as

Based on von Mises stress yield criterion, a similar scalar quantity of equivalent deviatoric fabric – that is, von Mises fabric (F_vM) – can be defined as

Figure 4(a) shows four pairs of isotropic and K₀ consolidation paths for the same preconsolidation e for both triaxial and DEM simulations in the e–log(p′) space. Although the loose triaxial test pair (YS-CIU05 and YS-K₀U13) had approximately the same preconsolidation e of 0·670, the K₀ consolidation path showed higher compression. Similarly, the dense triaxial test pair (YS-CIU15 and YS-K₀U08) had approximately the same preconsolidation e of 0·590, the K₀ consolidation path showed slightly higher compression. However, the difference of compressibility for isotropic and K₀ consolidations was higher for higher preconsolidation e. Similarly, the loose (DE-CIU24 and DE-K₀U04) and dense (DE-CIU21 and DE-K₀U03) pairs from DEM simulations showed different compressibilities for different consolidation paths. The K₀ consolidation paths showed slightly higher compression than isotropic consolidation, which increases with e. The purpose of Figure 4 is not to compare quantitatively the results from experiments and DEM simulations, but to capture the change in compressibility for isotropic and K₀ consolidations and the trend of increasing compressibility with e for YS sand and other experimental studies (Butterfield and Marchi, 2017; Wong et al., 2017). Figure 4(b) shows the stress ratio (⁠$K=σ3′/σ1′$⁠) evolution paths of the triaxial tests and DEM simulations (presented in Figure 4(a)) in the $K−σ1′$ space. It shows that K = 1 was maintained for isotropic consolidation; however, the K quickly dropped from 1 at the early stage of K₀ consolidation (approximately $σ1′<100 kPa$⁠) and then gradually tends to reach a constant value at a higher stress level (approximately $σ1′=400 kPa$⁠). It seems that the constant K – that is, K₀ value – is dependent on the preconsolidation e of DEM simulation. For instance, the preconsolidation e values of DE-K₀U03 and DE-K₀U04 were 0·540 and 0·634, respectively, and the loose specimen showed an early drop of K at lower $σ1′$ than the dense specimen, but it achieved a higher K₀ than the dense specimen – that is, K₀ increases with e. Similarly, the loose triaxial specimen (YS-K₀U13) showed an early drop of K at lower $σ1′$⁠; however, the test was terminated at a smaller $σ1′$ (to maintain the same p′ of YS-CIU05) and could not confirm that K₀ increases with e. However, in general, it was found for triaxial tests that the value of K₀ increases with preconsolidation e. This observation differs from the widely used Jaky (1944) equation – that is, K₀ = 1 – sin θ_m, which suggests that K₀ only depends on the mobilised angle of internal friction (θ_m).

The micromechanical quantities – that is, C_N and F_vM from the DEM simulations – also depicted the difference in consolidation history as shown in Figure 5. Starting at the same value, the C_N evolved to a much lower value under K₀ condition than the C_N under isotropic condition (Figure 5(a)). Note that higher C_N implies stable state. Since K₀ has lower C_N, it shows higher compressibility in Figure 4(a) than isotropic consolidation. Figure 5(b) shows that F_vM fluctuates around zero during isotropic condition. On the other hand, F_vM quickly increases at a low $σ1′$ for K₀ condition and then reaches equilibrium at a higher $σ1′$⁠. It is noted that micromechanical quantities – for example, C_N and F_vM – tend to reach equilibrium at a higher $σ1′$⁠, which is very similar to macromechanical behaviour – that is, $K=σ3′/σ1′$ evolution in the $K−σ1′$ space in Figure 4(b). In other words, the micromechanical measures seem to be able to reflect the macromechanical quantities.

A series of undrained triaxial compression tests for 21 isotropic and 14 K₀ consolidated specimens was conducted. The range of e₀ values was 0·497–0·655 and $p0′$ varied from 50 to 500 kPa, where subscript ‘0’ denotes the start of shearing. The details of these tests are presented in Table 1. Another series of DEM simulations under undrained triaxial condition for 54 isotropic and 15 K₀ consolidated specimens was conducted (Table 2). The range of e₀ values was 0·395–0·729 and $p0′$ varied from 19 to 1000 kPa. Some of these undrained simulations were part of another project and previously published by Nguyen et al. (2017).

The F behaviours in q–p′ and q–ϵ _q spaces after isotropic and K ₀ consolidations of loose triaxial specimen and DEM simulation under triaxial conditions are presented in Figures 6(a) and 6(b) and Figures 6(c) and 6(d), respectively. The ESPs of two isotropic triaxial specimens (YS-CIU04 and YS-CIU05) started from q = 0 and raised to a peak q (q _peak). The strain softening or static liquefaction – that is, dq/dϵ _q < 0 – was triggered just after q _peak. The η _IS of YS-CIU04 and YS-CIU05 were 0·742 and 0·762, respectively, as shown by dotted lines in Figure 6(a). However, after K ₀ consolidation, YS-K₀U10 and YS-K₀U11 reached the initial q of 135 and 193 kPa, respectively. During undrained shearing, the q of these two triaxial tests quickly rose to a q _peak, and then strain softened to static liquefaction. The η _IS of YS-K₀U10 and YS-K₀U11 were 0·726 and 0·796, respectively. There was no obvious correlation between η _IS and consolidation type for these tests. This may be due to the fact that these tests had very different initial q, e and p′. The effect of these variables is discussed later. The stress–strain behaviour of these triaxial tests are shown in the q–ϵ _q space in Figure 6(b). Interestingly, q _peak occurred quickly for K ₀-consolidated specimens (approximately 0·5% strain), whereas the q _peak of isotropic consolidated specimen occurs at around 2–5% strain, as shown in Figure 6(b).

Similar to triaxial tests, the DEM simulation of loose assembly under triaxial condition showed q _peak for both isotropic and K ₀ consolidations as presented in Figure 6(c). However, η _IS of DE-K₀U05 and DE-K₀U11 were 0·805 and 0·849, respectively, which were higher than η _IS of 0·240 and 0·418 for DE-CIU13 and DE-CIU14, respectively. Note that η _IS values for K ₀ consolidated were significantly higher that isotropic consolidation. η, after K ₀ consolidated corresponding to initial q at the start of undrained shearing, was somewhat close to the CS stress ratio M. The stress–strain behaviour of these DEM simulations showed a quick drop from q _peak to q = 0 at CS as shown in Figure 6(d).

The NF behaviour in q–p′ and q–ϵ _q spaces after isotropic and K ₀ consolidations of dense triaxial specimen and DEM simulation under triaxial conditions are presented in Figures 7(a) and 7(b) and Figures 7(c) and 7(d), respectively. Usually, LF and NF responses exhibit the PT state at the knee of the ESP, which can be defined as the state of dp′ = 0 on the ESP in the q–p′ space (Lade and Ibsen, 1997). However, other researchers argued that the PT state separates the contractive to dilative tendency at dϵ _v = 0; therefore, the equivalence of this for the undrained test would be dΔu = 0. Hence, Sukumaran et al. (1996) define PT state as a state of dΔu = 0. Although the triaxial test and DEM simulations in Figure 7 did not show the clear knee of ESP, the latter definition allowed the identification of PT states at dΔu = 0. The stress ratio at PT (η _PT) on the ESP is considered as an important characteristic feature and has been used in constitutive formulation (Manzari and Dafalias, 1997; Wood et al., 1994). η _PT in triaxial tests occurred below the M-line, and after that, ESPs follow the M-line as shown in Figure 7(a). It appears that both K ₀ and isotropic consolidation tests tend to a unique CSL in the q–p′ space – that is, the M-line. Figure 7(c) shows that ESPs of DEM simulations evolved beyond the CSL (M-line) and then come back to CS at large strain (ϵ _q > 15%). η _PT was also higher than M. The stress–strain paths for both triaxial and DEM simulations showed very similar patterns of strain hardening as shown in Figures 7(b) and 7(d). Note that the Δu development in Figures 6 and 7 exhibits a similar pattern to what was observed in Figure 1(c), and Δu data for these tests are shown in Figures 1 and 2 in the online supplementary material.

To evaluate the effect of density index, the relative density (D _r) as defined in the following equation, was primarily used for triaxial specimen of YS sand.

where e _max and e _min are the maximum and minimum void ratios, respectively. However, the measurements of e _max and e _min are challenging and problematic in DEM; therefore, an alternative measure of D _r, called densification index (D _I), was adopted to capture the density of the assembly of ellipsoids. D _I was defined as the state of soil relative to the upper and lower boundary consolidation lines (see Figure 8) as proposed by Zhang et al. (2018). D _I can be expressed as follows.

where e _up and e _low are the corresponding void ratios on upper and lower consolidation lines, respectively. The upper consolidation line (CL_up) was the consolidation line of loosest assembly, whereas the consolidation line of densest assembly was the lower consolidation line (CL_low). So, the density indices for the triaxial tests and DEM simulation were D _r and D _I, respectively, in this study. Note that D _r and D _I are slightly different as D _I is dependent on $p0′$⁠, and therefore, its value for the same e can change with consolidation pressure p′. However, if the isotropic consolidation lines (ICLs) are parallel, the value of D _I does not change with p′.

The relationships between η_IS and the density indices – that is, D_r and D_I – are shown in Figure 9 for both triaxial and DEM simulations. In general, η_IS was increased with D_r for both K₀- and isotropically consolidated specimen in triaxial tests as shown in Figure 9(a). The root mean square deviation (RMSD) of the best-fit trend was 0·100. The D_r values of K₀ tests were comparatively higher than that of isotropic tests, as e after K₀ consolidation were mostly lower than e after isotropic consolidation (as shown in Figure 4(a)). Similarly, η_IS was increased with D_I0 in the η_IS–D_I0 space for both K₀ and isotropically consolidated specimen in DEM simulation as shown in Figure 9(b), where D_I0 is the densification index at the start of shearing. The RMSD of the best-fit trend was 0·268. The η_IS of K₀ simulations seem to be higher than that of isotropic simulation. The C_N at instability (C_NIS) increased with D_I0, but it did not exhibit any dependency on the consolidation type (the RMSD of this relation was 0·416) as shown in Figure 9(c). The F_vM at instability (F_vM,IS) also increased with D_I0, and F_vM,IS of K₀ simulation seems to be higher than that of the isotropic simulation (the RMSD of this relation was 0·021) as shown in Figure 9(d). The pattern of the F_vM,IS–D_I0 relation in Figure 9(d) is somewhat similar to the η_IS–D_I0 relation in Figure 9(b).

Figure 10(a) shows that the η_PT decreases with D_r in the η_PT–D_r space for triaxial tests. The best-fit trend for η_PT and D_r was higher for K₀ consolidation than that of for isotropic consolidation. These two trends appear to meet at D_r = 0·90. However, the η_PT increased gradually with D_I0 in the η_PT–D_I0 space for DEM simulation as shown in Figure 10(b). Interestingly, the best-fit trend did not have any obvious dependency on consolidation type for the DEM simulations of this study, which is different from what was observed for triaxial data in Figure 10(a).

The flow potential (u_F), as defined in Equation 12, is often used to characterise granular material behaviour.

where $pPT′$ is the mean effective confining stress at the PT state. Figure 11(a) shows that the u_F decreases with D_r in the u_F–D_r space for triaxial tests. The best-fit line has an RMSD of 0·361. Note that the u_F values can be negative for $pPT′$ higher than $p0′$ in Equation 12 – that is, the PT occurs at a higher stress level (p′) than the initial state of $p0′$⁠. Unlike the η_PT–D_I0 relation in Figure 10(b), u_F was decreased with D_I0 in the u_F–D_I0 space, and two distinct relations were observed for isotropic and K₀ consolidations for DEM simulation (see Figure 11(b)). Noticeably, most u_F values of K₀ simulations are relatively small, and some tend to negative.

The above investigation suggests that the overall qualitative behaviour of YS sand in triaxial testing can be captured by DEM simulation for triaxial testing with a simple linear contact model for ellipsoid particles. However, there are disagreements between the experiment and DEM data for some characteristic features, which require further investigation and development. Since D_r and D_I were not equivalent density indices, this study investigated the correlation of these characteristic features within the CSSM framework.

This study revisited the classical CSSM framework and investigated the CS behaviour of the YS sands and DEM ellipsoids.

Figure 12(a) shows a pair of ESPs for DEM simulations with almost the same e ₀ of 0·592–0·599 and $p0′$ of 500 kPa after K ₀ and isotropic consolidations. As previously shown in Figure 7(c), the ESPs for DE-CIU71 and DE-K₀U04 evolved beyond the CSL (M-line) and then both the ESPs came back/approached towards the same CS (η = M) at the end of simulation (see Figure 12(a)). Note that the stress state of DE-K₀U04 after K ₀ consolidation or at the start of undrained shearing was higher than M as indicated by an arrow in Figure 12(a). The F _vM paths for these simulations evolved to a narrow zone (see Figure 12(b)). This zone was termed as ‘CS zone’, which was in line with the concept of CS fabric as discussed by Dafalias (2016). Interestingly, both the ESP and F _vM evolution paths were similar – that is, they go beyond the CS and come back to CS or a narrow zone – which suggests that F _vM may be a good micromechanical parameter for characterising the macromechanical behaviour.

Most triaxial tests and DEM simulations in this study reached clear CS. However, in few cases in which a clear CS was not reached, an extrapolation method, as discussed by Rahman and Lo (2014), was used. In this method, the rate of change of excess pore pressure (dΔu/dϵ _q) is plotted against η as shown in Figure 13 for both experimental and DEM data. The value of η at CS (M) can be estimated by extrapolating dΔu/dϵ _q to 0. Then, p′ at the CS is estimated by extrapolating η to M in an η–p′ plot. Such an extrapolation approach was also employed in other studies presented in the literature (Carrera et al., 2011; Murthy et al., 2007; Zhang et al., 2018). In the case of dense specimen in DEM simulation, dΔu/dϵ _q = 0 was achieved twice: at the PT state and at the CS as shown in Figure 13(b). However, for the loose specimen in triaxial test, dΔu/dϵ _q becomes 0 at the CS (Figure 13(a)) – that is, PT and CS may coincide.

The CS data points from triaxial experiments formed two different CSLs in the e–log(p′) space for K ₀ and isotropically consolidated specimens (Figure 14(a)). The finding of non-unique CSLs for YS sand for different consolidation conditions was somewhat unexpected and undesirable for the CS framework. This finding did not coincide with the authors’ earlier findings of unique CSL under triaxial compression, irrespective of drainage conditions, for Sydney sand (Rahman et al., 2014a), pond ash (Zhang et al., 2018) and Hostun sand (Goudarzy et al., 2017). Therefore, it has also been verified for different drainage conditions and confirmed. The details of these verifications are published in a separate article (Rabbi et al., 2018). Note that the non-unique CSLs for different consolidations is not new. Finno and Rechenmacher (2003) observed non-unique CSLs for undrained loading after anisotropic consolidations. Fourie and Tshabalala (2005) observed slightly different CSLs for undrained loading after isotropic and anisotropic consolidations (see Figure 15). However, the CS data points from DEM simulations of this study formed single CSL regardless of consolidation condition. Other DEM studies – for example, that of Zhou et al. (2017) – observed different critical state stress ratios (q/p′) for the drained condition after isotropic and anisotropic consolidations. However, they have not reported CSL in the e–log(p′) space. This suggested that the assumption of a presumptive unique CSL from triaxial experiment or DEM simulation is not necessarily true for all granular materials. Further research is required to understand the material properties that might lead to such deviation from unique CSL.

Figure 16(a) shows that the C _N at CS forms unique relationships with p′. Figure 16(b) shows that F _vM at CS in the F _vM–p′ space also forms a unique narrow CS zone for both K ₀ and isotropic simulations. Although it was previously reported that the fabric evolved during undrained shearing (see Figure 12(b)), Figure 16(b) suggested that all F _vM values at CS formed a unique trend with the corresponding p′. Hence, the fabric anisotropy, which represents the internal structure of the specimen, evolved from an initial value towards a CS value. This supported the fabric evolution and CS fabric concept of the anisotropic critical state theory (ACST) by Li and Dafalias (2012), and Dafalias (2016). The evolution of F _vM in this study was in fact a presentation of the fabric evolution in ACST. Although a fabric evolution model is already suggested for ACST, this study presents evidence for the evolution of fabric (see Figures 12(a) and 12(b)) and its relation with characteristic features to adopt necessary corrections on such model.

The CSL in the e–log(p′) space has been used as the reference state for soil behaviour. For example, a soil state above the CSL exhibits F behaviour, whereas a soil state below the CSL exhibits NF behaviour. A soil state with respect to CSL can be presented by a state index. One of the widely used state indices is the state parameter (ψ). Been and Jefferies (1985) defined ψ as the difference between the current e and the corresponding e on CSL at the same p′. A positive ψ – that is, e above CSL – exhibits F behaviour, whereas a negative ψ – that is, e below CSL – exhibits NF behaviour. The relation between the state parameter (ψ) and characteristic features on undrained behaviour are explored in the following subsection.

The η_IS and ψ₀ data points for triaxial testing after K₀ and isotropic consolidations showed a good relation, where η_IS decreased with increasing ψ₀ (see Figure 17(a)). The overall trend is consistent with earlier studies for isotropic consolidated sands (Baki et al., 2012; Imam et al., 2002; Rahman and Lo, 2012; Rahman et al., 2014b; Yang, 2002). The best-fit line in this study has an RMSD of 0·100. However, one may notice that the η_IS data points for K₀ consolidation were slightly higher than those for isotropic consolidation for higher ψ₀. This behaviour was clearly observed for the η_IS–ψ₀ relation in DEM simulation – that is, η_IS for K₀ consolidation was significantly higher than that for isotropic consolidation for higher ψ₀ (RMSDs for the relations after isotropic and K₀ consolidations are 0·313 and 0·084, respectively), as shown in Figure 17(b). The unique relationship between C_N at η_IS and ψ₀ with RMSD of 0·364 is presented in Figure 17(c), which suggests that the difference between η_IS–ψ₀ relations for K₀ and isotropic consolidations are not dependent on C_N. The F_vM and ψ₀ data points are presented in Figure 17(d), which shows that F_vM–ψ₀ for K₀ consolidation was significantly higher than that for isotropic consolidation for higher ψ₀ (RMSDs for the relations after isotropic and K₀ consolidations are 0·021 and 0·010, respectively). This is similar to macromechanical behaviour – that is, the η_IS–ψ₀ relation in Figure 17(b). This again suggests that the fabric quantity (F_vM) may represent the macromechanical behaviour better than the number of contacts (C_N).

Figure 18(a) shows that η at the PT state (η_PT) form a unique relationship with ψ₀ in the η_PT–ψ₀ space for both isotropic and K₀ triaxial experiments. The best-fit trend has an RMSD of 0·117. One may recall that these data exhibited separate trends in the η_PT–D_r space in Figure 10(a) for isotropic and K₀-consolidated specimens. It appeared that ψ₀, which captures the effect of e and p′, is a better state index than D_r. Figure 18(b) also shows the unique relationship between η_PT and ψ₀ for both DEM isotropic and K₀ simulations. The best-fit trend has an RMSD of 0·081. However, one may note that η_PT increased with ψ₀ for triaxial data, whereas η_PT decreased with ψ₀ for DEM data. While these are observed behaviours from this study, authors were cautious to draw a general conclusion that triaxial data are more reliable/acceptable than DEM. This is because (a) both triaxial and DEM simulations did not show a clear knee on the ESP at dp′ = 0, and an alternative approach of dΔu = 0 was used to identify the PT state; (b) Zhang et al.’s (2018) triaxial data also show a similar trend as observed for DEM simulation. Although the authors suggest further study on the shape of η_PT–ψ₀ curve, ψ appeared as a better state of soil.

The u_F–ψ₀ relations for triaxial tests after isotropic and K₀ consolidations are presented in Figure 19(a), which shows that u_F increased with ψ₀. These data were not dependent on consolidation type. The best-fit line for these data has an RMSDs of 0·280. The u_F–ψ₀ relations for DEM simulations of isotropic and K₀ consolidations are presented in Figure 19(b). Unlike for the triaxial tests, these data showed two distinct relations for isotropic and K₀ consolidations. The best-fit trend of the u_F–ψ₀ relation for isotropic consolidation has an RMSD of 0·046 and the best-fit line for K₀ consolidation has an RMSD of 0·004.

This study combined triaxial tests on YS sand and 3D DEM cubical assembly consisting of ellipsoids under triaxial condition to investigate the effect of consolidation condition on the overall undrained behaviour and its characteristic features. Both density indices and the CSSM framework were used to develop the correlations for characteristic features of undrained behaviour (often referred to as macromechanical behaviour). The micromechanical measures such as coordination number (C_N) and fabric anisotropy (F_vM) were evaluated to understand how the micromechanical behaviour affect the macromechanical behaviour. The major findings from this study are listed in the following.

• For the same preconsolidation void ratio (e), the consolidation paths have a significant influence on the compressibility and post-consolidation void ratio (e₀) for both triaxial and DEM simulations. The K₀ consolidation path achieved higher compression than that of corresponding isotropic consolidation path. The DEM simulation showed that for isotropic loading – that is, $σ1′=σ2′=σ3′$ – the C_N increased with $σ1′$ to achieve a stable state and exhibited lower compressibility, whereas C_N initially decreased with $σ1′$ for K₀ consolidation, which may have contributed to higher compressibility. Similarly, the fabric anisotropy (F_vM) evolved significantly during K₀ consolidation, whereas F_vM reduced during isotropic consolidation.

• For loose specimen, the instability or static liquefaction was triggered just after peak deviatoric stress (q_peak) at instability stress ratio (η_IS). The K₀-consolidated triaxial test reached the triggering of static liquefaction – that is, η_IS at a much smaller deviatoric strain (ϵ_q = 0·5%) than the deviatoric strain isotropic consolidated tests (ϵ_q = 2–3%). A similar behaviour was observed for the ESP in the DEM simulation. However, there was no conclusive evidence from the DEM simulation that static liquefaction would trigger at a much smaller ϵ_q for K₀ consolidation than isotropic consolidation.

• For a dense specimen, the ESP for both triaxial and DEM simulations after K₀ and isotropic consolidations reached critical state stress ratio (M). However, the ESP for DEM simulation rose above the M-line and then turned back to M towards the end of the simulation. This behaviour is not very common for experimental studies.

• The η_IS for loose specimen increases with relative density (D_r) for triaxial tests. Similarly, η_IS increased with density index (D_I) for DEM simulation. The trend was not dependent on K₀ or isotropic consolidation for both triaxial and DEM simulations. The micromechanical measures – for example, C_N and F_vM at instability – also increased with D_I. Similarly, the characteristic features of the PT state also showed correlations with D_r and D_I for triaxial and DEM simulations, respectively.

• The critical state lines (CSLs) for YS sand were different in the e–log(p′) space for isotropic and K₀ consolidations, where the CSL of K₀ data was higher than that of isotropic data. However, the CSL for DEM simulation for isotropic and K₀ consolidations was unique in the e–log(p′) space. The micromechanical measures – for example, C_N and F_vM at CS – also showed a unique relation with p′.

• The η_IS for loose specimen decreased with ψ₀ for both experiments and DEM simulations. However, the η_IS–ψ₀ relation was different for isotropic and K₀ consolidation types. The unique relation was observed between C_NIS and ψ₀ irrespective of consolidation type. The F_vM–ψ₀ relation was dependent on the consolidation type as observed for the η_IS–ψ₀ relation. This suggests that the difference in the η_IS–ψ₀ relation for different consolidation was not due to the change of C_NIS but for F_vM.

• The η_IS and u_F showed good correlation with density indices (D_r and D_I) and ψ₀.

This study was based on the experimental data for YS sand and DEM simulations of the 3D cubical assembly of ellipsoid particles with linear contact model. There were assumptions made on estimating PT states, extrapolation for CS and others. Therefore, reader should be cautious in generalising the outcome of this article.

The authors would like acknowledge Professor Matthew R. Kuhn of University of Portland, USA, for his discrete-element method software Oval. The second and third author would like to acknowledge the financial support of the University Presidents Scholarship and the University of South Australia Postgraduate Research Award from University of South Australia.