Discrete element modelling of material non‐coaxiality in simple shear flows

We investigate the quasi‐static simple shear flow of a two‐dimensional assembly of cohesionless particles using discrete element method (DEM) simulations. We focus on the unsteady flow regime where the solid would experience significant evolution of stresses, mobilised shear strength and dilation. We construct the DEM model using a discretised‐wall confined granular cell where the apparent boundary is allowed to dilate or contract synchronously with the confined solid. A rather uniform simple shear field is achieved across the whole assembly, which benefits rheological studies in generalising constitutive laws for continuum methods. We examine two aspects of the simple shear behaviour: macroscopic stress and strain rate evolution, particularly the non‐coaxiality between the principal directions of the two; and micromechanics such as evolution of fabric. For an initially anisotropic specimen sheared under constant normal pressure, the direction of principal stress rotates towards that of the principal strain rate, gradually reducing the degree of non‐coaxiality from about 45° to fluctuating around 0°. The rate in approaching coaxiality is slower in samples with larger initial porosity, stress ratio and mean stress. Generally, a faster rate in approaching coaxiality in simple shear is observed in a more dilatant sample, which often shows a larger degree of mobilised fabric anisotropy, suggesting the possible important role of instantaneous internal friction angle. The evolution of principal fabric direction resembles that of the principal stress direction. © 2013 The Authors. International Journal for Numerical and Analytical Methods in Geomechanics published by John Wiley & Sons, Ltd.


INTRODUCTION
Loading conditions encountered in geotechnical engineering invariably involve principal stress rotation [1], particularly when cyclic loading is involved, for example, under earthquakes, vehicular traffic and sea waves. Upon continuous principal stress rotation, permanent strains may accumulate steadily even if the magnitudes of the effective principal stresses are kept constant [1]. It is therefore important to examine the effects of principal stress rotation on the deformation of soils and to incorporate such effects into constitutive models if reliable predictions are to be made.
Extensive evidence from laboratory has shown that the principal axes of stress and strain rates are generally not coincident during stress rotation (e.g. [2][3][4][5][6][7][8][9][10][11]). This 'non-coaxiality' phenomenon is an important part of the elastoplastic theory as the postulate of coaxiality made in coaxial plasticity [12] is only valid for isotropic media and problematic when applied to anisotropic media. Therefore, the constitutive relationship cannot be sufficiently formulated in the principal stress space unless non-coaxiality is taken into account [13]. Some numerical evidence of the importance of *Correspondence to: Jun Ai, Nottingham Centre for Geomechanics, The University of Nottingham, University Park, Nottingham NG7 2RD, U.K.
including non-coaxiality in engineering applications such as soil-structure interaction and solid flow in silos have been reported recently (e.g. [14,15]).
Despite numerous experimental and numerical studies, the understanding of the non-coaxial behaviour, particularly its micro-mechanism, remains limited, and the non-coaxial plasticity theory has not yet been satisfactorily developed [16]. The non-coaxial soil stress-strain response can be studied in the laboratory such as with simple shear apparatus [2,8], directional shear cell [5,6], torsional hollow cylinder [17] and two-dimensional (2D) general stress apparatus [18,19]. These tests allow various degrees of control of the rotation of the principal axes of stress or strain. As well known, they suffer various amounts of inhomogeneity in achieved stress and strain fields [13,19,20]. From the viewpoint of rheological study of granular materials using a representative volume element, obviously, the global uniformity of the geometrical and mechanical properties across the whole element is pursued. Another inherent problem of the laboratory studies is that they cannot create accurate replicates of the physical systems on site.
Complementary numerical experiments of stress rotation by discrete approaches such as discrete element method (DEM) simulations, though limited, have been carried out (e.g. [21][22][23]). It is expected that the control of the boundary conditions could be easier in the idealised numerical systems, and the richness of output data, particularly micromechanical information, may offer further insights that help construct more physically realistic, such as fabric-based, non-coaxial plasticity models [16]. At present, most existing non-coaxial constitutive models are based on continuum mechanics and mostly heuristic assumptions (e.g. [24][25][26]). This paper presents a DEM examination of non-coaxial behaviour in a granular cell undergoing simple shear. Simple shear is one of the most common deformation patterns of particulate systems. For example, simple shear is generally appreciated as the typical deformation mode in localised failure zones (shear bands) of granular materials and in shaking level grounds under seismic shear waves. It is also the dominant deformation mode in plane shear flows of granular solids [27]. Most of the earlier DEM studies of simple shear mainly focused on steady state (or critical state) behaviours (e.g. [27][28][29][30]). However, in such large strain regimes, the volumetric strain has been fully mobilised, and the directions of principal stress and strain rates have already developed to be largely coincident. In contrast, the relatively small shear strain regime is more important in geomechanics where deformation and modulus are important factors for practical design. Previous DEM simulations investigating non-steady regime of simple shear are largely limited to strength and volumetric behaviour (e.g. [31][32][33][34]). Only very few studies consider the non-coaxiality [22,35,36].
In past DEM simulations, the shear cell commonly consisted of four geometrically rigid line-walls forming a parallelogram, which is similar to that of the laboratory simple shear devices where the soil specimens are enclosed by rigid platens. Consequently, some inherent limitation of the laboratory device also exits in the numerical model; for example, the state of stress and strain in the cell is seldom homogeneous because of the restraints of solid dilation or compaction near the wall boundaries.
We report a systematic investigation of simple shear behaviour focusing on non-coaxiality, using the proposed 'discretised-wall confined granular cell'. In this cell, the boundary may deform strictly conforming to the applied strain rates as a result of its discrete nature, which avoids any spurious boundary shear stress that distorts the contraction or dilation of the confined granular solid. Rather uniform stress and strain fields have been achieved, which facilitates a more reliable homogenisation for rheological study. The effects of initial porosity, stress ratio, lateral pressure ratio, mean stress on the simple shear responses are explored. Both macro and micro-mechanical information are analysed.
The results presented here are based on 2D systems with idealised circular particles. The threedimensional case and the effects of particle shape on the non-coaxiality development, and that of partially restraining particle rolling via introducing rolling resistance [37][38][39][40][41] are also studied, but the results will be presented elsewhere.

