Micromechanical investigation of liquefaction of granular media by cyclic 3D DEM tests

Liquefaction is a very common natural phenomenon that occurs in saturated loose sands during an earthquake, causing very significant destruction to the built environment. Since the two main earthquakes in 1964 (i.e. Niigata earthquake, Japan and Alaska earthquake, USA), much research has been carried out to investigate the physical conditions leading to liquefaction using laboratory experiments, model testing and analytical/numerical methods. As regards the numerical methods, the analysis of the undrained response is traditionally based on a continuum approach. This approach, however, is unable to provide information on the physics occurring at the grain scale. On the contrary, the discrete-element method (DEM) (Cundall & Strack, 1979) provides detailed information of the micromechanical behaviour at the grain scale. Research on undrained conditions (especially for loose sand) has been performed since the works of Bishop (1966, 1971) and Castro (1969). However, significant knowledge gaps remain for the liquefaction of granular materials induced by cyclic loadings.

Liquefaction of gravels is a phenomenon observed in nature, but much less frequently than for sands and silts. The behaviour of gravel under cyclic loading has been less studied than sands and silts due to experimental difficulties involved in testing materials with large particles. The small amount of laboratory information available (e.g. Kokusho, 2007) suggests that gravels are more dilatant than sand and silt for the same level of compaction.

As is well summarised in Carrera et al. (2011), in the geotechnical community there are two definitions of liquefaction: for some it simply indicates strain softening and compressive undrained behaviour, so that a sample under stress control would fail uncontrollably at the peak deviator stress (e.g. Sladen et al., 1985); for others it occurs only when the granular material is subject to zero effective stress and zero strength (e.g. Yamamuro & Covert, 2001). The first definition stems from the practical necessity of geotechnical engineers to design mitigation measures to prevent the occurrence of liquefaction: a criterion in terms of state of stresses for the onset of liquefaction is useful for the design of such measures. But from a phenomenological point of view, reaching a stress state that leads the sample to experience uncontrollable failure cannot be considered liquefaction because it can be argued that the sample is still in a solid-like state, whereas the word liquefaction suggests the sample having experienced a state change: from solid-like to liquid-like. With regard to the second definition, it is noted that even in the liquefied state there are interparticle collisions, which, together with the stress contribution because of the fluctuating kinetic energy density, so-called Reynolds stresses, leads to a state of stress that is never actually zero but simply a rather small value in the context of traditional soil mechanics (Gong et al., 2012). In the granular physics community more attention is paid to the phenomena occurring at the microscale, so liquefaction (the word unjamming tends to be used instead) indicates the transition from a solid-like behaviour to a liquid-like behaviour. A solid-like behaviour is identified with persistent load-bearing contact force chains, whereas liquid-like behaviour is identified with the absence of such persistent load-bearing force chains (Campbell, 2002; Chialvo & Sundaresan, 2013). In this paper the phrase ‘onset of liquefaction’ is employed to indicate the mean effective stress becoming transiently equal to zero, which occurs with the medium still exhibiting a solid-like behaviour – that is persistent load-bearing contact force chains are present, as will be shown later on, and full liquefaction to indicate the completion of the transition process from solid-like to liquid-like behaviour that requires some significant extra straining to be imposed (see the later section entitled ‘Results and analysis’).

In the paper a micromechanical investigation of liquefaction performed by way of three-dimensional (3D) DEM simulations of cyclic undrained triaxial tests in a periodic prismatic cell is illustrated. The purpose of the investigation is to shed light on the micromechanical signatures of the onset of liquefaction and the establishment of full liquefaction for granular media subject to seismic conditions.

The undrained behaviour of loose samples was first investigated by way of the DEM by Thornton & Barnes (1986), Ng & Dobry (1994), Sitharam et al. (2002) and Zhang (2003), who all employed dry samples loaded under constant volume – that is isochoric tests. Bonilla (2004) made a comparison between DEM isochoric simulations without pore fluid and DEM undrained simulations with fluid coupling, using two-dimensional assemblies of elliptical particles, reaching the conclusion that ‘both methods provide similar results’. More recently, Gong et al. (2012) performed monotonic 3D DEM axisymmetric undrained compression triaxial tests in a prismatic periodic cell maintained at constant volume, measuring effective stress paths that were qualitatively similar to published physical experimental results. Also, several recent works (Sitharam & Dinesh, 2003; Sitharam & Govinda Raju, 2007; Sitharam et al., 2009; Soroush & Ferdowsi, 2011; Xu et al., 2015) have successfully simulated undrained conditions by 3D DEM analyses without including the fluid in the model. The tests presented in this paper are of this type.

