Bifurcation analysis of shear band in sand under true triaxial conditions with hypoplasticity

Predicting the onset of shear band is of significance in understanding the failure of geomaterials. The prediction accuracy is dictated by the constitutive model used for the description of the pre‐bifurcation behaviour. In this study, we first modify a recently proposed hypoplastic constitutive model by incorporating a general strength criterion and a stiffness function. We proceed to consider the onset of shear band in sand under true triaxial conditions. We demonstrate that our analyses capture the pre‐bifurcation stress–strain relationship at different values of intermediate principal stress and predict fairly well the onset of shear band. The acoustic tensor criterion generally adopted in elastoplastic approaches is inadequate for hypoplastic approaches. No special non‐coaxial treatment is required for the present approach to yield a reasonable variation trend of bifurcation strain with intermediate principal stress ratio b .

investigated, including boundary condition, sample size, the sand's density and particle size, confining pressure, imperfection and heterogeneity, inherent and stress-induced anisotropy, and so forth.
Accompanied by the experimental studies on the shear band, some theoretical studies have been conducted to detect the onset of shear band, among which the most notable is that by Rudnicki and Rice. 17 The onset of shear band is mathematically treated as a bifurcation problem, that is, the constitutive relations may allow the homogeneous deformation pattern to lead to a bifurcation point, at which an inhomogeneous deformation pattern with the jump of strain rate in a planar band can be admissible. The singularity of the determinant of acoustic tensor, which is associated with the elastoplastic tangential modulus and the inclination of shear band, denotes the onset of shear band. With this bifurcation criterion, the prediction of the onset of shear band is greatly dependent on the pre-bifurcation constitutive model. Some studies have been carried out by using different elastoplastic constitutive models. It has been found that the constitutive models based on classical plasticity theory often provide poor predictions of the onset of shear band. Therefore, efforts have been made to improve the prediction capability by incorporating more versatile elastoplastic constitutive models, including nonnormality, noncoaxial, sub-loading surface model and other extensions, 1,[18][19][20][21][22][23][24][25][26][27] of which the noncoaxial plastic models may be the most successful. [20][21][22][23]27 Hypoplastic constitutive models, as an alternative to elastoplastic models, have aroused increasing interests recently. [28][29][30][31][32][33][34][35][36][37][38][39][40][41] Some beneficial features of hypoplastic constitutive models include simple formulation, fewer parameters, non-linear behaviour before failure and smooth transition from contractancy to dilatancy. In terms of the bifurcation analysis of shear band, the studies based on hypoplastic models are rare in comparison to those based on elastoplastic models. As hypoplastic models are incremental non-linear, the bifurcation criterion should be different. [42][43][44] Tamagnini et al. 42 presented a critical review of results obtained based on two particular classes of hypoplastic constitutive models, which are the most popular ones at that time. Wu 45 proposed a bifurcation criterion for non-linear analysis with no additional assumptions. The criterion provides a lower bound solution and is more general than the incremental linear ones. Huang 46 proposed general bifurcation conditions for a non-polar hypoplastic model and its micro-polar continuum extension. Huang 47 investigated the possibility of localized failure in true triaxial tests by examining the bifurcation condition with hypoplasticity in element tests following different stress paths. As mentioned above, the predicting accuracy of the onset of shear band is greatly determined by the pre-bifurcation constitutive model. With the development of hypoplastic constitutive models, the older ones have been gradually replaced by novel ones (e.g. Refs. [28,48]). However, the predicting capability of the onset of shear band of these novel models remains unclear.
This study aims at predicting the onset of shear band in sand under true triaxial conditions within the framework of a recently proposed hypoplastic basic model for sand. 48 In order to well capture the pre-bifurcation stress-strain behaviour of sand under true triaxial conditions, the basic model is extended to describe the sand behaviour under true triaxial conditions by incorporating a general strength criterion proposed by Yao et al. 49 Furthermore, the constitutive formulation is revised by introducing a stiffness function to better fit the tangential stiffness of sand under true triaxial conditions. With the novel hypoplastic model constructed, the bifurcation criterion proposed by Wu 45 is reviewed and further discussed. Emphasis is put on the comparison with the widely accepted acoustic tensor criterion proposed by Rudnicki and Rice. 17 The model responses, in the aspects of both describing the homogeneous mechanical behaviour and predicting the onset of shear band are demonstrated. Consequently, the numerical results by the proposed approach are compared with the laboratory observations of sand under true triaxial conditions.
Notations and conventions: Compact or index tensorial notation is used throughout. Vectors and second-order tensors are denoted with bold letters (e.g. , ) and fourth-order tensors are denoted with calligraphic bold letters (e.g. ). Symbols '⋅' and '∶' between tensors of various orders denote inner product with single and double contractions, respectively. The dyadic product of two tensors is indicted by '⊗'. ‖̇‖ = √ tr(̇2) denotes the Euclidean norm and the arrow operator is defined as ⃗̇=̇∕ ‖̇‖. Following the convention of soil mechanics, compression stress and contraction strain are taken to be positive.

