Prediction of Shear Localisation in Granular Materials based on a Critical State Non-coaxial Model

Experimental evidence indicates that the shear localisation acts as a precursor to the failure in biaxial compression tests of granular materials. Once formed they are persistent and lead to progressive failure of most geotechnical structures. It is generally accepted that the primary mode of deformation within these shear bands is simple shear which is accompanied by rotation of principal axes. Hence, the conventional plasticity theories based on the assumption of coaxility is not sufficient to describe the behaviour within those shear bands. This paper highlights the influence of the non-coaxility on the initiation and orientation of shear bands in both drained and undrained sand. The con-coaxial plasticity theory is integrated into a critical state constitutive model enriched with the state parameter concept. The model is capable of taking account of the variation of lode angle under plane strain condition. Numerical plane strain biaxial compression tests are conducted to observe the effect of non-coaxility on shear localisation. Bifurcation criteria based on the acoustic tensor are checked to predict the onset and inclination of the shear band. Predictions from the non-coaxial model are compared with those of coaxial model. The influence of the initial void ratio for the formation of shear bands is explored. Results are compared qualitatively with experimental


Introduction
The mathematical modelling of the strain localisation consists of two main phases: detection of the onset of material instability and capturing the post localisation. Both of which are important to get a clear view of the failure of geo-materials. The former identifies the emergence and inclination of the localised band, whereas the latter determines a physically meaningful representation of the material during a localised deformation.
The objective of this study is to evaluate the role of the constitutive model to capture the onset of shear localisation in granular material. The influence of non-coaxial flow rule for the detection of localisation is explored. The influence of the initial state on the formation of shear band is explored for both drained and undrained conditions.

Literature Review
According to the literature, the prediction of onset and inclination of the shear band is governed by pre localised constitutive relation. It is reported that classical plastic idealisations often poorly predict the onset of shear band [1][2][3][4][5][6]. Smooth yield surfaces and normality condition render a symmetric elasto-plastic stiffness matrix. In this case, the criteria for failure, the diffused instability and the localisation coincide with the peak stress ratio. This is not * e-mail: hem42@cam.ac.uk consistent with experimental findings which state that instability and/or localisation in a single phase material tend to appear in the hardening regime. Therefore, both nonnormality of plastic strain increments and vertex like yield surfaces have been adopted in the past to induce the destabilising effect [2]. According to [1], the lack of symmetry in elasto-plastic stiffness itself is sufficient to the loss of material stability. [2], [3], [4] reported that even with non-associative flow rule, there is a significant difference between theoretical and experimental predictions of onset of bifurcation. They attributed this observation to a major deficiency in the ordinary flow rule which is built on the assumption of coaxility between the direction of principal plastic deformation rate and direction of principal stress. This assumption leads to an overestimation of stress increment in the direction tangential to the yield surface.
As a solution, [2] introduced the yield vertex theory in order to bring a destabilising effect to the smooth yield surface. It is reported that the vertex effect significantly alters the predictions of localisation of models that entails smooth yield surfaces. Following the same footsteps, [4] extended the critical state Cam-Clay model (with associative flow rule) to incorporate a non-coaxial term. They reported that the non-coaxial term facilitates an easier access to the elliptic/hyperbolic boundary. Further, they state that non-coaxial models are more inclined to the instability by localisation than coaxial ones even without any kinematic constraints. This is mainly due to the decrease of the in-stantaneous shear moduli of the shear stress. [4] further state that the existence of microscopic localisation bands may not result from the coaxility assumption. [3] also discovered that the shear band orientation and its inception are influenced by the degree of non-coaxility according to the yield vertex theory. Further, they report that the strong non-coaxility at low confining pressure leads to the diffused bifurcation whereas decreased noncoaxility at higher confining pressures favours the shear localisation.
Tangential plasticity concept by [7] follows fundamentally the same notion. It states that the direction of the plastic strain rate depends on the direction of stress rates. He reports that the tangential stress rate (caused by noncoaxial plastic strain) increases the shear band inclination angle. [8] mentions that the excessive plastic volume expansion stems from the associative plasticity can be circumvented by incorporating the tangential plasticity alone.
Recent bifurcation studies on the drained and undrained sand in plane strain condition were conducted by [5,6]. They developed a non-associative and noncoaxial Mohr-Coulomb model. They report that the stress strain relationships are hardly changed, but the prediction of localisations (singularity of the determinant of acoustic tensor) is significantly improved from the non-coaxial theory. However, they further state that the non-coaxial flow rule does not alter the prediction of diffuse instability (zero second order work). Their non-coaxial predictions of shear localisations are consistent with experimental results of both dry and saturated undrained sand whereas coaxial models tend to underpredict the strain at the initiation of the shear band. It should be noted that in the medium dense sand, the localisation tends to occur in the hardening regime (before the peak effective stress ratio) although the deviatoric stress is decreasing.
[9] also adopted a non-associative, non-coaxial flow rule to predict the shear band formation in the undrained dense sand. He reported that the non-coaxial bifurcation point is close to experimentally observed localisation. However, in this case, the shear band appears after the peak stress ratio when the deviatoric stress is still hardening. In summary, above bifurcation studies agree that noncoaxial models facilitate localisation compared to coaxial models.