To perform the simulations described here, the open source software YADE (Šmilauer et al., 2015) was employed. Several arrangements of particles of different grain size distributions (GSDs) were considered and, for each one of them, cyclic undrained triaxial test simulations were performed at different porosities. The DEM simulations were performed on dry granular assemblies – that is the fluid phase was ignored and all stresses were calculated from the orientational distributions of forces at the contacts between particles. Consequently, the effective stresses were calculated directly. This contrasts with laboratory experiments in which the total stresses and the pore water pressure are measured, and the effective stresses are obtained indirectly using Terzaghi's effective stress equation. This paper focuses on one such arrangement.

The work presented in this paper stemmed from the need to assess the liquefaction potential of a real gravelly deposit to be employed as the foundation of a concrete-faced rockfill dam, called Punta Negra, located in the province of San Juan, Argentina. DEM simulations were initially performed for four different GSDs to study the influence of the sand content on the liquefaction resistance when subjected to cyclic loading under undrained conditions. Since the qualitative trends exhibited by these samples did not differ, here only the simulations run for one GSD are reported. The sample consists of 10 000 elastic spheres with an initial porosity of 0·416 and mean stress of 250 kPa; the GSD is shown in Fig. 1.

The contact law used in the tests presented here is the Hertz–Mindlin (HM) no micro-slip solution (Mindlin, 1949) which was implemented in YADE by Modenese et al. (2012). The model considers non-linear compression, linear rotational resistance, non-linear tangential shearing with a Mohr–Coulomb slip criterion, and non-linear viscous damping in both the normal and tangential directions. The normal force is computed at each time step from the current overlap between two grains in contact. The tangential force, on the other hand, is updated incrementally.

The elastic normal force is expressed by the Hertzian solution (Johnson, 1985), as

where $E∗$ is the equivalent Young's modulus (equation (2)); $R∗$ is the equivalent particle radius (equation (3)); and $uN$ is the normal relative displacement of the particles at the contact.

where $νa$ and $νb$⁠, $Ea$ and $Eb$⁠, $Ra$ and $Rb$ are the Poisson ratios, elastic moduli and radii of the two particles in contact, respectively. The tangential force is a function of the relative tangential displacement increment $Δδ$ at the contact (Mindlin, 1949), and it is calculated incrementally, as shown by equations (4) and (5).

where $FTeliandFTeli−1$ are the elastic tangential forces at time steps i and i − 1, respectively; $kTi$ and $kTi−1$ are the tangential stiffnesses at time steps i and i − 1, respectively; and $Δδ$ is the relative tangential displacement within the interval.

Finally, the normal and tangential contact stiffnesses are given by

where $G∗$ is the equivalent shear modulus

where $Ga$ and $Gb$ are the shear moduli of the two particles in contact. Sliding is also included in the model, limiting the maximum shear force with the Mohr–Coulomb  failure criterion – for example, $FSMAX=μFN$⁠, where $μ$ is the coefficient of interparticle friction. A non-linear viscous damping term, already implemented in YADE (Modenese et al., 2012), is included at the contacts. The damping coefficients used are

where $m∗$ is an equivalent mass computed as the harmonic average of the masses of the two particles in contact, equation (11); and $ξN$⁠, $ξT$ are fractions of the critical damping coefficients of a linear spring–dashpot system of stiffnesses $KN$ and $KT$⁠, respectively.

Finally, the normal and tangential contact viscous forces are computed

The total normal force, $FN=FN,el+FN,ν$⁠, is then employed to compute the maximum tangential force allowed by the Mohr–Coulomb criterion. The normal viscous component could exceed the value of the elastic repulsive force $FN,el$⁠, producing a negative normal force. This would be a non-physical, meaningless situation; hence in this case $FN$ is set to zero, as suggested by Schwager & Pöschel (2007) for a linear spring–dashpot model.

To account for the effect of particle shape, rolling resistance is applied for both rolling and twisting moments incrementally, as for the tangential force