HYPOPLASTIC CONSTITUTIVE MODEL
In a hypoplastic constitutive model, the stress rate in response to material deformation is expressed as a non-linear tensorvalued function of Cauchy stress and strain rate. A general hypoplastic framework for describing rate-independent material behaviour can be written as 39,50̊= where the terms  and , respectively, denote the linear and non-linear components, anḋ, respectively, stand for the Cauchy stress tensor and the stretching strain rate tensor. The Jaumann stress rate is defined bẙ where denotes the spin tensor. The stretching tensor and spin tensor are related to the velocity gradient tensor througḣ where is the velocity and ∇ is the gradient operator. Within the above framework, a novel basic hypoplastic constitutive model for sand has been recently proposed by Wu et al. 48 In the novel basic model, there are three linear terms and one non-linear term: = 1 (tr )̇+ 2 (tṙ) + 3 tr(̇) tr in which ( = 1, 2, 3, 4) are dimensionless parameters. * denotes the deviatoric stress tensor defined by * = − (tr ) ∕3 with being the Kronecker delta tensor. The above basic model has a failure surface of a conical shape, which is similar to that of the Drucker-Prager model. In order to describe the soil behaviour under true triaxial conditions, a general strength criterion proposed by Yao et al. 49 is incorporated in this study. According to the criterion, noting that the yield locus on the deviatoric plane of various materials are generally inside the Mises circle and outside the spatially mobilized plane (SMP), 51 the radius of the yield locus for a specific Lode angle can be taken to be the weighted average of those of the Mises criterion and the SMP criterion. Following this general strength criterion, the ratio between the radius of failure surface on the deviatoric plane for a specific Lode angle and the radius under triaxial compression state can be calculated by, in which is a weighted factor lying between 0 and 1, with = 0 corresponding to the Mises criterion and = 1 corresponding to the SMP criterion. By adjusting the weight value between 0 and 1, a continuous transition from the Mises criterion to SMP criterion can be obtained. Given an arbitrary stress state, and are the corresponding length of failure surface on the deviatoric plane under triaxial compression stress state for the Mises criterion and the SMP criterion, respectively. and can be calculated by Ref. [49], where 1 , 2 and 3 are the stress invariants. Note that the specific formula is slightly different from that in Ref. [49], as the so-called transformed stress space is used in their work. In order to obtain a failure surface with the above strength criterion in hypoplasticity, let us consider the spherical and deviatoric parts of the constitutive equation (4), respectively. The spherical part is tr̊= 1 (tr )(tṙ) + 2 (tr )(tṙ) Note that in the fourth term tr * vanishes since tr * = 0. The deviatoric part is * = 1 (tr )̇ * + 2 * (tṙ) + 3 tr(̇) tr * + 2 4 * ‖̇‖.
Noting that tṙ= 0 at critical state, by letting tr̊= 0, Equation (7) can be rewritten as Noting that tṙ= 0 and tr(̇) = tr( * ̇) at critical state, by letting̊ * = 0, Equation (8) can be rewritten as According to Equations (9) and (10), noting that tr in the deviatoric plane remains a constant, the deviatoric stress at failure surface for a specific Lode angle can be adjusted by dividing the fourth term by and dividing the third terms by 2 . With this simple approach, * of the failure is multiplied by which is associated with the Lode angle, while tr remains unchanged. The radius on the deviatoric plane of failure surface is thus adjusted to fulfil the general strength criterion. Therefore, the extension of the basic hypoplastic model considering the general strength criterion can be written as, It has been reported that in some cases the hypoplastic models predict improperly the development of stiffness with shear strain. A common case is that the model predicts a sharp bend of the stress-strain curve before the critical state is achieved, whereas the experimental results show a more gradual stiffness decrease. 28,40 To overcome this shortcoming, Wu and Kolymbus 40 suggested multiplying a stiffness function in the model. In this study, the stiffness function is modified to better fit the true triaxial case, which is calculated as where denotes the stress ratio ‖ * ‖∕(tr ), is the ratio defined in Equation (5) and is a parameter to reflect the effect degree of the ratio . Note that multiplying the whole constitutive by a scaling function does not affect the failure surface. A more sophisticated constitutive model can be constructed by considering the evolution of void ratio, which can be straightforward. 48 However, the authors select to keep the constitutive model simple because insufficient experimental data on shear band are accessible. Consequently, the refined constitutive model proposed in this study can be written in the following complete form:̊= There are totally six parameters in the model. The four parameters 1 ∼ 4 are the same with the original model, which can be determined using a standard procedure described by Wu and Bauer. 39 The additional two parameters and are related to strength and stiffness, respectively, which offer more flexibility to fit model prediction to experimental data under true triaxial conditions. It should be pointed out that the rate dependence 52,53 and anisotropy [54][55][56][57][58] is out of the scope of this study.
The numerical integration of the hypoplastic constitutive equation (13) is more straightforward when compared with the elastoplastic model, due to the fact that the hypoplastic model only has a single equation. 35 Since hypoplastic model does not have a yield surface, the stress return mapping algorithms common for elastoplastic models [59][60][61] are not needed. In this study, a simple explicit forward Euler scheme with constant step sizes 62,63 is adopted. This scheme can achieve satisfactory accuracy as long as the step size is small enough, which can be determined by a trial and error approach. Since only one stress point is considered here, the time cost for numerical integration is low. More complicated numerical integration schemes for hypoplastic model are available in Ref. [35]. Wu 45 proposed a bifurcation criterion for hypoplastic materials. In this section, the bifurcation criterion is briefly reviewed. Then some further discussions are made. Emphasis is put on the comparison with the acoustic criterion, which has been widely accepted as a general bifurcation criterion for elastoplastic materials (e.g. Refs. [1, 20-23, 25, 27]).