Non-coaxial plasticity
In conventional elasto-plastic theory, the total strain rate d i j is composed of plastic d p i j and elastic d e i j components (Equation 1a). The traditional plastic flow theory assumes that the plastic flow depends only on the current state of stress σ i j : the direction of the plastic flow is normal to the plastic potential function g as given in Equation 1b where dλ is scalar plastic multiplier.
Representation of stress path and plastic strain increment path in deviatoric plane, principal stress direction = α, principal plastic strain increment direction = α , angle of non-coaxility = ψ, [10] This assumption fails to describe the phenomenon observed in stress probes [10][11][12], hollow cylinder [13] or simple shear tests [14]. These laboratory experiments suggest that the magnitude and direction of plastic deformation are dependent on the direction of stress rates. Hence the principal stress direction does not coincide with the principal plastic strain increment direction during shear deformations or non proportional loading. This phenomenon is displayed in Figure 1. The degree of noncoaxility is greater in the beginning, but principal stress and principal plastic strain rate directions coincide near the critical state [10][11][12]14].
Yield vertex theory modifies the plastic strain rate in Equation 1a to include both coaxial d pc i j and non-coaxial d pn i j parts.
The coaxial part d pc i j can be calculated by the conventional plastic flow theory in Equation 1b whereas the non-coaxial plastic flow rate d pn i j is mathematically determined by equation 3a and 3b in which h nc is the noncoaxial plastic modulus, G is shear modulus (G = Ap r ) and n i j represents the vector of non-coaxial plastic flow direction. N i jkl in Equation 3b denotes the flow direction matrix.
Original yield vertex theory by [2] suggests that noncoaxial part n i j is normal to the unit vector along the deviatoric stress (Equation 4 and Figure 2). This has been adopted for simple constitutive models with circular yield surfaces such as Drucker-Prager and Mohr-Coulomb with only isotropic hardening.
where s i j = σ i j − σ kk 3 is the deviatoric stress tensor and τ = 0.5s i j s i j .
In this paper, a generalised yield vertex theory (Equation 5 and Figure 3) developed by [7] is utilised. This is  Illustration of coaxial and non-coaxial plastic strain rate for non-circular yield surface f assuming the associative plasticity [7] valid for non-circular yield surfaces f with any type of hardening. Equation 5 coincides with the original yield vertex theory when the yield surface is circular in the deviatoric plane and only isotropic hardening is adopted.
The final constitutive relationship is given by where D e and D e p are the elastic and elasto-plastic stiffness matrices respectively. The initial two terms indicate the coaxial effect and the last term deflects the non-coaxial effect. A smaller non-coaxial modulus leads to a higher non-coaxial influence.