where $MROLL$ and $MTWIST$ are the rolling and twisting moments; $KR$ and $KTWIST$ are the linear rolling and twisting coefficients of the interaction; and $θR$ and $θTWIST$ are the relative rolling and twisting rotations between particles in contact. Relative rotation–moment relationships were first proposed by Iwashita & Oda (1998). Here they are employed to account for the fact that the directions of the normal forces exchanged between real gravel (non-spherical) particles do not intersect the particle centres, thus generating torques that do not occur in the case of spherical particles (Ai et al., 2011). The rolling moment between two contacting particles is limited to a maximum value given by $μR∗Rmean$⁠, whereas no limit is applied to the twisting moment.

The triaxial simulations were run under quasi-static conditions. To check the presence of quasi-static conditions, the so-called inertial number, I (Da Cruz et al., 2005), was employed

where $ε˙$ is the shear strain rate applied to the periodic cell; $d$ is the average particle diameter; $ρ$ is the particle density; and p′ is the effective pressure applied to the arrangement. The density of all the grains is set to a value of 2600 kg/m³ and no gravity is applied.

According to Radjai & Dubois (2011), the inertial number must be lower than 10⁻³ to guarantee that simulations are in the quasi-static regime. Consequently, the maximum deviator strain rate initially applied in the tests was 0·1 s⁻¹. However, Radjai & Dubois (2011) employed an almost uniform particle size distribution in their tests, so the I value to be considered as threshold for quasi-static regime needs to be ascertained afresh for the specific PSD considered here. Hence monotonic undrained (isochoric) triaxial tests were performed on four different samples, with the deviator strain rate varied from 0·016 s⁻¹ to 8·192 s⁻¹. In Fig. 2 the maximum stress against the deviator strain rate applied is plotted for the tests performed. From Fig. 2 it emerges that a strain rate lower than 0·1 s⁻¹ is enough to avoid the influence of inertial effects on the stress–strain behaviour during monotonic loading. The same is expected for cyclic deformation.

The total stress tensor of a soil volume can be decomposed into the sum of the effective stress tensor due to the forces at the contacts between the particles and the pressure tensor of the fluid within the soil voids. The effective stresses can be calculated from the summation of the contact forces and the fluctuating part of the particle velocities (Thornton, 2000)

where $m$ is the mass of each particle; $u~i$ and $u~j$ are the fluctuating velocities of each particle; V is the volume of the sample; l is the distance between the centres of the two contacting spheres; and F is the contact force. The first term in equation (17) represents the Cauchy effective stress tensor, $σ′ijC$ while the second term corresponds to the Reynolds effective stress tensor, $σ′ijR$⁠. Under quasi-static conditions, Reynolds stresses are negligible, but they have been calculated in order to check whether they remain negligible following the onset of liquefaction (see Fig. 10 later in the paper).

Following Salot et al. (2009) and Modenese (2013), samples were prepared according to a three-step procedure: (a) generation of non-overlapping particles randomly in space with an intergranular friction angle of 40°; (b) compression of the particles up to the development of a cell pressure of 250 kPa; (c) gradual decrease of the friction angle at constant pressure until the target porosity of 0·416 was reached. The intergranular friction angle is then returned to the value assigned at particle generation. To ensure quasi-static conditions throughout the preparation phase, the ratio of the average out-of-balance force of all the particles to the average force at all the contacts was maintained lower than 0·01.

Then, the assemblies were subjected to different values of cyclic deviator strain, ranging from ±0·2% to ±1%. The volume was kept constant by maintaining a constant value of the axial strain rate of 0·1 s⁻¹, and half that value (0·05 s⁻¹) in the two other orthogonal directions. Since it is not possible to measure total stresses during an isochoric (constant volume) DEM test, the total stress path during cyclic shearing is unknown and therefore the pore pressure too. The assembly was considered to be liquefied when the ratio of the mean effective stress of the assembly, p′, to the initial confining pressure applied, p₀ = 250 kPa, fell below 0·01.

The main parameters used for all the DEM simulations are: Poisson ratio, $ν=0⋅1$⁠; Young's modulus, $E=1×109N/m$⁠; friction coefficient, $μ=tanϕ$⁠, with ϕ = 40°; coefficient of rolling stiffness, $KR=7500J/rad$⁠; coefficient of twisting stiffness, $KTWIST=5000J/rad$⁠; coefficient of rolling friction, $μR=0⋅8$⁠; normal coefficient of restitution, $eN=0⋅9$⁠; and tangential coefficient of restitution, $eT=0⋅9$⁠. Given that the values of Poisson ratio and Young's modulus are the same for all the particles within the arrangement, the ratio between the tangential and normal contact stiffness is 0·947 (⁠$KT=0⋅947KN$⁠).