F I G U R E 1 Schematic interpretation of shear band
The schematic interpretation of shear band is shown in Figure 1. The development of a shear band bifurcation in a homogeneously deforming material can be completely defined with the following kinematic conditionṡ in which the superscripts +' and −' refer to values inside and outside the shear band, respectively, the superscript } ' denotes the jump of the value of the shear band, is the unit normal vector to the shear band and is a vector defining the velocity gradient jump in the direction . In addition, the requirement of continuing equilibrium at the onset of shear band implies the following static condition: in whicḣ+ anḋ− are the Cauchy stress rates inside and outside the band, respectively. Note that the stress itself is continuous before bifurcation, and thus the jump of the Jaumann stress rate can be obtained according to Equation (2) as̊=̊+ Meanwhile, the stress rates inside and outside shear band can be obtained from the constitutive Equation (1) F I G U R E 2 Geometrical meaning of the jump of the strain rate across shear band The basic equation for bifurcation analysis can be obtained by incorporating Equations (16) and (17) into Equation (15), in which is the acoustic tensor calculated by = , where is the component of  in Equation (1); is a second-order tensor associated with finite strain and is calculated by while is calculated by By introducing the orientation angle and inclination angle of shear band, as defined in Figure 1, the components of the unit normal vector of the shear band can be calculated as follows: 1 = cos cos , 2 = cos sin , where and are, respectively, the orientation and inclination angles of shear band .
The onset of shear band means that a solution other than ≡ 0 exists for Equation (18). To resolve the problem we need to have a closer look at the jump of strain rate across the shear band. As shown in Figure 2, the jump of strain rate tensor must fulfil the triangle inequality, which reads, where is termed as the degree of bifurcation. The bifurcation problem can be solved with the regula-falsi method. Specifically, the unit normal vector (i.e. the orientation and inclination) of shear band is first assumed and then the degree of bifurcation d can be obtained according to Equation (18) and Equation (22). The maximum value of d for all possible unit normal vector of shear band is greater than or equal to 1 indicates the bifurcation and the corresponding unit normal vector is considered to be admissible.
Note that the solution of Equation (18) is proportional to the magnitude jump of the strain rate ‖̇+‖ − ‖̇−‖. As a result, d is independent of the specific value of ‖̇+‖ − ‖̇−‖ and thus a unit value can be assumed during the calculation.
It is of interest to further discuss the above bifurcation criterion. There are three possibilities:  Figure 3(A). This is a general case. As discussed above, the same bifurcation condition can be drawn with Ref. [45], that is, the degree of bifurcation calculated by Equation (18) and Equation (22) is greater than or equal to unity.
(2) Case B: ( + ) is invertible and the directions oḟ+ anḋ− are the same, as shown in Figure 3(B). In this case, the degree of bifurcation is 1. Let us first recast constitutive equation (1) by virtue of Euler's theorem for homogeneous functions as follows: ( ) ⊗ ⃗̇) ∶̇=  ′ ( , ⃗̇) ∶̇ (23) in which ⃗̇=̇∕ ‖̇‖ denotes the direction of the strain rate;  ′ = ( ) + ( ) ⊗ ⃗̇c orresponds to the tangential constitutive tensor in the classic plasticity. As the directions oḟ+ anḋ− are the same in this case, the jump of stress rate, which is calculated by Equation (17), can now be rewritten as With Equation (17) replaced by Equation (24), the bifurcation condition can be derived as in which ′ is the acoustic tensor calculated by ′ = ′ , where ′ is the components of  ′ in Equation (23); Apparently, the bifurcation condition is ( ′ + ) = 0. The same bifurcation criterion as the one widely accepted for elastoplastic materials 17 is obtained. It is clear that the criterion proposed by Rudnicki and Rice 17 is a special case of that proposed by Wu. 45 (3) Case C: ( + ) is not invertible. The bifurcation condition is ( + ) = 0. Meanwhile, ‖̇+‖ should equal to ‖̇−‖ while the directions oḟ+ anḋ− should be the opposite, as shown in Figure 3(C). Note that in this case, the acoustic tensor is derived from the linear part of the constitutive equation (23), that is, . Moreover, the degree of bifurcation approaches infinity.