Nor-Sand (NS) model
Aforementioned yield vertex theory can be utilised with any constitutive model. In this paper, the Nor-Sand model (NS) by [16] based on fundamental axioms of the critical state theory has been used. It is an elasto-plastic model which incorporates associative flow rule with a bounding surface concept. The state parameter concept integrates the effect of both void ratio and the mean pressure. Hence a single set of parameters is required for a given sand. NS yield surface is given by where p, q and M θ are three stress invariants namely effective mean pressure, deviatoric stress and critical stress ratio. M θ is a function of lode angle θ and calculated by Matsuoka-Nakai failure criterion [17]. N is the dilatancy parameter which specifies the shape of the yield surface. The size of the bounding yield surface is determined through the image pressure p i . At image state the dilation is zero, but the rate of dilation is non-zero. Hence this model can decouple two criteria of the critical state. NS assumes an isotropic hardening rule driven solely by the plastic deviatoric strain p q .
where h is the hardening or softening modulus and p i,max is the maximum image pressure which is a function of state parameter at the image state Ψ i .
The critical void ratio at image state (e c ) i is defined by the relationship where e max and e min are maximum and minimum void ratios respectively. The hardening of the yield surface is restricted to match the maximum dilatancy of the real sand. According to experimental data, a linear relationship can be obtained between the maximum dilatancy D max and Ψ i .
where χ tc and M tc are the maximum dilatancy coefficient and critical stress ratio obtained from triaxial compression data. This maximum dilatancy is used to size the yield surface so that the real soil behaviour matches with the dilatancy from the normality. Hence, the limiting hardness or the maximum yield surface size can be calculated from Due to this limiting hardness concept, the NS is capable of properly capturing the soil behaviour in dry side of the critical state even with an associative flow rule. In line with [18]'s concept of sub-loading surface, the current NS model is elasto-plastic from the beginning. Further information about the implementation of the NS model can be found in [17].

Bifurcation analysis
The bifurcation criterion states that the determinant of the acoustic tensor to be zero at the onset of localisation. In other words, the singularity of the acoustic tensor is the necessary condition for the loss of ellipticity of governing partial differential equations.
where the acoustic tensor is defined as where n j is unit normal to the shear band. The derivation of this criterion is given below. The direction of the shear band θ sb is determined during the process of calculating the minima of det [Q]. A change in sign of the minima of the polynomial f (θ sb ) signals the onset of localisation. In two dimensional plane strain condition (x-y plane) it can be calculated as = a 4 tan 4 θ sb + a 3 tan 3 θ sb + a 2 tan 2 θ sb + a 1 tan θ sb + a 0 (16)

Biaxial compression tests
A series of single element drained and undrained biaxial compression tests are conducted to detect the bifurcation based on the constitutive behaviour. A plane strain fournode element with unit height and width is subjected to the biaxial compression analysis. The bottom left node is constrained in two directions whereas the bottom right node is assumed to be roller supported. The top nodes are unconstrained. During the initial step, all boundary constraints and initial void ratio are specified. Additionally, a constant volume constraint is applied only for the undrained analysis. A confining pressure of 100 kPa is applied to the side and top surfaces during the second step. A vertical downward displacement of 0.5 m is applied to the top nodes at the last step.
Constitutive formulations are implemented in ABAQUS finite element software as a user defined material model. NS material parameters are given in Table 1. Three initial void ratios are considered to explore the behaviour of dense (e-0.55), medium dense (e-0.75) and loose (e-0.85) sand with both drained and undrained boundary conditions. Both onsets of the bifurcation and inclination angles are calculated based on the criterion of the acoustic tensor. In the undrained tests, the acoustic tensor is calculated based on the effective stress.