Figure 3 shows the observed macroscopic behaviour of the tested sample under cyclic shear strain. The figure shows, in Fig. 3(a), that the amplitude of the Cauchy stress component oscillations gradually reduces to zero when the normalised effective pressure, shown in Fig. 3(b), becomes null. Fig. 3(c) shows the deviator stress plotted against the deviator strain. It can be seen that the inclination of the hysteresis loops decreases as the simulation progresses as a consequence of the decrease in p′, as shown in Fig. 3(d).

The variables analysed are the mechanical coordination number, the geometrical coordination number, the redundancy index, number of contacts, size and number of clusters, Reynolds stresses and inertial number.

The coordination number is defined as the average number of contacts per particle. Thornton (2000) introduced the so-called mechanical coordination number Z_m to exclude particles that do not contribute to load bearing – that is floaters and rattlers, as

where C is the total number of contacts within the particle arrangement; $N$ is the total number of particles; $N1$ is the number of particles with only one contact; and $N0$ is the number of particles with no contacts. When the mechanical coordination number was introduced by Thornton (2000) it was thought that, since the particles with only one contact were transient collisions, they did not significantly contribute to the stress tensor and therefore should be ignored. However, when considering isostaticity at any moment in time, these contacts cannot be ignored. Consequently, it is more relevant to consider the geometrical coordination number, Z_g (Thornton, 2015), which is defined as

If the total number of constraints is equal to the total number of degrees of freedom, then the system is said to be isostatic on average. In three dimensions, the number of degrees of freedom per particle is six (three translations and three rotations). In the simulations reported by Thornton (2000), the number of constraints per contact was three (the normal force and the two components of the tangential force). In the current simulations, three additional constraints are present due to the two components of the rolling resistance and the twisting resistance; therefore, the system is isostatic on average if 6C = 6(N − N₀) and substituting into equation (19) the threshold coordination number becomes Z_g-iso = 2.

An alternative, and perhaps a more reliable, parameter that can be used to assess the stability of the system is the redundancy index (Gong et al., 2012) or redundancy ratio (Zhou et al., 2017), redundancy factor (Kruyt & Rothenburg, 2009) or constraint ratio (Cundall & Strack, 1983), which takes into account the number of contacts that are sliding and/or rolling at any moment in time. The redundancy index is defined as the ratio of the total number of constraints to the total number of degrees of freedom of the active particles. In the current simulations no limit was specified for the twisting moment, therefore the redundancy index is here defined as

where $CNSR$ is the number of contacts that are neither sliding nor rolling (six constraints given by three forces and three  moments); $CS$ is the number of contacts that are sliding (four constraints given by the normal force and three moments); $CR$ is the number of contacts that are rolling (four constraints given by the three force components and one twisting moment); $CSR$ is the number of contacts that are sliding and rolling simultaneously (two constraints given by the normal force and the twisting moment); $N$ is the total number of particles; and $N0$ is the number of particles with no contacts. The condition $IR=1$ identifies isostatic conditions on average for the granular system. Fig. 4 shows that the onset of liquefaction from a stress state point of view – that is, the occurrence of $p′/p0=0$ – coincides with the redundancy index reaching unity, which also confirms the correctness of the definition of I_R.

The calculation of stresses is based on contact forces and particle velocities (see equation (17)). It is important to note that the number of particles employed in the simulations is limited in comparison with nature, 10 000 rather than billions. This implies that the stresses measured in the periodic cell exhibit some small oscillations, stemming from the limited number of particles of the cell, which would not be seen in a much larger sample. Therefore, the stress threshold for the onset of liquefaction p′/p₀ = 0 is considered to be reached when p′/p₀ < 0·01. Using this criterion, the onset of liquefaction is identified to take place at about 5 × 10⁶ time steps (see the dotted and dashed vertical line in Figs 4 and 5). After this point the effective pressure increases due to the cyclic nature of the loading, until after a few more cycles p′/p₀ remains stable (null). From Fig. 4 it emerges that at the identified onset of liquefaction the geometric coordination number is very close to 2, indicating the sample to be isostatic on average. This is a rather remarkable finding since it indicates that the onset of liquefaction from a stress state point of view – that is, the occurrence of p′/p₀ = 0 – is associated with the geometrical coordination number approaching the value of the threshold geometrical coordination number expressing the isostatic condition for the granular system which has not previously been shown for assemblies of granular particles with rolling resistance.

