Predictingtheonsetofshearbandisofsignificanceinunderstandingthefailure of geomaterials. The prediction accuracy is dictated by the constitutive model usedforthedescriptionofthepre-bifurcationbehaviour.Inthisstudy,wefirst modify a recently proposed hypoplastic constitutive model by incorporating a general strength criterion and a stiffness function. We proceed to consider the onsetofshearbandinsandundertruetriaxialconditions.Wedemonstratethat 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 trendofbifurcationstrainwithintermediateprincipalstressratio 𝑏. KEYWORDS bifurcation,hypoplastic,localization,sand,shearband,truetriaxial 1INTRODUCTION Shearbandisoneofthemostcommonlyobservedfailuremodesingeomaterials.Thedefiningfeatureofshearbandis the transition from a homogeneous deformation pattern to an inhomogeneous one with the emergence of a system of localizedzonesorbandswithexcessiveshearstrainwithinthem.1Theonsetofshearbandisgenerallyfollowedbythe overall strength reduction of geomaterials. Meanwhile, the existing research has revealed that in some common cases theshearbandcanoccurduringthehardeningstageofgeomaterials,thatis,beforethepeakstrengthofhomogeneous deformationpatternisachieved.2–8Therefore,itisofimportanceinunderstandingthefailureofgeomaterialsaswellas geotechnicalapplicationstoaccuratelypredicatetheonsetofshearband.9 Duringthelastdecades,therehavebeensomeexperimentalstudiesontheshearbandformationingeomaterials,espe- cially in sands,10–16in which various key factors that influence the strain localization and the onset of shear band are Thisisanopenaccessarticleunderthetermsofthe CreativeCommonsAttribution License,whichpermitsuse,distributionandreproductioninanymedium,providedthe originalworkisproperlycited. ©2021TheAuthors. InternationalJournalforNumericalandAnalyticalMethodsinGeomechanics publishedbyJohnWiley&SonsLtd. IntJNumerAnalMethodsGeomech. 2021;1–16. wileyonlinelibrary.com/journal/nag 1 2 ZHANGetal. investigated,includingboundarycondition,samplesize,thesand’sdensityandparticlesize,confiningpressure,imper- fectionandheterogeneity,inherentandstress-inducedanisotropy,andsoforth. Accompanied by the experimental studies on the shear band, some theoretical studies have been conducted to detecttheonsetofshearband,amongwhichthemostnotableisthatbyRudnickiandRice.17Theonsetofshearband 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 strainrateinaplanarbandcanbeadmissible.Thesingularityofthedeterminantofacoustictensor,whichisassociated withtheelastoplastictangentialmodulusandtheinclinationofshearband,denotestheonsetofshearband.Withthis 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 constitutivemodels,includingnonnormality,noncoaxial,sub-loadingsurfacemodelandotherextensions,1,18–27ofwhich thenoncoaxialplasticmodelsmaybethemostsuccessful.20–23,27 Hypoplastic constitutive models, as an alternative to elastoplastic models, have aroused increasing interests recently.28–41Some beneficial features of hypoplastic constitutive models include simple formulation, fewer param- eters, 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 elastoplasticmodels.Ashypoplasticmodelsareincrementalnon-linear,thebifurcationcriterionshouldbedifferent.42–44 Tamagninietal.42presentedacriticalreviewofresultsobtainedbasedontwoparticularclassesofhypoplasticconstitutive models,whicharethemostpopularonesatthattime.Wu45proposedabifurcationcriterionfornon-linearanalysiswith no additional assumptions. The criterion provides a lower bound solution and is more general than the incremental linear ones. Huang46proposed general bifurcation conditions for a non-polar hypoplastic model and its micro-polar continuum extension. Huang47investigated 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 predictingaccuracyoftheonsetofshearbandisgreatlydeterminedbythepre-bifurcationconstitutivemodel.Withthe development of hypoplastic constitutive models, the older ones have been gradually replaced by novel ones (e.g. Refs. [28,48]).However,thepredictingcapabilityoftheonsetofshearbandofthesenovelmodelsremainsunclear. Thisstudyaimsatpredictingtheonsetofshearbandinsandundertruetriaxialconditionswithintheframeworkofa recentlyproposedhypoplasticbasicmodelforsand.48Inordertowellcapturethepre-bifurcationstress–strainbehaviour of sand under true triaxial conditions, the basic model is extended to describe the sand behaviour under true triaxial conditionsbyincorporatingageneralstrengthcriterionproposedbyYaoetal.49Furthermore,theconstitutiveformulation isrevisedbyintroducingastiffnessfunctiontobetterfitthetangentialstiffnessofsandundertruetriaxialconditions.With thenovelhypoplasticmodelconstructed,thebifurcationcriterionproposedbyWu45isreviewedandfurtherdiscussed. EmphasisisputonthecomparisonwiththewidelyacceptedacoustictensorcriterionproposedbyRudnickiandRice.17 Themodelresponses,intheaspectsofbothdescribingthehomogeneousmechanicalbehaviourandpredictingtheonset ofshearbandaredemonstrated.Consequently,thenumericalresultsbytheproposedapproacharecomparedwiththe laboratoryobservationsofsandundertruetriaxialconditions. Notationsandconventions :Compactorindextensorialnotationisusedthroughout.Vectorsandsecond-ordertensors aredenotedwithboldletters(e.g. 𝝈,𝒗)andfourth-ordertensorsaredenotedwithcalligraphicboldletters(e.g. ).Symbols ‘⋅’and‘∶’betweentensorsofvariousordersdenoteinnerproductwithsingleanddoublecontractions,respectively.The dyadic product of two tensors is indicted by ‘ ⊗’.‖̇𝜺‖=√ tr(̇𝜺2)denotes the Euclidean norm and the arrow operator is definedas ⃗̇𝜺=̇𝜺∕‖̇𝜺‖.Followingtheconventionofsoilmechanics,compressionstressandcontractionstrainaretakento bepositive. 2HYPOPLASTICCONSTITUTIVEMODEL Inahypoplasticconstitutivemodel,thestressrateinresponsetomaterialdeformationisexpressedasanon-lineartensor- valuedfunctionofCauchystressandstrainrate.Ageneralhypoplasticframeworkfordescribingrate-independentmate- rialbehaviourcanbewrittenas39,50 ̊𝝈=(𝝈) ∶ ̇𝜺+𝐍 ( 𝝈 ) ‖̇𝜺‖, (1) ZHANGetal. 3 wheretheterms and𝐍,respectively,denotethelinearandnon-linearcomponents, 𝝈anḋ𝜺,respectively,standforthe Cauchystresstensorandthestretchingstrainratetensor.TheJaumannstressrateisdefinedby ̊𝝈=̇𝝈+𝝈⋅ ̇𝒘− ̇𝒘⋅𝝈, (2) where𝒘denotesthespintensor.Thestretchingtensorandspintensorarerelatedtothevelocitygradienttensorthrough ̇𝜺=1 2(∇𝐯 +𝐯∇), ̇ 𝒘=1 2(∇𝐯 −𝐯∇), (3) where𝐯isthevelocityand ∇isthegradientoperator. Withintheaboveframework,anovelbasichypoplasticconstitutivemodelforsandhasbeenrecentlyproposedbyWu etal.48Inthenovelbasicmodel,therearethreelineartermsandonenon-linearterm: ̊𝝈=𝑐1(tr𝝈)̇𝜺+𝑐2𝝈(tṙ𝜺)+𝑐3tr(𝝈̇𝜺) tr𝝈𝝈 +𝑐4(𝝈+𝝈∗)‖̇𝜺‖,(4) in which 𝑐𝑖(𝑖=1,2,3,4) are dimensionless parameters. 𝝈∗denotes the deviatoric stress tensor defined by 𝝈∗=𝝈− (tr𝝈)𝜹∕3with𝜹beingtheKroneckerdeltatensor. Theabovebasicmodelhasafailuresurfaceofaconicalshape,whichissimilartothatoftheDrucker–Pragermodel. Inordertodescribethesoilbehaviourundertruetriaxialconditions,ageneralstrengthcriterionproposedbyYaoetal.49 is incorporated in this study. According to the criterion, noting that the yield locus on the deviatoric plane of various materialsaregenerallyinsidetheMisescircleandoutsidethespatiallymobilizedplane(SMP),51theradiusoftheyield locusforaspecificLodeanglecanbetakentobetheweightedaverageofthoseoftheMisescriterionandtheSMPcriterion. Following this general strength criterion, the ratio between the radius of failure surface on the deviatoric plane for a specificLodeangleandtheradiusundertriaxialcompressionstatecanbecalculatedby, 𝑟=𝛼 𝑞𝑚∕𝑞𝑠+1−𝛼 (5) inwhich 𝛼isaweightedfactorlyingbetween0and1,with 𝛼=0correspondingtotheMisescriterionand 𝛼=1corre- spondingtotheSMPcriterion.Byadjustingtheweightvalue 𝛼between0and1,acontinuoustransitionfromtheMises criteriontoSMPcriterioncanbeobtained.Givenanarbitrarystressstate, 𝑞𝑚and𝑞𝑠arethecorrespondinglengthoffail- uresurfaceonthedeviatoricplaneundertriaxialcompressionstressstatefortheMisescriterionandtheSMPcriterion, respectively. 𝑞𝑚and𝑞𝑠canbecalculatedbyRef.[ 49], 𝑞𝑚=√ 𝐼2 1−3𝐼2, 𝑞𝑠=2𝐼1 3√ (𝐼1𝐼2−𝐼3)∕(𝐼1𝐼2−9𝐼3)−1,(6) where𝐼1,𝐼2and𝐼3arethestressinvariants.NotethatthespecificformulaisslightlydifferentfromthatinRef.[ 49],asthe so-calledtransformedstressspaceisusedintheirwork. Inordertoobtainafailuresurfacewiththeabovestrengthcriterioninhypoplasticity,letusconsiderthesphericaland deviatoricpartsoftheconstitutiveequation( 4),respectively.Thesphericalpartis tr̊𝝈=𝑐1(tr𝝈)(tṙ𝜺)+𝑐2(tr𝝈)(tṙ𝜺) +𝑐3tr(𝝈̇𝜺)+𝑐4tr𝝈‖̇𝜺‖.(7) Notethatinthefourthterm tr𝝈∗vanishessince tr𝝈∗=0.Thedeviatoricpartis ̊𝝈∗=𝑐1(tr𝝈)̇𝜺∗+𝑐2𝝈∗(tṙ𝜺)+𝑐3tr(𝝈̇𝜺) tr𝝈𝝈∗ +2𝑐4𝝈∗‖̇𝜺‖.(8) Notingthattr ̇𝜺=0atcriticalstate,bylettingtr ̊𝝈=0,Equation( 7)canberewrittenas 𝑐3tr(𝝈∗̇𝜺)+𝑐4(tr𝝈)‖̇𝜺‖=0 . (9) 4 ZHANGetal. Notingthattr ̇𝜺=0andtr(𝝈̇𝜺)=tr(𝝈∗̇𝜺)atcriticalstate,byletting ̊𝝈∗=0,Equation( 8)canberewrittenas 𝑐1(tr𝝈)̇𝜺∗+𝑐3tr(𝝈∗̇𝜺) tr𝝈𝝈∗+2𝑐4𝝈∗‖̇𝜺‖=0 . (10) AccordingtoEquations( 9)and(10),notingthattr 𝝈inthedeviatoricplaneremainsaconstant,thedeviatoricstressat failuresurfaceforaspecificLodeanglecanbeadjustedbydividingthefourthtermby 𝑟anddividingthethirdtermsby 𝑟2. Withthissimpleapproach, 𝝈∗ofthefailureismultipliedby 𝑟whichisassociatedwiththeLodeangle,whiletr 𝝈remains unchanged.Theradiusonthedeviatoricplaneoffailuresurfaceisthusadjustedtofulfilthegeneralstrengthcriterion. Therefore,theextensionofthebasichypoplasticmodelconsideringthegeneralstrengthcriterioncanbewrittenas, ̊𝝈=𝑐1(tr𝝈)̇𝜺+𝑐2𝝈(tṙ𝜺)+𝑐3 𝑟2tr(𝝈̇𝜺) tr𝝈𝝈 +𝑐4 𝑟(𝝈+𝝈∗)‖̇𝜺‖.(11) Ithasbeenreportedthatinsomecasesthehypoplasticmodelspredictimproperlythedevelopmentofstiffnesswith shearstrain.Acommoncaseisthatthemodelpredictsasharpbendofthestress–straincurvebeforethecriticalstateis achieved,whereastheexperimentalresultsshowamoregradualstiffnessdecrease.28,40Toovercomethisshortcoming,Wu andKolymbus40suggestedmultiplyingastiffnessfunctioninthemodel.Inthisstudy,thestiffnessfunctionismodified tobetterfitthetruetriaxialcase,whichiscalculatedas 𝐼𝑠𝑒=𝑟𝛽 (1+𝑟𝑠)2, (12) where𝑟𝑠denotesthestressratio ‖𝝈∗‖∕(tr𝝈),𝑟istheratiodefinedinEquation( 5)and𝛽isaparametertoreflecttheeffect degreeoftheratio 𝑟. Notethatmultiplyingthewholeconstitutivebyascalingfunctiondoesnotaffectthefailuresurface.Amoresophis- ticatedconstitutivemodelcanbeconstructedbyconsideringtheevolutionofvoidratio,whichcanbestraightforward.48 However,theauthorsselecttokeeptheconstitutivemodelsimplebecauseinsufficientexperimentaldataonshearband areaccessible.Consequently,therefinedconstitutivemodelproposedinthisstudycanbewritteninthefollowingcom- pleteform: ̊𝝈=𝐼𝑠𝑒[𝑐1(tr𝝈)̇𝜺+𝑐2𝝈(tṙ𝜺)+𝑐3 𝑟2tr(𝝈̇𝜺) tr𝝈𝝈 +𝑐4 𝑟(𝝈+𝝈∗)‖̇𝜺‖].(13) Therearetotallysixparametersinthemodel.Thefourparameters 𝑐1∼𝑐4arethesamewiththeoriginalmodel,which can be determined using a standard procedure described by Wu and Bauer.39The additional two parameters 𝛼and𝛽 arerelatedtostrengthandstiffness,respectively,whichoffermoreflexibilitytofitmodelpredictiontoexperimentaldata undertruetriaxialconditions.Itshouldbepointedoutthattheratedependence52,53andanisotropy54–58isoutofthescope ofthisstudy. Thenumericalintegrationofthehypoplasticconstitutiveequation( 13)ismorestraightforwardwhencomparedwith theelastoplasticmodel,duetothefactthatthehypoplasticmodelonlyhasasingleequation.35Sincehypoplasticmodel doesnothaveayieldsurface,thestressreturnmappingalgorithmscommonforelastoplasticmodels59–61arenotneeded. In thisstudy,a simple explicit forwardEulerschemewith constantstep sizes62,63is adopted. Thisschemecan achieve satisfactory accuracy as long as the step size is small enough, which can be determined by a trial and error approach. Sinceonlyonestresspointisconsideredhere,thetimecostfornumericalintegrationislow.Morecomplicatednumerical integrationschemesforhypoplasticmodelareavailableinRef.[ 35]. 3BIFURCATIONCRITERIAFORTHEONSETOFSHEARBAND Wu45proposedabifurcationcriterionforhypoplasticmaterials.Inthissection,thebifurcationcriterionisbrieflyreviewed. Thensomefurtherdiscussionsaremade.Emphasisisputonthecomparisonwiththeacousticcriterion,whichhasbeen widelyacceptedasageneralbifurcationcriterionforelastoplasticmaterials(e.g.Refs.[ 1,20–23,25,27 ]). ZHANGetal. 5 FIGURE 1 Schematicinterpretationofshearband The schematic interpretation of shear band is shown in Figure 1. The development of a shear band bifurcation in a homogeneouslydeformingmaterialcanbecompletelydefinedwiththefollowingkinematicconditions ̇𝜺+=̇𝜺−+̇𝜺𝑗=̇𝜺−+1 2(𝐠⊗𝐧+𝐧⊗𝐠), ̇𝒘+=̇𝒘−+̇𝒘𝑗=̇𝒘−+1 2(𝐠⊗𝐧−𝐧⊗𝐠),(14) in which the superscripts +’a n d−’ refer to values inside and outside the shear band, respectively, the superscript }𝑗’ denotesthejumpofthevalueoftheshearband, 𝐧istheunitnormalvectortotheshearbandand 𝐠isavectordefining thevelocitygradientjumpinthedirection 𝐧. Inaddition,therequirementofcontinuingequilibriumattheonsetofshearbandimpliesthefollowingstaticcondition: ̇𝝈+𝐧−̇𝝈−𝐧=0 (15) inwhich ̇𝝈+anḋ𝝈−aretheCauchystressratesinsideandoutsidetheband,respectively. Notethatthestressitselfiscontinuousbeforebifurcation,andthusthejumpoftheJaumannstressratecanbeobtained accordingtoEquation( 2)as ̊𝝈𝑗=̊𝝈+−̊𝝈−=̇𝝈+−̇𝝈−+𝝈⋅ ̇𝒘𝑗−̇𝒘𝑗⋅𝝈. (16) Meanwhile,thestressratesinsideandoutsideshearbandcanbeobtainedfromtheconstitutiveEquation( 1) ̊𝝈+=(𝝈) ∶ ̇𝜺++𝐍(𝝈)‖̇𝜺+‖, ̊𝝈−=(𝝈) ∶ ̇𝜺−+𝐍(𝝈)‖̇𝜺−‖.(17) 6 ZHANGetal. FIGURE 2 Geometricalmeaningofthejumpofthestrainrateacrossshearband ThebasicequationforbifurcationanalysiscanbeobtainedbyincorporatingEquations( 16)and(17)intoEquation( 15), (𝐋+𝐀)𝐠=−( 𝐍(‖̇𝜺+‖−‖̇𝜺−‖) (18) inwhich 𝐋istheacoustictensorcalculatedby 𝐿𝑗𝑘=𝑛𝑖𝐿𝑖𝑗𝑘𝑙𝑛𝑙,where𝐿𝑖𝑗𝑘𝑙isthecomponentof inEquation( 1);𝐀isa second-ordertensorassociatedwithfinitestrainandiscalculatedby 𝐴𝑗𝑘=1 2(−𝑛𝑗𝑛𝑖𝜎𝑖𝑘+𝑛𝑖𝜎𝑖𝑟𝑛𝑟𝛿𝑗𝑘+𝑛𝑖𝜎𝑖𝑗𝑛𝑘−𝜎𝑘𝑗) (19) while𝐍iscalculatedby 𝐍=𝑛𝑖𝑁𝑖𝑗. (20) Byintroducingtheorientationangleandinclinationangleofshearband,asdefinedinFigure 1,thecomponentsofthe unitnormalvector 𝑛𝑖oftheshearbandcanbecalculatedasfollows: 𝑛1= cos𝛽cos𝛼, 𝑛2=c o s 𝛽 s i n 𝛼 , 𝑛3=s i n 𝛽 ,(21) where𝛼and𝛽are,respectively,theorientationandinclinationanglesofshearband. Theonsetofshearbandmeansthatasolutionotherthan 𝐠≡0existsforEquation( 18).Toresolvetheproblemweneed tohaveacloserlookatthejumpofstrainrateacrosstheshearband.AsshowninFigure 2,thejumpofstrainratetensor mustfulfilthetriangleinequality,whichreads, 𝑑=‖̇𝜺+−̇𝜺−‖ abs(‖̇𝜺+‖−‖̇𝜺−‖)≥1, (22) where𝑑istermedasthedegreeofbifurcation. Thebifurcationproblemcanbesolvedwiththeregula-falsimethod.Specifically,theunitnormalvector(i.e.theorien- tationandinclination)ofshearbandisfirstassumedandthenthedegreeofbifurcation dcanbeobtainedaccordingto Equation(18)andEquation( 22).Themaximumvalueof dforallpossibleunitnormalvectorofshearbandisgreaterthan orequalto1indicatesthebifurcationandthecorrespondingunitnormalvectorisconsideredtobeadmissible. NotethatthesolutionofEquation( 18)isproportionaltothemagnitudejumpofthestrainrate ‖̇𝜺+‖−‖̇𝜺−‖.Asaresult, disindependentofthespecificvalueof ‖̇𝜺+‖−‖̇𝜺−‖andthusaunitvaluecanbeassumedduringthecalculation. Itisofinteresttofurtherdiscusstheabovebifurcationcriterion.Therearethreepossibilities: ZHANGetal. 7 FIGURE 3 Geometricalmeaningofthejumpofthestrainrateacrossshearbandinterface:(A)CaseA;(B)CaseB;(C)CaseC (1) Case A: (𝐋+𝐀)is invertible and the directions of ̇𝜺+anḋ𝜺−are not the same, as shown in Figure 3(A).T h i si s a general case. As discussed above, the same bifurcation condition can be drawn with Ref. [ 45], that is, the degree of bifurcationcalculatedbyEquation( 18)andEquation( 22)isgreaterthanorequaltounity. (2) Case B: (𝐋+𝐀)is invertible and the directionsof ̇𝜺+anḋ𝜺−are the same, as shown in Figure 3(B). In this case, thedegreeofbifurcation 𝑑is1.Letusfirstrecastconstitutiveequation( 1)byvirtueofEuler’stheoremforhomogeneous functionsasfollows: 𝐍(𝝈)⊗⃗̇𝜺)∶ ̇𝜺=′(𝝈,⃗̇𝜺)∶ ̇𝜺 (23) inwhich⃗̇𝜺=̇𝜺∕‖̇𝜺‖denotesthedirectionofthestrainrate; ′=(𝝈)+𝐍(𝝈)⊗ ⃗̇𝜺correspondstothetangentialconstitu- tivetensorintheclassicplasticity.Asthedirectionsof ̇𝜺+anḋ𝜺−arethesameinthiscase,thejumpofstressrate,which iscalculatedbyEquation( 17),cannowberewrittenas ̊𝝈𝑗=((𝝈)+𝐍(𝝈)⊗⃗ 𝜺)∶̇𝜺𝑗=′(𝝈,⃗𝜺)∶̇𝜺𝑗. (24) WithEquation( 17)replacedbyEquation( 24),thebifurcationconditioncanbederivedas (𝐋′+𝐀)𝒈=0 (25) in which 𝐋′is the acoustic tensor calculated by 𝐿′𝑗𝑘=𝑛𝑖𝐿′ 𝑖𝑗𝑘𝑙𝑛𝑙, where𝐿′ 𝑖𝑗𝑘𝑙is the componentsof ′in Equation ( 23); Apparently,thebifurcationconditionis 𝑑𝑒𝑡(𝐋′+𝐀)=0.Thesamebifurcationcriterionastheonewidelyacceptedfor elastoplasticmaterials17isobtained.ItisclearthatthecriterionproposedbyRudnickiandRice17isaspecialcaseofthat proposedbyWu.45 (3) Case C: (𝐋+𝐀)is not invertible. The bifurcation condition is 𝑑𝑒𝑡(𝐋+𝐀)=0 . Meanwhile, ‖̇𝜺+‖should equal to ‖̇𝜺−‖whilethedirectionsof ̇𝜺+anḋ𝜺−shouldbetheopposite,asshowninFigure 3(C).Notethatinthiscase,theacoustic tensorisderivedfromthelinearpartoftheconstitutiveequation( 23),thatis, .Moreover,thedegreeofbifurcation 𝑑 approachesinfinity. 4MODELRESPONSE 4.1Failuresurfaceinthree-dimensionalstressspace Stressprobeinthree-dimensionalstressspaceiscarriedouttoobtainthefailuresurfaceoftheproposedmodel.Tofacilitate this,thethreeprincipalstressesareexpressedbythethreeinvariants, 𝑡𝑟𝝈,𝜃and𝑙: 𝜎𝑖=√ 2 3𝑙sin(𝜃+2𝜋(2−𝑖) 3)+tr𝝈 3, (26) 8 ZHANGetal. FIGURE 4 Failuresurfacesinthree-dimensionalstressspace where 𝑙=1√ 3√ (𝜎1−𝜎2)2+(𝜎2−𝜎3)2+(𝜎3−𝜎1)2 𝜃 = arctan(1√ 32𝜎2−𝜎1−𝜎3 𝜎1−𝜎3)(27) inwhich𝑙istheradiusofthestressvectoronthedeviatoricplaneand 𝜃istheLodeangle. Duringthestressprobeprocess,themeanprincipalstresstr 𝝈∕3isincreasedgraduallyfrom1kPato100kPa.Foreach mean principal stress, the Lode angle 𝜃is increased gradually from 0 to 2𝜋. For each Lode angle, the length of stress vectoronthedeviatoricplane 𝑙isincreasedgraduallyfromzero,whichcorrespondstoahydrostaticstressstate,untilthe followingfailurecriterionisfulfilled39: (𝐋−1𝐍)𝑇𝐋−1𝐍−1=0, (28) where𝐋and𝐍arethematrixandvectorformoflinearandnon-linearpartsinthemodel. Thefollowingmaterialparametersareconsideredfirst: 𝐸𝑖∕𝜎𝑐= 170,𝜈𝑖=0 . 2,𝜙𝑐𝑟𝑡=3 5◦and𝜓=0◦,where𝐸𝑖isthe initialtangentialmodulus, 𝜎𝑐isthereferenceinitialmeanstress, 𝜈𝑖istheinitialPoisson’sratio, 𝜙𝑐𝑟𝑡isthefrictionalangle at critical state and 𝜓is the dilatancy angle. With these four parameters, the parameters 𝑐1∼𝑐4in Equation ( 13)a r e determinedusingastandardproceduredescribedbyWuandBauer.39Theparameter 𝛼inEquation( 5)isfirsttakentobe zero,whichcorrespondstotheoriginalmodelproposedbyWuetal.,48andthentakentobe0.5,whichcorrespondstothe newmodelwithageneralstrengthcriterion.Thestiffnessfunction 𝐼𝑠𝑒istakentobe1.Theobtainedfailuresurfacesin theprincipalstressspaceareplottedinFigure 4.Itisclearthatthefailuresurfaceoftheoriginalmodelisaconicalshape, whichissimilartothatoftheDrucker–Pragermodel,whilethefailuresurfaceofthenewmodelinheritsthepropertiesof thegeneralstrengthcriterion.Thedeviatoriccross-sectionswithdifferent 𝜙𝑐𝑟𝑡areplottedinFigure 5(A).Thedeviatoric cross-section approaches a circular shape with 𝜙𝑐𝑟𝑡approaching 0◦and a triangular shape with 𝜙𝑐𝑟𝑡approaching 90◦, whichisconsistentwiththeSMPfailurecriterion.Thedeviatoriccross-sectionswithdifferent 𝛼for𝜙𝑐𝑟𝑡=3 5◦areplotted inFigure5(B).Itisclearthatwiththenewparameter 𝛼introduced,agradualtransitionbetweentheMisescriterionand theSMPcriterioncanbeobtained,whichisusefultodescribethesoilbehaviourinlaboratorytestsmoreprecisely. 4.2Truetriaxialtestswithdifferentvaluesof 𝒃 Numericalsimulationsarecarriedoutontruetriaxialtestswithdifferentvaluesofintermediateprincipalstressratio 𝑏, whichisdefinedby 𝑏=( 𝜎2−𝜎3)∕(𝜎1−𝜎3),inordertofurtherdemonstratetheperformanceoftheproposedapproach. ZHANGetal. 9 FIGURE 5 (A)failuresurfaceonthedeviatoricplanewithvarying 𝜙𝑐𝑟𝑡;(B)failuresurfaceonthedeviatoricplanewithvarying 𝛼(𝜙𝑐𝑟𝑡= 35◦) FIGURE 6 (A) Stress–strain curve with different values of intermediate principal stress ratio 𝑏(𝛼=0 . 0); (B) volumetric strain versus axialstrain( 𝛼=0 . 0);(C)Stress–straincurvewithdifferentvaluesof 𝑏(𝛼=0 . 5);(D)volumetricstrainversusaxialstrain( 𝛼=0 . 5) Thefollowingmaterialparametersareused: 𝐸𝑖∕𝜎𝑐= 170,𝜈𝑖=0 . 2,𝜙𝑐𝑟𝑡=3 5◦,𝜓=0◦and𝛼=0 . 5.Meanwhile,theparam- eter𝛼isalsotakentobezerotoobtainresultscorrespondingtotheoriginalmodelforcomparison.Againthestiffness function𝐼𝑠𝑒is taken to be 1. The obtained stress–strain and volumetric curves are shown in Figure 6. By comparison between results with 𝛼=0 . 0and𝛼=0 . 5, it is clear that more reasonable peak stress ratios can be obtained with the incorporationofthegeneralstrengthcriterion. 10 ZHANGetal. Asdescribedintheabovesection,thegeneralbifurcationcriterionforhypoplasticmaterialsproposedbyWu45canbe furtherdividedintothreeparticularcases:(1)CriterionA:thesamewiththegeneralbifurcationcriterion;(2)Criterion B:acoustictensorcriterionwithtangentialconstitutivetensor;(3)CriterionC:acoustictensorcriterionwithlinearcon- stitutivetensor.Thesethreeparticularcriteriaareused,respectively,topredicttheonsetofshearband.Itisfoundthat forCriterionCnoshearbandcanbedetected.Thatistosay, (𝑳+𝑨)inEquation( 18)isalwaysinvertible.ForCriteriaA andB,shearbandscanbedetectedatsomevaluesof 𝑏.TheresultsareplottedinFigure 6.Asexpected,theonsetofshear bandpredictedbyCriterionBisgenerallylaterthanthatpredictedbyCriterionA.Thereasonisthatfewerassumptions areusedinCriterionA.ItwillbefurtherclarifiedinthenextsectionthatCriterionAismoreapplicableforfittingthe experimentaldata.Theacoustictensorcriterionwithtangentialconstitutivetensor,17whichisgenerallyadoptedforan elastoplasticmodel,isinadequateforahypoplasticmodel. 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 bifurca- tion, the discontinuous bifurcation has been studied as well. For example, Borja64studied the discontinuous bifurca- tion with the assumption of finite strain. Rice and Rudnicki65studied the discontinuous bifurcation with the assump- tion that elastic unloading occurs outside shear band while plastic loading occurs inside shear band. The discontin- uous 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 predic- tion of the onset of shear band, if an elastoplastic model is used, localization first occurs as a continuous bifurca- tion. As the deformation proceeds an infinitesimal amount past the continuous bifurcation point, elastic unloading occurs and discontinuous bifurcation may occur.65On 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 . 5is plotted in Figure 7(A).T h e correspondingdegreeofbifurcation 𝑑valueinEquation( 22)isalsoplotted.AscanbeseenfromFigure 7(A),thedegree ofbifurcation 𝑑increasesgraduallywiththeloadingprocess,and 𝑑exceedsunitywhentheaxialstrainexceedsapproxi- mately1.2%,indicatingtheonsetofbifurcationasmarkedbyafilledsquare.Afterthispoint,thedegreeofbifurcation 𝑑 furtherincreasesslightlyandthenapproachesaconstantvalue.Itisofinteresttofurtherstudytheadmissibleinclinations andorientationsofshearbandobtainedbytheproposedapproach,asdefinedinFigure 1andEquation( 21).Figure7(B) showsthevariationofadmissibleinclinationsofshearbandwithdifferentaxialstrains.Attheonsetofbifurcation,the admissible inclination of shear band is about 53◦. With the increase of axial strain, more values of the inclination of shearbandareadmissible.Thefinalrangeoftheinclinationofshearbandliesbetween 46◦and64◦.Figure7(C)shows the admissible orientations of shear band. A symmetric range is obtained, which lies in the range of −32◦∼3 2◦and 148◦∼ 212◦. 5COMPARISONAGAINSTTRUETRIAXIALEXPERIMENTALRESULTS 5.1Comparisonagainststrengthonthedeviatoricplane Goldscheider66conductedaseriesofstrain-controlledtruetriaxialtestsondenseKarlsruhesand.Thestrengthobtainedon thedeviatoricplaneisobtained,asshowninFigure 8.Theparametersusedinthenumericalsimulationare 𝐸𝑖∕𝜎𝑐= 170, 𝜈𝑖=0 . 2,𝜙𝑐𝑟𝑡=3 0◦,𝜓=0◦and𝛼=0 . 5.Theobtainedresultsoffailuresurfaceandbifurcationsurfaceonthedeviatoric planeareshowninFigure 8forcomparison.Notethatboththeexperimentalandnumericalresultsarenormalizedto makeacomparisonbetweenthem.AsshowninFigure 8,shearbandingoccurspriortohomogeneousfailureforalarger rangeofLodeangles,whilethehomogeneousfailureoccursfirstforanarrowerrangeneartriaxialcompression,which isconsistentwiththeexistingexperimentalfindings.4,5,20 22,45Whenshearbandingoccurspriortohomogeneousfailure, thebifurcationsurfaceagreeswellwiththeexperimentalresults.Theauthorsbelievethatitisnotacoincidencebutan indicationofshearbandformation.Aftertheonsetofshearband,thestrain-softeningoccurs,andthepeakstrengthis confirmed.AsshowninFigure 8(B),theadmissibleinclinationsofshearbandwithdifferentLodeanglesliebetween 52◦ and54◦. ZHANGetal. 11 FIGURE 7 Numericalresultsat 𝑏=0 . 5:(A)stressandbifurcationdegree 𝑑;(B)inclinationofshearband;(C)orientationofshearband 5.2Comparisonagainsttruetriaxialtestswithdifferentvaluesof 𝒃 Theproposedapproachisalsoappliedtosimulatethetruetriaxialtestsondensesandwithdifferentvaluesofintermediate principal stress ratio 𝑏. True triaxial tests on Santa Monica Beach sand were conducted by Wang and Lade.5The tall rectangularprismaticspecimenswithawidthof76mm,alengthof76mmandaheightof188mmareutilized,which made the shear bands develop more freely and softening behaviour more pronounced. A constant effective confining pressureof49kPawasusedinalltests. 12 ZHANGetal. FIGURE 8 Resultsonthedeviatoricplane:(A)failuresurface,bifurcationsurfaceandexperimentalresults66;(B)variationofshearband inclinationwithLodeangle 𝜓 Theparametersusedinthenumericalsimulationsare 𝐸𝑖∕𝜎𝑐= 900,𝜈𝑖=0 . 4,𝜙𝑐𝑟𝑡=4 2◦,𝜓=24◦,𝛼=0.5and𝛽=1 2. Theobtainednumericalresultswithtruetriaxialtestsatvariousvaluesof 𝑏(𝑏 = 0.0,0.20,0.41,0.60,0.80,1.0 )areshown inFigures9(A)to9(F).Thecomparisonsagainsttheexperimentalresultsshowthattheproposedhypoplasticmodelcan wellcapturethehomogeneousstress–strainbehaviourbeforeshearbandformation.Meanwhile,theonsetofshearband can also be well predicted. For the true triaxial test at 𝑏=0 . 0, which corresponds to the triaxial compression test, no bifurcationpoint ispredicted,whichcoincideswiththeexperimentalfindings.Fortruetriaxial testsatother 𝑏values, thebifurcationpointsobtainedbytheproposedapproachgenerallyfitwellwiththepeakpointinthelaboratorytests. Notethatbifurcationanalysisprovidesonlythenecessaryconditionforshearbanding.Inotherwords,theshearband mayinitiatebeyondbutnotpriortothebifurcationpoint.Thismayexplainwhythemodelgenerallypredictsaslightly earlier onset of bifurcation compared to the experimental results. As shown in Figure 10, the axial strains at the onset ofshearbandpredictedbytheproposedapproacharecomparedwiththoseinthelaboratorytests.Theresultsobtained byelastoplasticapproachesproposedbyHuangetal.20arealsoplotted.Experimentalresultsrevealthattheaxialstrains at the onset of the shear band first decrease rapidly and then tend to be stable with increasing value of intermediate principalstressratio 𝑏.AsreportedbyHuangetal.,20fortheelastoplasticapproach,thenon-coaxialbehaviourshouldbe incorporatedtoobtainasimilartrend.However,fortheproposedhypoplasticapproach,asthenon-coaxialbehaviouris impliedintheconstitutivemodel,nospecialnon-coaxialtreatmentisrequired.Usingtheproposedhypoplasticapproach, theaxialstrainsattheonsetoftheshearbandarepredictedwithsatisfactoryaccuracy. 6CONCLUSIONS An approach based on hypoplasticity is proposed to predict the onset of shear band in sand under true triaxial condi- tions.First,arecentlyproposedbasicmodelisrefinedtobetterdescribethepre-bifurcationbehaviourofsandundertrue triaxialconditions.Therevisionsincludetwoparts.Ageneralstrengthcriterionisincorporatedtoconsiderthestrength variationwithdifferentLodeangles,whileanewstiffnessfunctionisproposedtobetterdescribethestiffnessvariation duringtheloadingprocess.Then,ageneralbifurcationcriterionforhypoplasticmaterialsisderivedandfurtherdiscussed, basedonwhichthebifurcationanalysisapproachiscompletedforpredictingtheonsetofshearbandinsandundertrue triaxialconditions. The performance of the proposed approach is evaluated. Special attention is given to the comparison with classical elastoplasticapproaches.Itisfoundthattheacoustictensorcriterionwithtangentialstiffnesstensor,whichisgenerally adoptedforanelastoplasticmodel,isinadequateforahypoplasticmodel,whereasthegeneralbifurcationcriterionfor 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 𝑏.F o rm o s t 𝑏values except for 𝑏near zero, the onset of shear ZHANGetal. 13 FIGURE 9 Comparisonbetweenexperimental5andnumericalresultsoftruetriaxialtestsonSantaMonicaBeachsandwithdifferent valuesofintermediateprincipalstressratio 𝑏 14 ZHANGetal. FIGURE 10 Axialstrainattheonsetofshearbandversusintermediateprincipalstressratio 𝑏 bandoccurspriortothehomogeneouspeakstrength,whichfurtherindicatestheimportanceofbifurcationanalysisof shearband. ACKNOWLEDGEMENTS TheresearchissupportedbytheNaturalScienceFoundationofGuangdongProvince(GrantNo.2018A030310346),the WaterConservancyScienceandTechnologyInnovationProjectofGuangdong(2020-11),theH2020MarieSkłodowska- CurieActionsRISE2017HERCULES(778360)andFRAMED(734485),theErasmus +KA2projectRe-built(2018-1-RO01- KA203-049214),theNazarbayevUniversityResearchFund(SOE2017001)andtheNaturalScienceFoundationofChina (GrantNo.41807223). DATA AVAILABILITY STATEMENT Thedatathatsupportthefindingsofthisstudyareavailablefromthecorrespondingauthoruponreasonablerequest. ORCID WeiZhang https://orcid.org/0000-0002-3933-0109 WeiWuhttps://orcid.org/0000-0002-0286-0720