Results and Discussion
The shear localisation is only a subset of instabilities observed in geomaterials. Although the emergence of a shear band is predicted based on the material behaviour itself, other external parameters also significantly affect it. Apart from initial states (void ratio and confining pressure), the load path (boundary and drainage conditions) as well as controlling variable (load or displacement) govern the initiation and the type of instability. The focus of this study is only about the loss of ellipticity or shear localisation signalled by the singularity of the acoustic tensor. Diffused instability and liquefaction conditions are not considered.
Within the general instability study, two limiting stress states are important in addition to the bifurcation point: Tresca limit state at which the deviatoric stress is maximum and the Coulomb limit state at which the effective stress ratio is maximum. Depending on the type of constitutive model and the flow rule adopted these limiting states may occur simultaneously with the bifurcation point. In traditional geomechanics, it is conventional to take either of these limiting states as the failure criterion. Table 2 and 3 present basic results of bifurcation analysis: deviatoric strain at the onset of localisation and orientation of the shear bands (degrees) to the horizontal. Tresca and Coulomb limit states are also included for the comparison. Results are discussed qualitatively with experimental observations.
The non-coaxility is triggered by the principal axis rotation. Hence the mechanical response of the soil in the plane strain condition is not much affected as reported by [6]. However, slight differences can be seen near the peak due to the rotation of intermediate principal stress. The stress strain relationships are only shown for dense drained and undrained condition in Figure 4, due to space limitation. The influence of varying lode angle during plane strain deformation activates the non-coaxility. Therefore, predictions on the emergence and orientation of shear bands are different between coaxial and noncoaxial models.
For the drained analysis, the coaxial NS model predicts localisation only for dense sand. According to the non-coaxial NS model, shear bands emerge in drained sand irrespective of the initial density. This is in line with [19,20]'s experimental predictions which state that shear bands can occur in both loose and dense sand in the biaxial displacement controlled tests.
As for drained dense sand, the non-coaxility delays the localisation and increases its inclination. Further, the  shear band occurs in the hardening regime of the mobilised friction (before maximum stress ratio). These are consistent with the findings of [19]. On the other hand, in the coaxial NS model, the bifurcation point coincides with the Coulomb and Tresca states. This is because of the symmetric stiffness matrix resulted by the associative flow rule. It is also evident that under the coaxial flow rule, the localisation is possible only for the soil with negative state parameters. The initial void ratio has a direct impact on both the emergence and inclination of bands. The initiation of localisation is delayed by the increase of void ratio. The shear band angles to the horizontal axis tend to decline as the void ratio increases.
Both coaxial and non-coaxial NS models do not predict bifurcation for the undrained dense sand. On the contrary, undrained loose sand satisfies the bifurcation criterion at a very early stage of the deformation. However, unlike dry sand the singularity of acoustic tensor based on effective stress does not necessarily indicate the localised instability in saturated sand. According to [21] for the formation of shear band both drained and undrained acoustic tensor should be zero. This study only calculated the initial bifurcation condition.
The non-coaxial flow rule leads to a delayed bifurcation and increased inclination angles. It can be concluded that the irrespective of the flow rule, bifurcation under isochoric constraint is possible only for soil with positive initial state parameters.
When the drainage is suppressed, the pore volume must theoretically remain constant since the fluid is considered as incompressible. Thus the granular skeleton is constrained to an isochoric deformation. When the material deforms without a volume change, its tendency to dilate or contract induces negative and positive pore pressure respectively. In the case of the loose specimen where the positive excess pore pressure is generated, the effective stress can drop down to zero leading to a liquefaction. On the contrary, in the undrained dense granular material pore pressure is reduced and the effective stress is increased. This delays the localisation or any other kind of instabilities and it is popularly known as dilative hardening.
However, there is no proper consensus regarding the experimental findings of localisation in the undrained sand. [22] claimed that shear bands are observed in undrained loose sand whereas [23] denied it. [22] concluded that localisation in the saturated undrained dense sand is precluded unless the isochoric constraint is relaxed by cavitation. On the contrary, [23] reported localised bands near the Coulomb state of the undrained dense sand. Further [23] reported that in the medium dense sand, the shear band was observed during the hardening regime of effective stress ratio after the Tresca state. Predictions of this study are consistent with [22] for dense and loose sand and with [23] for medium dense sand.
The bifurcation of saturated sand is much influenced by the soil and pore fluid coupling. Therefore, in laboratory experiments, the rate of pore fluid dissipation inside the sample influences the overall stability. Although globally undrained, the internal water movements make the bifurcation problem subjective of time and scale. Therefore, it is difficult to assess the potency of rate independent constitutive relations to capture the shear localisation in the undrained soil.
It is mentioned that the non-coaxial model delays the bifurcation in both drained and undrained sand. This delay is sensitive to the magnitude of the tangential modulus. The lower its value, the shear band is detained longer.

Conclusion
• For the drained condition, the coaxial NS model predicts bifurcation only for soil with negative initial state parameter, but non-coaxial NS model detects bifurcation regardless of the initial state of soil • For the undrained soil with initial negative state parameters, both models do not predict any instability • Both coaxial and non-coaxial models detect bifurcation for undrained soil with positive state parameters.
• Non-coaxiality delays the onset of localisation and increases its orientation • Non-coaxial model always predicts the bifurcation before the maximum effective stress ratio (in the hardening regime) • Bifurcation in the saturated soil is governed by the soil and pore fluid coupling. Hence, the local drainage condition can differ from the global boundary condition. Therefore, experimentally observed localisation in globally undrained dense sand cannot be captured by constitutive models