Failure surface in three-dimensional stress space
Stress probe in three-dimensional stress space is carried out to obtain the failure surface of the proposed model. To facilitate this, the three principal stresses are expressed by the three invariants, , and : in which is the radius of the stress vector on the deviatoric plane and is the Lode angle. During the stress probe process, the mean principal stress tr ∕3 is increased gradually from 1 kPa to 100 kPa. For each mean principal stress, the Lode angle is increased gradually from 0 to 2 . For each Lode angle, the length of stress vector on the deviatoric plane is increased gradually from zero, which corresponds to a hydrostatic stress state, until the following failure criterion is fulfilled 39 : where and are the matrix and vector form of linear and non-linear parts in the model. The following material parameters are considered first: ∕ = 170, = 0.2, = 35 • and = 0 • , where is the initial tangential modulus, is the reference initial mean stress, is the initial Poisson's ratio, is the frictional angle at critical state and is the dilatancy angle. With these four parameters, the parameters 1 ∼ 4 in Equation (13) are determined using a standard procedure described by Wu and Bauer. 39 The parameter in Equation (5) is first taken to be zero, which corresponds to the original model proposed by Wu et al., 48 and then taken to be 0.5, which corresponds to the new model with a general strength criterion. The stiffness function is taken to be 1. The obtained failure surfaces in the principal stress space are plotted in Figure 4. It is clear that the failure surface of the original model is a conical shape, which is similar to that of the Drucker-Prager model, while the failure surface of the new model inherits the properties of the general strength criterion. The deviatoric cross-sections with different are plotted in Figure 5(A). The deviatoric cross-section approaches a circular shape with approaching 0 • and a triangular shape with approaching 90 • , which is consistent with the SMP failure criterion. The deviatoric cross-sections with different for = 35 • are plotted in Figure 5(B). It is clear that with the new parameter introduced, a gradual transition between the Mises criterion and the SMP criterion can be obtained, which is useful to describe the soil behaviour in laboratory tests more precisely.  The following material parameters are used: ∕ = 170, = 0.2, = 35 • , = 0 • and = 0.5. Meanwhile, the parameter is also taken to be zero to obtain results corresponding to the original model for comparison. Again the stiffness function is taken to be 1. The obtained stress-strain and volumetric curves are shown in Figure 6. By comparison between results with = 0.0 and = 0.5, it is clear that more reasonable peak stress ratios can be obtained with the incorporation of the general strength criterion.