Owing to cyclic straining, fluctuations of p′/p₀ continue to occur for a while, as has already been reported by Ng & Dobry (1994). Also from both Figs 4 and 5 it is clear that not only p′/p₀ but also the geometrical coordination number and redundancy index keep varying after the onset of liquefaction for some significant further shearing. This indicates that liquefaction does not occur instantaneously but as a process whereby the behaviour of the granular medium transitions from solid-like to liquid-like. The key question left open is: when does the transition end, or, in other words when is full liquefaction established?

Looking at Figs 4 and 5 it is possible to conclude that, once the fully liquefied condition is established, Z_m, Z_g and I_R remain constant at approximately 2, 1 and 0·35, respectively, which occurs at some point after 6 × 10⁶ time steps. However, it is not straightforward to establish when they reach their steady value owing to the strong oscillations exhibited. In the quest to clarify the solid-like to liquid-like transition, it is useful to analyse the evolution of contacts and clustering of the particles.

In Fig. 6 the total number of contacts is plotted together with the type of contact – that is the fraction of rolling, sliding and rolling–sliding contacts. It is apparent that the number of contacts decreases during the application of cyclic deformation and, following the onset of liquefaction, a sudden decrease of the contacts to a value oscillating between 200 and 2500 is visible. Then contacts continue to decrease to a value lower than 100 after more than 6 × 10⁶ time steps. In correspondence with the abrupt lowering of the number of contacts, a sudden increase in the fraction of particles that are either sliding, rolling, or sliding and rolling at the same time occurs.

So far observations on the number of contacts (Fig. 6) and redundancy index suggest full liquefaction is established after 6 × 10⁶ time steps. To elucidate further the transition from solid-like to liquid-like behaviour, an analysis of how particles cluster together was undertaken.

A cluster is a group of particles in contact at any given time step throughout the simulation. The cluster size is defined by the number of particles in the cluster, and the number of clusters includes singlets. Fig. 7 shows the evolution of the number of clusters and the size of the largest cluster. Initially, the largest cluster consists of 5000 particles and there are 5001 clusters. This means that all the particles that are not in the largest cluster are singlets – that is floaters and rattlers. As the mean effective stress decreases, the magnitude of the normal contact forces reduces. With continued decrease in the mean effective stress, an increasing number of contact normal forces reduce to zero and contacts are lost and the maximum cluster size gradually reduces at an increasing rate.

To identify the establishment of full liquefaction, in Fig. 8 the evolution of the maximum cluster size, the number of particles per cluster, is plotted rescaled after the onset of liquefaction – that is after 5 million time steps. Looking at this figure and its inset, it can be concluded that full liquefaction is established after 6·6 × 10⁶ time steps with the maximum cluster size then oscillating between three and 20 particles (see Fig. 8 inset). This shows that after full liquefaction takes place, only a few particles are in contact, with the geometrical coordination number and redundancy index oscillating around stable averages of 1 and 0·35, respectively (Figs 4 and 5, respectively).

Figure 9 shows the particle assembly at the onset of liquefaction – that is after 5 million time steps, when the mean effective pressure is momentarily null and the redundancy index, I_R = 1. The complete arrangement of particles, the largest remaining cluster (1229 spheres) and the force chain network within the largest cluster are plotted in Figs 9(a), 9(b) and 9(c), respectively. Interestingly, the figure suggests that the loss of contacts is uniformly distributed throughout the assembly and that there is no break-up of the large cluster into smaller clusters. In addition, this appears to be a good example of diffuse failure (Darve, 1996). It emerges that, during degradation of the initial large cluster, only singlets are created, increasing the number of ‘rattlers’.