Granular cell configuration
2.1.1. Previous line-wall confined cell. A typical simple shear cell is of a parallelepiped configuration (parallelogrammic in 2D, cf. Figure 1b). For a uniform stress state to occur under simple shear mode, equilibrium demands that complementary shear stresses be developed on the vertical sides of the element normal to the plane of deformation. This consequently requires frictional boundary conditions. The configuration should not prohibit the enclosed solid from any dilation or contraction under shearing. To simultaneously satisfy both conditions remains a big challenge in the laboratory; thus, certain simplifications have to be made. For example, in the Cambridge Simple Shear Apparatus, the two side platens are made smooth, and in the Direct Simple Shear Apparatus, the sample is laterally wrapped by a flexible membrane. Consequently, it has been shown that the state of stress or strain in neither of these apparatuses is homogeneous [20,[44][45][46]. It is worth noting that the 1γ2ε apparatus developed by Joer et al. [19] significantly alleviates such limitations, in which the confining platens are novelly replaced by extendable articulated arms instead of rigid platens. The numerical hexagonal configuration used by Li et al. [47] also shows some improvement in alleviating the arching effect near the corners. Figure 1 shows examples of force chain network distributions in a granular sample enclosed by four frictional line-walls as in earlier DEM studies. The line-walls are geometrically non-extendable, hence cannot dilate or compact together with the solid. This may result in spurious relative movement between solid and boundaries, for example, producing either excess or insufficient shear stress near frictional boundaries (e.g. Figure 1a). Regarding simple shear simulations, there appears no report of any successful realisation of uniform distribution of force chains. Typically, much weaker force chains are observed near the two acute corners than in the rest of the element (cf. Figure 1b). Also, using too large shear rates is another potential cause of non-uniformity. Moreover, incompatible motions between two pairs of parallel walls due to inaccurate kinematic control may also induce incorrect resultant shear stresses at the boundary.

2.1.2.
Discretised-wall confined cell. In the proposed configuration, we replace the line-walls by pointwalls, that is, enclosing the granular sample by individual rigid points with inter-point gaps small enough to prevent any particle from escaping. Hence, the boundary is made fully discrete in nature. The motion of each boundary point is governed by the target strain rate according to its location: where x i is the coordinate of the boundary point and ε ij is the target strain rate tensor. Therefore, the target strain rate field is strictly met at the boundary, whereas that within the assembly is dependent upon material behaviour. Any general stress, strain or mixed loading paths may be applied to the cell. Prescribed strain rate components may be applied directly, whereas prescribed stress components may be indirectly applied via strain rate components using any servo-control algorithm (e.g. [48]). The volume enclosed by the discrete boundary points may take any shape in principle, no need to maintain parallelepipeds/ parallelograms. Because any corners or narrow regions may enhance local frustrations, it is always advantageous to choose a smooth and convex shape, such as the ellipse used in our 2D simulations ( Figure 2). As the relative motions between the boundary constituents are strictly governed by the prescribed loading modes, this type of boundary configuration falls into a 'rigid' boundary scenario, in contrast to the 'flexible' boundary scenarios such as pressure-controlled boundary [49]. Figure 2 shows the snapshots of our granular cell at various loading stages. The rather uniform distribution of force network observed across the cell suggests a well-achieved representative element possessing largely uniform stress and strain rate fields. It is noted that the improved cell configuration employed here does not completely eliminate the development of nonuniform strain fields. It only aims to mitigate those induced by inappropriate boundary control. Non-uniform strain distribution due to strain localisation is an inherent material behaviour that occurs spontaneously at appropriate states, hence should not (and cannot) be artificially suppressed [22].

Numerical samples
Our numerical samples are composed of about 4200 disks in a 2D system. The particle sizes are uniformly distributed from 2.4 to 3.6 mm, with the small ±20% polydispersity modelled to prevent crystallisation. Hertz-Mindlin model is adopted as the contact law. The initial samples were generated by radius expansion method [43]. Specifically, we first generated the particle clouds within the space enclosed by the boundary wall points. The radii of the particles were scaled down during the generation process. Once the particles were all generated, we restored their radii back to the true sizes, allowing the particles to engage contacts and settle down under certain numerical damping. By varying the contact friction coefficient and the confining stresses during the process of particle settlement, various levels of porosity can be achieved.
Using identical group of particles, we prepared three specimens with different porosities at a low isotropic stress level p = 50 kPa, namely specimen S1 (dense), S2 (mediate dense) and S3 (loose). Each sample was then subjected to isotropic and anisotropic consolidations, with their loading paths depicted by dashed and solid lines, respectively in Figure 3. For example, the cell status shown in Figure 2a corresponds to a stage on the horizontal dashline, whereas that of Figure 2b to a vertical solid line. A group of derivative samples were therefore prepared from each specimen, with their stress states denoted as solid disks in Figure 3. These derivative samples differ in either mean stress p or stress ratio η = q/p, where q = (σ 1 À σ 2 )/2 is the deviatoric stress and p = (σ 1 + σ 2 )/2 is the mean stress; σ 1 and σ 2 are major and minor principal stress, respectively. Simple shear tests were then carried out upon these derivative samples. This parameterisation allows us to investigate the effect of initial mean stress and stress ratio on the behaviour of subsequent simple shear.
The stress tensor σ ij across the whole cell is obtained by averaging the stress tensor of each particle σ ij k ð Þ over the N p particles contained in the cell [43] where V and V (k) are the total volume of the cell and the volume of particle k, respectively. The stress tensor of each particle was calculated from the contact forces f c ð Þ i acting on the particle  i are the locations of contact point and particle centroid, respectively. The aforementioned method (Eq. 2) is practically equivalent to that given by Drescher and de Josselin de Jong [3] and Cundall et al. [50].
Our study focuses on the quasi-static flow regime. The local damping [43] is adopted to sufficiently dissipate the kinetic energy of the system. The simple shear rate is controlled to limit the inertial number I to prohibit potential inertial effects. The inertial number is a dimensionless measure defined as the ratio between a micro timescale T P and a macro timescale T γ ( [27], Figure 4). It is defined in our study as where d, ρ, P and γ are mean particle radius, assembly bulk density, mean stress and simple shear strain rate, respectively. A flow system approaches quasi-static as I ≪ 1. In our simple shear simulations, the shear rate starts with a small magnitude and gradually increases until the inertial number I reaches an assigned limit. Similar to some other studies (e.g. [30]), a limit inertial number I lim = 10 À3 appears sufficiently small and is hence adopted. For isotropic and anisotropic compression procedures, the loading rate is instead monitored via limiting the unbalanced force ratio ξ = h f unbal i/h f c i defined as the ratio of mean unbalanced force h f unbal i over the mean contact force h f c i. A small magnitude of ξ = 10 À4 was shown to be adequate. The applied strain rate is reduced by multiplying a factor once the unbalanced force ratio exceeds the assigned threshold and increased again if ξ falls below the threshold. Such an adaptive rate control offers a balance between computational cost and accuracy.
The principal DEM parameters used in the simulations are summarised in Table I.

NUMERICAL RESULTS
A granular assembly is a highly nonlinear elastic-plastic system and its behaviour strongly depends on its past loading history. Apart from mean stress and stress ratio, our prepared derivative samples generally also differ in other properties such as porosity. It is well known that soil behaviour such as strength and dilation strongly depend on the initial density even at same stress state. Hence, the pre-simple shear deformation process must be specified for any reliable conclusions on simple shear behaviour.

Isotropic and anisotropic consolidation
The initial pre-consolidated isotropic specimens S1, S2 and S3 at a mean stress level p = 50 kPa have porosities of n = 0.175, 0.183 and 0.201, respectively. After the samples are subjected to isotropic consolidation from p = 50 to 1000 kPa (the dashline path in Figure 3), the samples reach porosities of n = 0.123, 0.128 and 0.141, respectively. It is of interest to note the dimensionless mean stress, p norm = p/G, increases from 2.86 × 10 À6 to 5.72 × 10 À5 within the specified range of mean stress from 50 to 1000 kPa. It has been shown that a system may be regarded as close to the 'rigid-particle limit' if p norm is smaller than the order of 10 À5 [51,52]. It suggests that our systems gradually evolve from a nearly 'rigid-particle limit' to soft-particle systems along the prescribed isotropic compression procedure. Consequently, the reduction in porosity at large stress states should be partially attributed to the overlaps between particles.
At selected levels of isotropic consolidation, the samples are subjected to anisotropic compression with constant mean stress (solid-line paths in Figure 3). The biaxial shearing was terminated once the specified maximum stress ratio η = 0.3 was reached. We observed that a more contractive behaviour is associated with a looser sample or a sample with larger mean stress, which is consistent with the experimental observations [53]. We note that the relative order of the samples with respect to their porosities remains unchanged throughout the biaxial shearing. This maintained sequence sets out the basis for the parametric studies of simple shear behaviour presented in the following.

Simple shear scenario
The prepared derivative samples of specimen S1, S2 and S3 at different levels of mean stress and stress ratio were subjected to simple shear. The vertical pressure σ y was kept constant during the shear, allowing vertical dilation. Figure 5a shows the evolution of normal stress σ x and shear stress σ xy on the shear plane in samples with initial mean stress p 0 = 200 kPa and stress ratio η 0 = 0.2. The accumulated simple shear strain γ is calculated by where ε yx is the shear strain rate in the shear direction and Δt is the size of the timestep. For the initially denser sample derived from specimen S1, the shear stress σ xy curve (thick line) shows a peak at γ = 0.15 before descending to the ultimate value at sufficiently large shear strain, whereas the relatively looser samples show no peak but a gradual increase towards the ultimate value. This behaviour is same to that in biaxial tests or other tests where the mobilised strength depends on the soil density [53]. The normal stress σ x shows a similar trend to that of shear stress. At large shear strain, σ x is largely equal to the vertical stress σ y . As shown in Figure 5b, the curve of stress ratio η exhibits a highly synchronised evolution to that of the shear stress σ xy in Figure 5a. For example, they share the identical positions of jumps, despite that they differ in their starting values: the former always starts from zero, whereas the latter has an initial value η 0 = 0.2. The ultimate magnitude of both the coulomb stress ratio σ xy /σ y and the stress ratio η are around 0.3. The effect of initial density on the volumetric strain (ε v = ε 1 + ε 2 ) is similar to that of other types of shearing tests; that is, a larger initial density results in less contraction at early shearing stage and larger dilation at larger shear strain.
The orientations of major principal stress θ σ and strain rate θ ε are compared in Figure 6, which clearly illustrates the phenomenon of non-coaxiality: a large difference exists between the two in the early stage of shearing. It should be noted that in non-coaxial plasticity theories, 'noncoaxiality' refers to plastic strain rates instead of the total strain rates as presented here. According to the hollow cylinder test reported by Gutierrez et al. [10], the degree of non-coaxiality is slightly smaller during stress rotation if plastic strain rates are considered, whereas the difference diminishes at high stress ratios. A detailed and quantitative evaluation of the effect of the elastic strains requires further study. The two orientations θ σ and θ ε are calculated by where ε x (zero-valued in simple shear), ε y and ε yx are the strain rates measured from the incremental deformation of the boundary. Alternatively, strain rate can also be evaluated from the particle velocity field using least square fittings. Note that the velocity gradient tensor is not symmetric in simple shear mode; hence, a global rigid rotation exists. However, the rigid rotation in a small shearing increment is very small so that it is omitted in Eq. 5b. The data presented in Figure 6 is from a derivative sample of specimen S3, which is a relatively loose sample. When we ignore the spikes, the curve shows that the major principal strain rate orientation θ ε starts from a value slightly below 45°, then gradually goes slightly above and then falls back to 45°. In contrast, the curve of major principal stress orientation θ σ is rather smooth. θ σ starts from zero and gradually approaches the strain rate orientation at large shear strain, which is consistent with earlier studies [22,36]. The predicted non-coaxiality angles are in general agreement with the experimental simple shear test data reported by Roscoe et al. [2]. The degree of non-coaxiality is large (>5°) within the shear strain range of γ < 0.05. Several factors that may affect the non-coaxiality are examined in the following sections. It is of interest to discern the source of the large spikes/fluctuations in the θ ε curves. Our tentative interpretation below is based on the idealised scenario of particle rearrangements sketched in Figure 4. During simple shearing, the particles from two neighbouring layers would roll and slide over each other. Such an inter-crossing process naturally involves numerous local instabilities. Given that the (constant) vertical stress component is servo-controlled, frequent changes in vertical sample dimension are therefore induced. According to the governing Eq. 5b, the orientation of the principal strain rate θ ε strongly deviates from 45°if the relative magnitude of ε yx over ε y is small. This suggests that a sufficiently small shear rate would inherently lead to large fluctuations in θ ε , whereas a large shear rate may on the contrary reduce the fluctuations.
We sheared the same sample with various levels of shear rate by assigning a smaller or larger limiting inertia number I lim . The result with I lim = 0.05, that is, 50 times larger than that in Figure 6a, is shown in Figure 6b. We see that the fluctuations are much smaller under the larger shear rate, which is consistent with the aforementioned corollary. Another characteristic of the fluctuations in the curve of θ ε is that size of spikes pointing towards zero is much larger than those pointing towards 90° (Figure 6a). This may be caused by that the climbing up (i.e. dilating) motions of the particles take longer than that of the snapping down (i.e. contracting) motions (T up p > T down p ; Figure 4). A too large shear rate produces a discrepancy between the strain rate measured from boundary (thin solid-line) and that from internal particle velocity field (thin dashline), as shown in Figure 6b. This indicates that the particles cannot catch up with the boundary when sheared too fast. With sufficiently small shear rates, the two measurements are largely indistinguishable.
The non-coaxiality (θ ε À θ σ ) in samples with different initial densities are compared in Figure 7. The three derivative samples respectively from specimen S1, S2 and S3 are in the same initial stress state of p 0 = 200 kPa and η 0 = 0.2 before shear. Their initial porosities are 0.161, 0.167 and 0.181, respectively. We observe a faster rate in approaching coaxiality in a sample with a larger density. As the strain rate orientations θ ε in the three samples are all largely 45°, the difference should mainly be due to the stress orientations θ σ . The governing Eq. 5a indicates that a larger rate in increasing σ xy and σ x results in a larger rate in increasing of θ σ . This is indeed confirmed by the evolutions of the two stress components shown in Figure 5a where the two stress components increase faster in a denser sample.
We now examine the fabrics of the contact network in the sheared samples. The polar distributions of contact normal direction n at three different loading stages are compared in Figure 8. The symbols denote the measured probabilities of the contact normal orientation P(n) in the samples. The curves are fitted second-order Fourier expansion of P(n) [54] where a c and θ c are the degree of fabric anisotropy and its principal orientation, respectively. The polar distribution is rather homogeneous in all orientations at the isotropic stress state (thin curve), and its degree of anisotropy slightly increases from a c = 0.010 to 0.037 at a stress ratio η = 0.2 (mediate-thick curve). This developed fabric anisotropy is rather weak but reasonable because the experienced deviatoric strain is only about ε d = ε 1 À ε 2 = 0.014, which is not sufficient to produce significant fabric changes. After a simple shear strain of γ = 0.3, the fabric anisotropy mobilises at a c = 0.18 and its resulting principal direction is θ c = 41.5°, which is a significant change from the initial state of simple shear.
The evolution of the fabric anisotropy and its principal orientation along simple shearing are shown in Figure 9. Instead of repeatedly fitting them from the measurements using Eq. 6a, it has been shown that the parameters a c and θ c can be more conveniently calculated via the related fabric tensor F ij in the 2D case [55,56]: where i and j designate the components in the reference frame, and N c is the total number of contacts in the cell volume V. The θ c is given by the principal orientation of the tensor F ij , and the anisotropy parameter is given by where F 1 and F 2 are the principal values of the fabric tensor F ij . Figure 9a shows that, before simple shear, the degrees of fabric anisotropy are a c = 0.019, 0.037 and 0.080 for the derivative samples from specimens S1, S2 and S3, respectively. This suggests higher fabric anisotropy develops in a sample with a higher initial porosity during the anisotropic compression, as a result of the larger magnitude of deviatoric strain ε d taken to reach the target stress ratio.
During the early stage of simple shear, the fabric anisotropy increases faster in a denser sample, though its initial degree of fabric anisotropy is smaller than that in a looser sample. Similar to the development of mobilised shear stress shown in Figure 5, the fabric anisotropy in the dense sample first reaches a peak then decreases to the ultimate value, whereas for loose samples, the fabric anisotropy gradually increases and then fluctuates around the ultimate value (Figure 9a). The evolution of principal fabric orientation θ c exhibits close similarity to that of the principal stress orientation. The curves start from about 0°and gradually approach the principal strain rate direction at around 45°, and a looser sample exhibits a slower approach rate (Figure 9b).

Effect of initial stress ratio η 0
In this section, we examine the behaviour of simple shear in samples under the same initial mean stress but different initial stress ratios. The derivative samples of specimen S2 with a mean stress p 0 = 200 kPa were studied.
We shall first examine the resultant η À p loading history. The isotropic and anisotropic compressions are of pure stress-loading paths; therefore their loading histories are readily prescribed as have been given in Figure 3. Simple shear mode is a mixed-loading path so that its η À p loading history is a resultant response. As shown in Figure 10a, in all four cases with η 0 ranging from 0 to 0.3, the stress ratio increases quickly from its initial value to the ultimate value around 0.3. The mean stress p develops to about equivalent to the servo-controlled vertical stress σ y (Figure 10a), which effectively diminishes the difference between the two normal stress components σ x~σy (Figure 10b). Figure 10b compares the evolution of stress components σ xy and σ x against the shear strain γ. The rate in mobilising the coulomb shear stress ratio, that is, σ xy /σ y , is slower if the initial stress ratio η 0 is higher.  Figure 11 shows that a larger initial stress ratio η 0 results in a slower rate in approaching coaxiality. This effect can be readily explained by examining the governing equation of the principal stress orientation Eq. 5a: a larger initial difference between σ x and σ y corresponds to a smaller θ σ .
We observe a more contractive behaviour in a sample with a larger initial stress ratio ( Figure 12). It has been generally shown that a larger initial porosity usually corresponds to a more contractive behaviour under shearing, for example, in biaxial shearing. Such observation appears not always true in simple shear. We note that the initial porosity of the simple shear sample (i.e. the resultant porosity from anisotropic consolidation) is actually lower at a higher stress ratio in the range η 0 = 0 to 0.2. This suggests that other fabric properties can be more important in affecting the simple shear volumetric behaviour than the porosity that is a scalar measure. Figure 12 also shows that the stress ratios η in four samples all increase quickly from their respective initial values to the plateau. After shear strain γ = 0.05, the stress ratios are indistinguishable among the four samples.
The influence of the initial stress ratio on the fabric evolution is shown in Figure 13. As a result of previous anisotropic compression, the fabric anisotropy a c before simple shear is larger in a sample with a larger initial stress ratio (Figure 13a). In samples with lower initial stress ratios (η 0 = 0, 0.1 and 0.2), anisotropy parameter a c increases quickly in the range γ < 0.1. For the sample with the highest initial stress ratio (η 0 = 0.3), a c does not increase until shear strain passes γ = 0.05, after which the curves of the four cases largely merge. The delayed development of fabric anisotropy in the more anisotropic sample reflects that the influence of the previous loading history needs to be largely erased before the fabric pattern conforms to the simple shear mode. The evolution of the principal fabric orientation (Figure 13b) further confirms this effect: in samples with smaller initial stress ratios (η 0 ≤ 0.3), the principal fabric orientation quickly rotates to about 45°, whereas it takes much longer for a sample with the largest initial stress ratio (η 0 = 0.3). It is worth noting that the initial principal fabric orientation θ c deviates significantly away from 0°for samples with small initial stress ratios (Figure 13b). This is because the deviatoric strains taken to reach the stress ratios are not sufficient to develop the fabric orientation to the loading direction.

Effect of initial lateral pressure ratio K 0
The simple shear samples presented in previous sections are all sheared with the servo-controlled normal stress σ y larger than σ x . This scenario corresponds to a laboratory test scenario with an initial lateral Figure 11. Effect of initial stress ratio η 0 on the evolution of non-coaxiality during simple shear. pressure ratio K 0 = σ x /σ y < 1. The samples can also be sheared with K 0 > 1, that is, with σ y < σ x . In our numerical simulations, this can be simply achieved by rotating the reference frame by 90°. The development of stress ratios (η) and dilations (ε v ) along simple shearing in the two different cases of K 0 are compared in Figure 14. For cases with initial stress ratio η 0 = 0.2, the sample hardly experiences any contraction before dilating in the case of K 0 > 1, and the overall dilation is much larger than that in the case K 0 < 1. As known from the previous section, the magnitude of the normal stress component in the simple shear direction will eventually approach that of the confining stress. Hence, the larger dilation in the case of K 0 > 1 than in the case of K 0 < 1 may be attributed to the lower final mean stress in the former than in the latter. For samples with initially isotropic stress, the curves from cases K 0 > 1 and K 0 < 1 are globally indistinguishable. The influence of K 0 on the development of stress ratio η is rather small, for both cases of η 0 = 0 and 0.2 ( Figure 14). Figure 15 shows that K 0 has a noticeable effect on the evolution of non-coaxiality. For the samples with initially isotropic stress (η 0 = 0), the non-coaxiality θ ε À θ σ is initially zero but deviates to a small non-zero value (thin solid-line and thick dashline in Figure 15a). The emergence of the slight noncoaxiality for initially isotropic sample is due to the sample dilation under shearing that deviates the principal strain rate orientation away from 45°. For the samples with η 0 = 0.3, the case of K 0 > 1 exhibits a faster rate in approaching coaxiality than the case of K 0 < 1. The difference should be  Figure 14. Evolution of stress ratio and volumetric strain during simple shear with initial lateral pressure ratio K 0 < 1 and K 0 > 1. Figure 15. Evolution of (a) non-coaxiality and (b) stress and strain rate directions during simple shear with initial lateral pressure ratio K 0 < 1 and K 0 > 1.
caused by different extents of dilation in the two cases. As also shown in Figure 15b, the direction of strain rate in the case of K 0 < 1 is initially smaller than 45°due to contraction, then increases to be larger than 45°due to dilation. In contrast, the strain rate direction in the case of K 0 > 1 is larger than 45°from nearly the beginning of the shearing due to much larger potential of dilation. Consequently, the variation in strain rate direction leads to the variation in rate in approaching coaxiality shown in Figure 15a. Figure 16 shows that the case of K 0 > 1 produces a higher degree of fabric anisotropy a c than the case of K 0 < 1. The global evolutions of principal fabric orientation θ c are generally not distinguishable within the discreteness of the samples, which is similar to that in the stress ratio developments shown in Figure 14.
3.5. Effect of initial mean stress p 0 As our systems are composed of relatively 'soft' particles, the influence of particle elasticity may be significant when the normalised confining stress is relatively large. Consequently, the role of mean stress on the simple shear behaviour may not be negligible, which is therefore explored in this section. This is also equivalent to testing the effect of particle stiffness.
We compare the results of samples with initial mean stress ranging from p 0 = 50 to 800 kPa in Figures 17 and 18. We observe that a larger mean stress results in a significantly slower rate in the development of normal stress σ x and a slightly slower rate in that of shear stress σ xy (Figure 17a). This dependence of shear stress mobilisation on the mean stress may be attributed to that of the potential of dilation as explained by Bolton [53]. Indeed Figure 17b shows that the sample is more contractive when with larger mean stress. Figure 17c suggests a rather strong effect of the mean stress on the evolution of non-coaxiality: larger mean stress leads to slower approach to coaxiality. Figure 18 shows that the fabric anisotropy parameter a c during simple shear is smaller if the mean stress is larger. The underlying mechanism may be attributed to the reduced dilation and increased degree of connectivity associated with the larger mean stress. The effect on the principal fabric orientation θ c is, however, insignificant.
Based on the aforementioned evidence, in some aspects, a sample with a larger mean stress appears to behave like a 'looser' sample, which is consistent with the observations in biaxial/triaxial tests. Nevertheless, the trend in the evolution of the stress ratio η appears to be an exception in which larger mean stress leads to larger peak strength (see Figure 17b in contrast to Figure 5b). This requires further investigation. Figure 16. Evolution of fabric anisotropy a c during simple shear with initial lateral pressure ratio K 0 < 1 and K 0 > 1 (specimen S2, p 0 = 200 kPa).

CONCLUSIONS
We have presented a DEM study of the quasi-static non-steady simple shear flows of a 2D assembly of cohesionless particles. We abandoned the rigid plane boundary configuration that had been popularly adopted in many earlier experimental and numerical studies, by developing a discretised-wall confined granular cell. In such configuration, the cell boundary is discrete in nature and each segment moves strictly conforming to the prescribed strain rate, allowing synchronised dilations between the boundary and the confined solid. Sufficiently uniform distribution of the stress-strain across the whole assembly has been achieved, which produces more reliable data for rheological studies to generalise constitutive models for continuum methods dealing with much larger systems. We have examined two aspects of the simple shear behaviour: macroscopic stress and strain rate evolution, particularly the non-coaxiality between the principal directions of the two, and micromechanics such as evolution of fabric. For an initially anisotropic specimen sheared under constant normal pressure condition, the direction of principal stress rotates towards that of the principal strain rate, gradually reducing the degree of non-coaxiality from about 45°to fluctuating around 0°. The rate in approaching coaxiality is slower in samples with larger initial porosity, stress ratio and mean stress. Generally, a faster approach to coaxiality in simple shear was observed in a more dilatant sample which often shows a larger degree of mobilised fabric anisotropy, suggesting the potentially important role of instantaneous internal friction angle. The evolution of principal fabric direction resembles that of the principal stress direction.