True triaxial tests with different values of
As described in the above section, the general bifurcation criterion for hypoplastic materials proposed by Wu 45 can be further divided into three particular cases: (1) Criterion A: the same with the general bifurcation criterion; (2) Criterion B: acoustic tensor criterion with tangential constitutive tensor; (3) Criterion C: acoustic tensor criterion with linear constitutive tensor. These three particular criteria are used, respectively, to predict the onset of shear band. It is found that for Criterion C no shear band can be detected. That is to say, ( + ) in Equation (18) is always invertible. For Criteria A and B, shear bands can be detected at some values of . The results are plotted in Figure 6. As expected, the onset of shear band predicted by Criterion B is generally later than that predicted by Criterion A. The reason is that fewer assumptions are used in Criterion A. It will be further clarified in the next section that Criterion A is more applicable for fitting the experimental data. The acoustic tensor criterion with tangential constitutive tensor, 17 which is generally adopted for an elastoplastic model, is inadequate for a hypoplastic model.
Note that in the elastoplastic approach, the acoustic tensor criterion with tangential constitutive tensor, which is analogous to Case B in the present study, is applicable for continuous bifurcation. Besides the continuous bifurcation, the discontinuous bifurcation has been studied as well. For example, Borja 64 studied the discontinuous bifurcation with the assumption of finite strain. Rice and Rudnicki 65 studied the discontinuous bifurcation with the assumption that elastic unloading occurs outside shear band while plastic loading occurs inside shear band. The discontinuous bifurcation is to some extent similar to Case A in the present study. However, there are some differences. The assumption that elastic unloading occurs outside shear band while plastic loading occurs inside shear band makes an elastoplastic model incrementally non-linear. Actually, a hypoplastic model is incrementally non-linear as long as the direction of strain changes, which is more general than the elastic unloading. Moreover, in terms of the prediction of the onset of shear band, if an elastoplastic model is used, localization first occurs as a continuous bifurcation. As the deformation proceeds an infinitesimal amount past the continuous bifurcation point, elastic unloading occurs and discontinuous bifurcation may occur. 65 On the contrary, our hypoplastic study shows that bifurcation in Case A occurs always earlier than that in Case B, that is, discontinuous bifurcation occurs always before continuous bifurcation.
With Criterion A, the stress-strain curve under homogeneous deformation at = 0.5 is plotted in Figure 7(A). The corresponding degree of bifurcation value in Equation (22) is also plotted. As can be seen from Figure 7(A), the degree of bifurcation increases gradually with the loading process, and exceeds unity when the axial strain exceeds approximately 1.2%, indicating the onset of bifurcation as marked by a filled square. After this point, the degree of bifurcation further increases slightly and then approaches a constant value. It is of interest to further study the admissible inclinations and orientations of shear band obtained by the proposed approach, as defined in Figure 1 and Equation (21). Figure 7 Figure 8 for comparison. Note that both the experimental and numerical results are normalized to make a comparison between them. As shown in Figure 8, shear banding occurs prior to homogeneous failure for a larger range of Lode angles, while the homogeneous failure occurs first for a narrower range near triaxial compression, which is consistent with the existing experimental findings. 4,5,20 22,45 When shear banding occurs prior to homogeneous failure, the bifurcation surface agrees well with the experimental results. The authors believe that it is not a coincidence but an indication of shear band formation. After the onset of shear band, the strain-softening occurs, and the peak strength is confirmed. As shown in Figure 8