In Fig. 10 both the mean Reynolds stress and the inertial number are plotted. It emerges that Reynolds stresses are negligible throughout the simulation, even after full liquefaction is established. With regard to the inertial number, this maintains at a value below 10⁻³, which is acknowledged as a threshold to indicate the presence of quasi-static conditions, approximately 5 × 10⁻⁴ approaching liquefaction. But following the onset of liquefaction, with p′/p₀ ∼ 0, a rather sudden increase is measured and then a linear increase over time to values well above the quasi-static threshold to reach ∼0·3, which is consistent with a collisional flow. However, it needs to be pointed out that after the onset of liquefaction, the inertial number may well cease to be a reliable indicator of the inertia of the system due to the pressure p′ approaching 0 (see equation (16)).

The paper has focused on the microscopic features of cyclic strain-induced liquefaction. A key finding is that the onset of liquefaction in the traditional geotechnical sense – that is, p′ = 0 – occurs when the granular medium becomes on average isostatic – that is, the geometrical coordination number is equal to 2 and the redundancy index is 1.

As the mean effective stress decreases, the magnitude of the contact normal forces reduces and an increasing number of contacts are lost, with the maximum cluster size gradually reducing at an increasing rate. The loss of contacts is uniformly distributed throughout the sample and there is no break-up of the large cluster into smaller clusters. This is different from monotonic undrained tests on loose specimens in which the isostatic state is reached soon after the deviator stress reaches a maximum value (Gong et al., 2012).

From the presented investigation it turns out that the transition from solid-like to liquid-like behaviour or, in other words, from the onset of liquefaction to full liquefaction, can be ascertained rather well by analysing the evolution of the geometrical coordination number, redundancy index, particle contacts and clusters. In particular, full liquefaction is associated with the point in time when all these variables start to oscillate steadily around constant values.

The authors gratefully acknowledge the contribution of the Marie Skłodowska-Curie Research and Innovation Staff Exchange (RISE), through the GEO-RAMP project (grant agreement no: 645665).

$C$

total number of contacts

$CNSR$

number of contacts that are neither sliding nor rolling

$CR$

number of rolling contacts

$CS$

number of sliding contacts

$CSR$

number of simultaneously sliding and rolling contacts

$d$

average particle diameter

$E∗$

equivalent Young's modulus

$Ea,Eb$

elastic Young's moduli of particles a and b

$eN$

normal coefficient of restitution

$eT$

tangential coefficient of restitution

F

contact force

$FN$

total normal force

$FN,el$

elastic normal force

$FN,ν$

viscous normal force

$FT,ν$

viscous tangential force

$FTMAX$

maximum tangential force

$FTeli,FTeli−1$

elastic tangential forces at time step i and i − 1

$G∗$

equivalent shear modulus

$Ga,Gb$

shear elastic modulus of particles a and b

$I$

inertial number

$IR$

redundancy index

$KN$

normal stiffness

$KR$

rolling coefficient

$KT$

tangential stiffness

$KTWIST$

twisting coefficient

$kTi,kTi−1$

tangential stiffnesses at time steps i and i − 1

l

distance between the centres of two contacting spheres

$MROLL$

rolling moment

$MTWIST$

twisting moment

$m∗$

equivalent mass

$ma,mb$

mass of particles a and b

$N$

total number of particles

$N0$

number of particles with no contacts

$N1$

number of particles with only one contact

p₀

initial confining pressure applied

p′

effective pressure

$R∗$

equivalent particle radius

$Ra,Rb$

radius of particles a and b

$Rmean$

mean radius between particles

$u~i,u~j$

fluctuating velocities of particles i and j

V

volume of sample

v_N,rel

relative normal velocity

v_T,rel

relative tangential velocity

$Zg$

geometrical coordination number

$Zg_iso$

threshold coordination number

$Zm$

mechanical coordination number

$Δδ$

relative tangential displacement

$ε˙$

shear strain rate

$ηN$

normal damping coefficient

$ηN,crit$

critical normal damping coefficient

$ηT$

tangential damping coefficient

$ηT,crit$

critical tangential damping coefficient

$θR$

relative rolling rotation

$θTWIST$

relative twisting rotation

$μ$

coefficient of interparticle friction

$μR$

coefficient of rolling friction

$νa,νb$

Poisson ratios of particles a and b

$ξN$

fraction of the critical normal damping coefficient

$ξT$

fraction of the critical tangential damping coefficient

$ρ$

particle density

$σ′ij$

total effective stress tensor

$σ′ijC$

Cauchy effective stress tensor

$σ′ijR$

Reynolds effective stress tensor

$ϕ$

friction angle