Comparison against true triaxial tests with different values of
The proposed approach is also applied to simulate the true triaxial tests on dense sand with different values of intermediate principal stress ratio . True triaxial tests on Santa Monica Beach sand were conducted by Wang and Lade. 5 The tall rectangular prismatic specimens with a width of 76 mm, a length of 76 mm and a height of 188 mm are utilized, which made the shear bands develop more freely and softening behaviour more pronounced. A constant effective confining pressure of 49 kPa was used in all tests. can also be well predicted. For the true triaxial test at = 0.0, which corresponds to the triaxial compression test, no bifurcation point is predicted, which coincides with the experimental findings. For true triaxial tests at other values, the bifurcation points obtained by the proposed approach generally fit well with the peak point in the laboratory tests. Note that bifurcation analysis provides only the necessary condition for shear banding. In other words, the shear band may initiate beyond but not prior to the bifurcation point. This may explain why the model generally predicts a slightly earlier onset of bifurcation compared to the experimental results. As shown in Figure 10, the axial strains at the onset of shear band predicted by the proposed approach are compared with those in the laboratory tests. The results obtained by elastoplastic approaches proposed by Huang et al. 20 are also plotted. Experimental results reveal that the axial strains at the onset of the shear band first decrease rapidly and then tend to be stable with increasing value of intermediate principal stress ratio . As reported by Huang et al., 20 for the elastoplastic approach, the non-coaxial behaviour should be incorporated to obtain a similar trend. However, for the proposed hypoplastic approach, as the non-coaxial behaviour is implied in the constitutive model, no special non-coaxial treatment is required. Using the proposed hypoplastic approach, the axial strains at the onset of the shear band are predicted with satisfactory accuracy.

CONCLUSIONS
An approach based on hypoplasticity is proposed to predict the onset of shear band in sand under true triaxial conditions. First, a recently proposed basic model is refined to better describe the pre-bifurcation behaviour of sand under true triaxial conditions. The revisions include two parts. A general strength criterion is incorporated to consider the strength variation with different Lode angles, while a new stiffness function is proposed to better describe the stiffness variation during the loading process. Then, a general bifurcation criterion for hypoplastic materials is derived and further discussed, based on which the bifurcation analysis approach is completed for predicting the onset of shear band in sand under true triaxial conditions. The performance of the proposed approach is evaluated. Special attention is given to the comparison with classical elastoplastic approaches. It is found that the acoustic tensor criterion with tangential stiffness tensor, which is generally adopted for an elastoplastic model, is inadequate for a hypoplastic model, whereas the general bifurcation criterion for hypoplastic materials is applicable. Comparisons with experimental results show that the proposed approach can not only well capture the pre-bifurcation stress-strain relationship at all values of intermediate principal stress ratio but also provide a satisfactory prediction of the onset of shear band. No special non-coaxial treatment is required to yield a reasonable variation trend of bifurcation strain with . For most values except for near zero, the onset of shear

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.