A non-orthogonal elastoplastic constitutive model for the sand with dilatancy is presented 4 in the characteristic stress space . Dilatancy of sand is represented by the direction of plastic flow , which 5 can be directly determined b y applying the non -orthogonal plastic flow rule to an improved elliptic yield 6 function . A new hardening paramet er is developed to de scribe the contractive and dilative volume change 7 during the shear process , which is coordinate d with the non-orthogonal plastic flow direction . The 8 combination of the non-orthogonal plastic flow rule and the proposed hardening parameter renders the 9 proposed model with the ability to reasonably describe t he stress -strain relationship of sand with 10 dilatancy . The model performance is evaluated by comparing with the experimental data of sand under 11 triaxial conditions. 12 Keywords : Sand; Dilatancy; Constitutive model; Non-orthogonal plastic flow rule ; Hardening parameter ; 13 Characteristic stress 14 15 16 1 Ph.D. Candidate. Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, 17 Beijing University of Technology, Beijing 100124, China. E -mail: liangjingyu@emials.bjut.edu.cn (J.Y. 18 Liang) . 19 2* Professor. Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, 20 Beijing University of Technology, Beijing 100124, China (Corresponding author). E -mail: 21 dechun@bjut.edu.cn (D.C. Lu). 22 3 Professor. Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, 23 Beijing University of Technology, Beijing 100124, China. E -mail: duxiuli@bjut.edu.cn (X.L. Du). 24 4 Profes sor. Institute of Geotechnical Engineering, University of Natural Resources and Life Sciences, 25 Vienna, 1180, Austria. E -mail: wei.wu@boku.ac.at (W. Wu). 26 5 Ph.D. Beijing Advanced Innovation Center for Future Urban Design, Beijing University of Civil 27 Engineering and Architecture, Beijing 102616, China. Email: machao@bucea.edu.cn (C. Ma). 28 This is a pre-print of the paper published in Computers and Geotechnics, Volume 118, February 2020, 103329 https://doi.org/10.1016/j.compgeo.2019.103329 2 1 In troduction 29 Dilatancy refer s to the volume change observed in granular materials during shear 30 deformation s [1-4]. It is an essential property of sand. Loose sand tends to show negative dilatancy 31 of volumetric contraction , while dense sand undergoes positive dilatancy [4, 5]. One of the 32 fundamental issues in modelling the mechanical behaviour of sand is to describe the dilatancy by 33 the ratio of the plastic volumetric strain increment to the plastic shear strain increment [6, 7]. Th e 34 ratio can be expressed by the direction of the plastic flow [8, 9]. Since the axis of the plastic strain 35 increment is assumed to be coaxial with that of stress , the plastic flow direction can be determined 36 by the plastic flow rule in the stress space [10-14]. The direction of the plastic strain increment is 37 usually determined based on the orthogonal gradient of the plastic potential function. Obviously, 38 this makes the plastic flow direction and its evolution rule inseparable from the plastic potential 39 surface. However, an appropriate plastic potential function is difficult to be directly constructed 40 based on experimental data, and some additional parameters may be require d [15-17]. In fact , the 41 purpose of constructing the plastic potential function is to produce the proper plastic flow direction . 42 As an alternative, the plastic flow direction can also be directly achieved by applying the non - 43 orthogo nal plastic flow rule to the yield function [18-20]. The non -orthogonal plastic flow rule has 44 been used to model the dilatancy behaviours of soil under triaxial compression conditions [20-23]. 45 However, for three -dimensional (3D) stress conditions, the non -orthogonal plastic flow rule has not 46 been well applied to simulate positive dilatancy [19]. H ow to use the non-or thogonal plastic flow rule 47 to describe the negative/ positive dilatancy of sand under 3D stress conditions and then establish a 48 3D non- or thogonal model for sand is worthy of further study. 49 Although the non-or thogonal plastic flow direction can be used to represent the stress -dilatancy 50 relationship for sand, it does not mean that the shear responses of the sand with dilatancy have been 51 totally captured. A proper hardening parameter is still needed in describing the contractive -dilative 52 behaviour of sand . I f hardening rule stems only from volume change , just as the models for normally 53 consolidated clay [21, 24-26], once the critical state of zero incremental dilation is reached, the ultimate 54 plasti c deformation occurs and hardening stops [17]. Obviously , the plastic volumetric strain is not 55 suitable for direct use as a hardening parameter to predict t he dilative behaviour of sand . Ther efore, 56 the hardening parameter should be reconstructed . Several hardening parameters and hardening 57 modulus have been proposed and used for modelling dilatancy behaviours of sand that is related to 58 3 density and stress state [27-30], anisotropy [31-33] and even fabric evolution [34-36]. However, they are 59 all formul ated based on orthogonal plastic flow conditions . As a key point to reflect the dilatancy of 60 sand, a n appropriate hardening parameter should be coordinated with the plastic flow direction . That 61 is, the hardening parameters derived for the orthogonal plastic flow rule cannot be directly 62 introduced into a non -orthogonal model. Therefore, a hardening parameter should be reconstructed 63 to match up with the non -or thogonal plastic flow direction . 64 In this paper, the non -orthogon al plastic flo w rule is applied to capture the dilatancy behaviour 65 of sand. The plastic flow direction under 3D stress conditions is directly determined by applying the 66 non-orthogonal plastic flow rule to the elliptic yield function in the characteristic stress space. 67 Furthermore, a hardening parameter for sand with dilatancy is constructed to match up with the non- 68 orthogonal plastic flow direction . The capa bility of the proposed 3D model in modelling the shear 69 response of sand with the negative and positive dilatancy is verified by predicting test results of 70 sand. 71 2 Fr amework of 3D n on-orthogonal elastoplastic model in the c ij space 72 Based on the elastoplastic theory, the total strain εij and its increment d εij can be decomposed 73 into the reversible elastic strain component and the irreversible plastic strain component: 74 e p e p, d d dij ij ij ij ij ij         (1) 75 For the p lastic strain increment , the non-orthogonal plastic flow rule proposed by Lu et al . [19, 76 37] is adopted and can be expressed as : 77 pdkl ij ij kl ijSfmS   (2) 78 where Λ is the plastic multiplier, and mij (i, j=1, 2, 3) is the plastic flow direction. Skl (k, l=1, 2, 3) is 79 stress invariant that construct s the yield function and is associated with σij. The non-orthogonal 80 plastic flow rule is illustrated in Fig. 1: f irst determine the non -orthogonal gradient kl fS of 81 the yield function by the fractional derivative in the Skl coordinate system , and then transform the 82 non-orthogonal gradient in the Skl coordinate system to the σij coordinate system by the coordinate 83 transformation kl ijS . The non-or thogonal plastic flow rule in tensor form is the kernel for 84 determin ing the general expression of plastic strain increment, which is essential for the model 85 implementation [26, 38, 39]. The plastic flow direction mij is non-orthogonal to the yield surface, which 86 4 can also degenerate into the normal direction nij when μ= 1. 87 88 Fig. 1 Schematic diagram of the non -orthogonal plastic flow rule . 89 In order to make the proposed non -orthogonal model be born with the potential to reflect the 90 3D properties of soil, the concept of the characteristic stress cij is introduced [19, 40]. More details on 91 the definition of the cij and the corresponding variables are provided in Appendix A . The general 92 form of the yield function in the cij space can thus be expressed as : 93  p,0ij rs f f c H  (3) 94 where H is the hardening parameter indicating the hardening degree of soil due to the occurrence of 95 plastic strain p rs (r, s=1, 2 , 3). In view of the superior abi lities of c ij in reflect i ng 3D properties of 96 soil only through two characteristic stress invariants cn and cs [19], the yield function can be written 97 as: 98  p ns, , 0rs f f c c H  (4) 99 where cn and cs are the mean characteristic principal stress and the characteristic shear stress, 100 respectively. 101 By applying the non -orthogonal plastic flow rule in Eq. (2) t o the yield function in Eq. (4), 102 the plastic flow direction m ij in tensor form can be der ived in the cij space : 103 s n nsij ij ijc cffmc c c c      (5) 104 The plastic multiplier Λ can be determined by the consistency condition , which is the total 105 differential of the yield function f=0 in the cij space : 106 O 1 2 3 ijm   3Sm   1Sm   2Sm   3m   1m A   2m 3S 1S 2S   Yield surface: ij klSf 5 p pd d d 0ij rs ij rsf f HfccH       (6) 107 Where p rs is the plastic strain used for constructing the hardening parameter H. By substituting 108 pdrs rsm from Eq. (2) into Eq. (6), Λ can thus be expressed as : 109 pdij ij rs rsfcc fHmH   (7) 110 By substituting the plastic flow direction mij in Eq. (5) a nd t he plastic multiplier Λ in Eq. (7) 111 into Eq. (2), the plastic strain increment in the tensor form can be written as: 112 p s n ns pd dij ij ij ij ij ij rs rsfcc c cffmfH c c c cmH                (8) 113 For the elastic strain increment in Eq. (1), we can make use of the elastic theory, e.g. H ooke's 114 law: 115 e1d d dij ij kk ijEE    (9) 116 where ν is the Poisson’s ratio. E is the Young’s modulus, which can be obtained from the e-lnp 117 relationship of the swelling line under the is otropic compression conditions: 118  013 1 2eEp  (10) 119 where κ is the swelling index, and e 0 is the initial void ratio at p=p0. 120 Finally, by substitut ing Eqs. (9) and (8) i nto Eq. (1), the total strain increment is thus 121 determined with i n the framework of the 3D non-orthogonal elastoplastic model . 122 3 F ormulation of 3D non -orthogonal model 123 Within the framework of the non -orthogonal plastic theory, t he yield function plays the key 124 role in determin ing Λ and mij. In this paper, the improved elliptic yield function in the cij space is 125 adopted. 126  2 2 2 s n n n 0x f c N c c c    (11) 127 where N is the ratio of the vertical axis to the horizontal ax is of the elliptic yield curve as shown in 128 6 Fig. 2, which can be determined by the method in Section 3.4 bas ed on the phase transition stress 129 ratio F. Th e equivalent characteristic consolidation pressure cnx is the intersection of f with the mean 130 stress axis and defines the position of the yield surface , w hich can be expressed by the current stress 131 state as follows : 132 2 2 2 sn n 2 nxc N ccNc (12) 133 134 Fig. 2 Yield curve in the characteristic stress space. 135 Note that our yield function can also reduce to the yield function of the modified Cam -clay 136 model [25] when β=1 and N=M. As shown by Lu et al. [19], the axisymmetric ellipsoidal yield surface 137 in the cij space corresponds to a smooth and convex triaxial symmetric yield surface in the σij space. 138 That is, in the σij space, the asymmetric curve with respect to the p-axi s in the meridian plane and 139 the curved triangle in the deviatoric plane can both be directly described by the simple expression 140 in the cij space (i.e., Eq. (11)). Therefore, the θ effect and stress -induced anisotropy of geomateri als 141 can be reasonably captured by the yield function . 142 3.1 3D non -orthogonal plastic flow direction 143 In this paper, the non-orthogonal plastic flow direction is expressed based on the Caputo 144 fractional derivative operator. The specific expression s and properties of the Caputo fractional 145 derivative operator involved in this paper are provided in Appendix B . 146 By applying the non -orthogonal plastic flow rule in Eq. (2) to t he yield function in Eq. (11), 147 and further incorporating Eq. (12) an d the recursive properties of the Gamma function 148 Γ(x+1)= x·Γ(x), the non- or thogonal gradient m =(mn, ms) of the yield function can be expressed as : 149       2 2 2 2 ns s ns n2 2,,33N c c cmmc      m (13) 150 nxc nc sc n2xc F N O Yield curve ( )ijfc fF 1 7 Furthermore, b y combining Eq. (13) with the coordinate transformation coefficient ∂cn/∂cij=δij/3 151 and ∂cs/∂cij=3/2·(cij−cnδij)/cs, the plastic flow direction tensor mij can thus be obtained : 152       2 2 2 2 n ns s ns3 2 2 3 3 3 2ij ij ij ijcc N c c cmcc         (14) 153 3.1.1 Plastic flow direction in the c ij space 154 The evolution of the plastic flow direction with the expansion of the yield surface will be 155 illustrated , and its capability to model the dilatancy of sand will be explored as follows. To gain 156 perspectives, let us con sider the following m aterial parameters for demonstration purpose: M=0.95, 157 Mf=1.65 and β=0.4, which are typical for dense sand . Therefore, the phase t ransition stress ratio 158 F=0.366 and the failure stress ratio Ff =0.664 are correspondingly obtained based on expression s 159 listed in Tab. A -1 ( see Appendix A ), and N=0.62 and μ=0.52 can be determined by the method in 160 Section 3.4. T ypical stress path s under true triaxial conditions under p=196kPa at di fferent stress 161 Lode angles θ (θ=0°, 15° , 30°, 45° and 60°) are used for illustrating the variation rule of the plastic 162 flow direction. 163 The simple expression of the yield function makes the yield surface in a simple ellipsoid al 164 shape in the cij space , an d the ellipsoidal yield surface expands only when the size change s w ithout 165 changing shape , as shown in Fig. 3. The elliptic yield curve in the meridian plane and the circular 166 yield curve in th e deviatoric plane are shown in Fig. 4(a) and (b), r espectively . With the expansion 167 of the ellipsoidal yield surface, the plastic flow directions for μ=0.52 and μ=1.00 correspondingly 168 change, as shown in Fig. 3 and Fig. 4. The plastic flow direction non-orthogonal to the yield su rface 169 when μ=0.52 is represented by the red arrow, and the plastic flow direction orthogonal to the yield 170 surface when μ=1.00 is represented by the black arrow. As the shear stress increases in the meridian 171 plane , the non -orthogonal plastic flow direction shown in Fig. 4(a) changes from horizontal to the 172 right at χ=0 to vertically upward at χ=F for describing the negative dilatancy, and finally turns to 173 the direction pointing to the left until χ reaches Ff for capturing the positive dilatancy . As for the 174 plastic flow direction in the deviatoric plane, it is independent of the fractional -order μ, since no 175 variable associated with the stress Lode angle con tains in the yield function. That is, the plastic flow 176 direction s in the deviatoric plane for μ=0.52 coincide with the plastic flow direction s for μ=1. Here, 177 the coincident plastic flow directions are represented by the red arrows, as shown in Fig. 4(b). The 178 plastic flow direction in the deviatoric plane is the radial direction along the current characteristic 179 8 stress Lode angle θβ defined in Tab. A -1 ( see Appendix A ). The θβ-values on the stress paths under 180 the condition that θ=0° and 60° are constant, while θβ-values on the characteristic stress paths at 181 θ=15°, 30° and 45° are changeable and increase with increasing shear stress. Therefore, the plastic 182 flow direction s along the radial direction change with θβ-v alues accordingly , as shown in Fig. 4(b) . 183 184 Fig. 3 Evolution of the n on-orthogonal plastic flow direction with the expanding yield surface in the cij 185 space. 186 0 5 10 15 20 251050510 Stress pathPlastic flow direction μ =1.00 μ =0.52F Ext.Comp.cs (kPa) cn (kPa)Ff (cnf, csf)=(7.70, 5.13) Stress pathθ β = 38.9°θ β = 23.2° θ β = 60°θ β = 50.3°c* z c* xc* y θβ = 0° (cn, cs)=(8.08, 2.96) 187 (a) (b) 188 Fig. 4 Yield curve and n on-orthogonal plastic flow direction in (a) the meridian plane of the cij space 189 and (b) the deviatoric plane of the cij space . 190 3.1.2 Plastic flow direction in the σ ij space 191 The expansion of the ellipsoidal yield surface in the cij space corresponds to the evolut ion of 192 the yield surface in the σij space , as shown in Fig. 5. The plastic flow direction determined in the cij 193 1c 2c 3c nc 1.00 0.52   9 space is applied to the corresponding stress state point on the yield surface in the σij space . 194 Corresponding t o the expansion of the elliptic yield curve shown in Fig. 4(a), t he evolution of the 195 unsymmetrical yield curves in the meridian plane of the σij space is shown in Fig. 6(a). A s the shear 196 stress increases along the stress path under triaxial compression conditions with a constant p, t he 197 plastic flow direction in Fig. 6(a) cha nges in the same way as that in Fig. 4(a). It changes f rom 198 horizontal at stress ratio η=0 to vertically upward at η=M and finally t urns to the direction pointing 199 to the left until η reaches Mf for capturing the dilatancy behaviour from negative to positive . In the 200 deviatoric plane of the σij space , the yield curve tends to change from an approximate circle to a 201 curved triangle with the increase of η, and the plastic flow directions at different θ-values 202 correspondingly change, as shown in Fig. 6(b). The plastic flow directions at θ=0° and 60° coincide 203 with the radial directions, while the plastic flow directions at θ=15°, 30° and 45° deviate leftward 204 from the radial direction, and the degree of the deviation increases with the increase of η. This means 205 the stress -induced anisotropy associated with the yield behaviour and plastic flow direction is 206 reflected. 207 208 Fig. 5 Evolution of the yield surface and the n on-orthogonal plastic flow direction in the σij space. 209 210 1 2 3 p 1.00 0.52   10 0 4 8 12 16630369 Stress pathPlastic flow direction μ =1.00 μ =0.52M Ext.Comp.q (102 kPa) p (102 kPa)Mf (p0, qf)=(196,323.4) Stress pathθ=60°θ=45°θ=30°θ=15°θ=0°* z * x* y(p0, q)=(196,186.2) 211 (a) (b) 212 Fig. 6 Yield curve and n on-orthogonal plastic flow direction in (a) the meridian plane of the σij space 213 and (b) the deviatoric plane of the σij space . 214 The above analysis illustrates that the non- orthogonal plastic flow rule [19] and the characteristic 215 stress [40] are complementary in capturing soil behaviours under 3D stress conditions. The non- 216 orthogonal plastic flow rule makes the plastic flow direction suitable for capturing the dilatancy of 217 sand from negative to positive , while the characteristic stress extends the ability to captur e dilatancy 218 to 3D stress conditions . The combination of these two concepts gives rise to a new approach to 219 capture the pl astic flow direction for the sand under 3D stress conditions. 220 3.2 Hardening parameter for sand with dilatancy 221 Hardening parameter is an internal variable reflecting the hardening behaviours during the 222 loading process . I t is closely related to the determination of the plastic multiplier Λ in Eq. (7). In 223 order to incorporate the dilatancy of sand from negative to positive into a non -orthogonal model, 224 the hardening parameter should be formulated to be coordinated with the plastic flow direction . 225 During the construct ion process , the hardening parameter should also be independent of the stress 226 path [8, 29]. O bviously, the plastic volumetric strain as a hardening parameter is effective for normally 227 consolidated clay . However, it is not suitable for describing initia l negative and subsequent positive 228 dilatancy of sand . By referring to the method proposed by Yao et al. [8], a stress -path-dependency 229 factor Ω is introduced and combine d w ith the volumetric strain to construct the hardening parameter 230 H, which aims to make H stress -path-independent and capable of capturing the dilatancy of sand. 231 232 11 p vd H (15) 233 For normally consolidated clay, p v is used as the hardening parameter , so Ω=1 is obt ained . 234 The relationship p v p n n0 lnx l c c can be derived based on the linear relationship between e and 235 ln(cn) under isotropic compression conditions [19], w here lp=(λ−κ)/[β·(1+e0)], and λ is the 236 compression index . Therefore, t he equivalent characteristic consolidation pressure can be expressed 237 as  p n n0 v p expxc c l  . For sand, H in Eq. (15) i s adopted to replace p v, cnx can thus be expressed 238 as: 239 p n n0 n0 v pp1exp exp dxHc c cll            (16) 240 The partial derivative of the yield function f with respect to p v is required for determining Λ in Eq. 241 (7) and can be obtained based on Eqs. (11) and (16): 242 2 n nn pp v n v px x xc ffN c ccl       (17) 243 In addition, t he following differential equation is also needed and can be deri ved from Eq. (11): 244  2 n n n s s d 2 d 2 dij x ijfc N c c c c cc   (18) 245 By substitu t ing Eqs. (17) and (18) together with Eq. (12) into Eq. (7), t he specific expression of 246 Λ is obtained : 247    2 2 2 n s n s n s p 2 2 2 n s n nd 2 d N c c c c c cl c c N c m   (19) 248 By substitut ing Λ in Eq. (19) and (mn, ms) in Eq. (13) into Eq. (8), together with the stress 249 ratio χ =cn/cs, the pl astic volumetric strain and the plastic shear strain can be expressed as: 250          22 ns p p vn 22 n 222 ns p p ss 22 22 nd 2 d d d 2 d 2d2N c clm cN N c clmN cN                   (20) 251 The two stress paths shown in Fig. 7, include the constant mean characteristic principal stress 252 path AB and the isotropic compression stress path AC, ar e used to derive the stress- path- dependent 253 12 factor Ω and an explicit expression of H. 254 255 Fig. 7 Stress path s for deriving the hardening parameter H. 256 3.2.1 Constant mean characteristic principal stress path AB (dcn=0) 257 Along the stress path AB, dcn=0 and cn is a constant. Therefore, p sd in Eq. (20) can be 258 simplified as: 259    3 p p s 2 2 2 24d 1d 2l NN         (21) 260 For the normally consolidated clay, Ω=1. p sd is expressed as : 261    3 p p s 2 2 2 24dd 2l NN         (22) 262 Similar to Eq. (22), the expression of p sd for sand can be assumed as: 263    3 p p s 2 2 2 2 ff4dd 2l NN         (23) 264 Here, ρ is an undetermined coefficient that reflects the ratio of the magnitude of p sd between clay 265 and sand along the stress path AB. N in Eq. (22) is changed to Nf that related to the failure stress 266 ratio Ff of sand. In fact, Ff is not constant in describing the shearing process with some complex 267 behaviours, but for the sake of simplicity, it i s assumed to be a constant . When sand reache s its 268 failure state at χ=Ff, the expression of p sd in Eq. (23) approach es infinity , which is similar to Eq. 269 (22) for clay . Thus, ff 21 NF  can be obtained. 270 Furthermore, b y combining Eq. (21) with Eq. (23), the stress -path-dependent factor Ω is 271 A nc sc F O Yield curve ( )ijfc fF B C 0f 1f 13 derived and can be expressed as : 272        2 2 2 2 ff 2 2 2 2221 2FF NN                (24) 273 Substituting Eq. (24) into Eq. (15) yields : 274        2 2 2 2 ff p v 2 2 2 2221d 2FF H NN                (25) 275 3.2.2 Isotropic compression stress path AC (dcs=0) 276 Under isotropic compression condition s (dcs=0), χ is equal to 0. Eq. (25) can thus be simplif ied 277 as follows : 278  2 4 f p v 2421dFHN (26) 279 Along stress path AC, only the plastic volumetric strain p v occurs for sand. That is, p vd H 280 can be used to depict the hardening rule of sand under isotropic compression condition s. The 281 expression of ρ can thus be derived by combining p vd H with Eq. (26). 282  24 2 4 f1 2N F   (27) 283 By substituting Eq. (27) i nto Eq. (24), the stress -path- dependent factor Ω can be expressed as : 284        2 2 2 24 ff 4 2 2 2 2 f2 2 2FF N FNN              (28) 285 Finally, the new hardening parameter is obtained by substituting Eq. (28) into Eq. (15) and can be 286 expressed as follows: 287        2 2 2 24 ff p v 4 2 2 2 2 f2 d2 2FF NHFNN             (29) 288 With the proposed hardening parameter H, the negative and positive dilatancy behaviours of 289 sand can be described. Typical results of sand behaviours under triaxial compression conditions 290 under p=constant are illustrated by the schematic diagram shown in Fig. 8(a), which is arranged in 291 terms of the stress ratio η=q/p and the volumetric strain εv vs. the shear strain εs. Since p=constant, 292 e0ij and p ij ij are obtained. The negative and positive dilatancy can be explain ed and 293 indicated by the following equation derived from Eq. (29). Mor eover, correspond ing to Fig. 8(a), 294 14 the following conclusions listed in Tab. 1 can be drawn from Eq. (30) based on the non -negative 295 dH. 296        2 2 2 2 4 f p v 4 2 2 2 2 ff2 2dd 2NNFHN FF               (30) 297 Tab. 1 V ariation rule of p vd with the increase of stress ratio η 298 Stress state Characteristic stress state dH d𝜀vp Current condition η=0 χ=0 dH= d𝜀vp d𝜀vp>0 Isotropic compression 0<η<M 0<χ<F dH>0 d𝜀vp>0 Negative dilatancy η=M χ=F dH>0 d𝜀vp=0 Phase transition state M<η≤Mf F<χ≤Ff dH>0 d𝜀vp<0 Positive dilatancy By adopting H with the properties indicated above, the shear responses transiting from 299 contraction to dilation are captured. When M f=M, Ff=F c an be obtained, H will degenerate to p v, 300 only negative dilatancy will be described , as shown in Fig. 8(b). For the condition that β=1, cij 301 becomes σij, F=M and Ff=Mf are obtained. When N=M, μ=1 can be derived , so H can degenerate 302 to the hardening parameter proposed by Yao et al. [8]. 303 304 (a) (b) 305 Fig. 8 Stress -strain relationship of soil with: (a) negative and positive dilatancy; and (b) negative 306 dilatancy . 307 O M  s fM v Transition state condition Negative dilatancy Positive dilatancy Failure stress condition O  s M v Negative dilatancy Failure stress condition f MM 15 3.2.3 Stress -path- independency of the hardening parameter 308 The stress -path- independence of the hardening parameter H can be interpreted the change of 309 H along different paths , as shown in Fig. 9 and Fig. 10. The yield curve passes through the stress 310 state point A0(cnx0, 0) which is the only intersection point of the yield surface f0 and the cn-axis. 311 According to Eq. (16), the unique intersection point corresponds to a uniqu e value of H0. That is, 312 the yield surface f0 is actually an equivalent surface of H0. As the loading continues along diff erent 313 paths to points B1, C1 a nd D1 that are all on the yield surface f1 as shown in Fig. 9(a), the yield surface 314 gradually expands, and the intersection moves along the cn-axis to the point A1 (cnx1, 0). Accordingly, 315 the value of the hardening parameter changes to H1, which is equivalent for any stress state on f1. 316 Therefore , the difference Δ H=H1-H0 can be calculated and is independent on the stress path . This is 317 exactly the basic properties that hardening parameters need to satisfy . The corresponding yield 318 surface and stress paths in the σij space are shown in Fig. 9(b). Fig. 10 show s that , as the sand is 319 loaded along different stress paths from f0 to f1, t he value of the H changes from H0 to H1. The results 320 of equal variation in H shown in F ig. 10 indicate s that H i s stress -path- independence . 321 322 (a) (b) 323 Fig. 9 Stress paths for verification of the stress -path-independence of H: (a) in the characteristic stress 324 space; (b) in the ordinary stress space. 325 1c 2c 3c nc 0A 1A 1B 1C 1D C D B 1 2 3 p C 0A 1D 1C 1B D 1A B 16 326 Fig. 10 Stre ss -path-independence of H. 327 3.3 Stress -strain relationship 328 With the newly derived stress -path-dependent factor Ω, the explicit expression of the plastic 329 multiplier Λ can be obtained based on Eq. (19): 330           2 2 2 2 2 2 2 4 n s n s n s f p 4 2 2 2 2 2 2 2 f f n s n n2 d 2 d 2 2N N N c c c c c c FlN F F c c N c m                   (31) 331 By substituting Eqs. (14) and (31) into Eq. (2), the plastic strain increment is thus derived: 332            4 2 2 2 2 1 f n s n p p 22 4 2 2 2 2 n f f2 2 d 2 d 3 d32 2ij ij ij ijF N N c c cc lN N c F F                               (32) 333 The expression of plastic strain increment is totally expressed by the characteristic stress, and thus 334 naturally, the β criterion is used as the shear failure condition of sand. When Mf=M, Ff=F can be 335 obtained, and Eq. (32) can thus be simplified to determine the plastic strain increment for the loose 336 sand with only negative dilatancy . 337 The total strain increment produced by the current stress increment can thus be obtained by Eq. 338 (1), where t he elastic strain increment for sand is expressed by Eq. (9) and the plastic strain increment 339 is expressed by Eq. (32). The parameters involved in the non-or thogonal elastoplastic model can be 340 easily determined, which will be discussed in the following section . 341 3.4 Parameter calibration 342 There are seven independent material parameters in the proposed model, i.e., λ, κ, ν, φ, φf, β 343 and ξ. All material parameters have clear physical meaning and can be easily determined by 344 conventional triaxial tests. The parameters λ, κ and ν are the same with the parameters of the 345 17 modified Cam -clay model . The parameters φ and φf are the internal friction angle s of the sand under 346 triaxial compression conditions at the phase transition state and the failure state , respectively . β is 347 the characteristic stress index, which can be determined based on the expression of the characteristic 348 failure stress ratio Ff listed in Tab. A-1 ( see Appendix A ). Actually, Ff can be determin ed by 349 combining the internal friction ang le φfb at any intermediate principal stress coefficient b (b=(σ2– 350 σ3)/(σ1–σ3) ) with β based on the expression of the β criterion in Eq. (39). Note t he value of Ff is 351 independent of b, i.e., 12ffbbFF , where the subscripts b1 and b2 are two different b-values. Based 352 on this, the parameter β can be determined by 1b and 2b , which reflects the relative proportional 353 relationship of shear strength at different b. In particular, φf at b1=0 and φfe at b2=1 are used to 354 determine β-value based on the method proposed by Lu et al. [19, 40]. P arameter ξ is called the spacing 355 ratio [21, 30]. It can also be used to describe the change of the void ratio ΔΠ of sand between the initial 356 state and the phase transition state under triaxial compressi on conditions under p=c onst . The change 357 of the void ratio has been derived and expressed as ΔΠ=(λ−κ)·ln(px/px0) [21], where px0 and px are the 358 equivalent mean principal stresses for the initial state and the phase transition state, respectively . 359 Here, px/px0=ξ is assumed, therefore, ξ can be determined by ξ=exp[ ΔΠ/(λ−κ)], and ξ=2 is a special 360 case derived from the modified Cam -clay model . 361 In addition, N and μ are t wo model parameters, which can both be easily determined based on 362 the above material parameters. Their determination is given below . 363 Since the equivalent mean principal stresses px0 and px are both defined for the is otropic 364 compression conditions , the equivalent mean characteristic principal stress es can be expressed as 365 cnx0= pr·(px0/pr)β and cnx= pr·(px/pr)β. Therefore, based on ξ=px/px0 mentioned above, we can get 366 cnx=ξβ·cnx0. On the other hand, substituting the phas e transition stress state ( cnt, cst)= (cnt, F·cnt) into 367 Eq. (12) yields: 368 2 n nt 21xFccN  (33) 369 By combining c nx=ξβ·cnx0 with Eq. (33), N can be derived as: 370 nt n 0 ntxcNFcc (34) 371 where cnt=(c1t+2c3t)/3, in which c1t=pr·(σ1t/pr)β and c3t=pr·(σ3t/pr)β can be obtained from Tab. A -1. σ1t 372 and σ3t are t he major and minor principal stresses at the phase transition state, which can be 373 18 calculated from σ1t=px0( 1+2M/3) and σ3t=px0( 1−M/3), respectively. Therefore, parameter N can be 374 determined by the following equation: 375        1 2 3 2 1 3 3 1 2 3 2 1 3MMNF MM         (35) 376 The shape of the yield surface is known when N is determined . Therefore, t he expression of 377 plastic flow direction is obtained , and the stress -dilatancy ratio d can be correspondingly derived : 378  22 p v p2 s2 d d2Nd    (36) 379 When the phase transition state reaches , χ reaches F, and d i n Eq. (36) i s equal to 0, then the 380 fractional -order μ can be determined by the following equation : 381 2 222F NF (37) 382 4 M odel verification 383 The performance of the proposed model can be evaluated and verified by comparing model 384 prediction with test results of sand from literature . W e consider the tests on saturated Toyoura sand 385 under conventional triaxial compression ( σ3=p0, p=p0 a nd σ1=p0) and extension ( σ3=p0, p=p0 and 386 σ1=p0) conditions and true triaxial conditions (θ=15°, 30° and 45° ) [29]. The mean diameter of sand 387 is 0.2mm, the uniformity coefficient is 1.3 and the specimens were prepared in medium dense with 388 an initial void ratio 0.68. The material parameters used for model prediction are listed in Tab. 2. The 389 corresponding model parameter s N=0.439 and μ=0.390 can thus be determined based on the 390 material parameters by the method mentioned in Section 3.4. 391 Tab. 2 Material p arameters used for model prediction. 392 λ/(1+e0) κ/(1+e0) ν φ φf β ξ 0.00403 0.00251 0.3 24.209 ° 40.475 ° 0.24 2.3 4.1.1 Toyoura sand under conventional triaxial conditions 393 Test data of saturated Toyoura sand under conventional triaxial compression (σ3=p0, p=p0 and 394 σ1=p0) and extension (σ3=p0, p=p0 and σ1=p0) c onditions are digitized from [8] and indicated by open 395 circles , as shown in Fig. 11. As shown in Fig. 11(a) and (b) , t he negative dilatancy of sand is obvious, 396 since p increases under triaxial compression and extension conditions when σ3=p0. Less negative 397 dilatancy can be found in Fig. 11(c) and (d), since a constant p i s applied. Moreover, due to the 398 19 decrease of p under triaxial compression and extension conditions when σ1=p0, the least negative 399 dilatancy is found in Fig. 11(e) and (f). The test results of Toyoura sand along six different stress 400 paths are predicted, as shown in Fig. 11(a)-(f). Due to the contributions of the non -orthogonal plastic 401 flow direction and the newly constructed H, different degrees of dilatancy from negative to positive 402 along six different stress paths are well predicted. Comparisons between the test results and the 403 prediction show that the relationship of ε1 vs. σ1/σ3 and ε1 vs. σ1/σ3 is also well captured by the 404 proposed model. On the other hand, benefit from the application of characteristic stress, the β 405 criterion is naturally used in the model, and thus the tendency that the shear strength of the sand 406 under extension conditions is greater than that under compression conditions can be well reflected. 407 σ1/σ3 1 -4 -2 0 2 435 ε3/%Toyoura sand prediction test data v /% ε1/%Compression σ3=196 kPa-4 -2 0 -4 -2 0 2 43 ε3/%Toyoura sand prediction test datav /% ε1/%σ1/σ3 Extension σ3=196 kPa 15 -4 -2 0 408 (a) Compression when σ3=p0 (b) Extension when σ3=p0 409 σ1/σ3 -4 -2 0 2 43Compression p=196 kPa ε3/%Toyoura sand prediction test datav /% ε1/% 15 -4 -2 0 -4 -2 0 2 4ε3/%Toyoura sand prediction test datav /% ε1/%σ1/σ3 Extension p=196 kPa 135 -4 -2 0 410 (c) Compression when p=p0 (d) Extension when p=p0 411 20 Compression σ1=196 kPa -4 -2 0 2 43 ε3/%Toyoura sand prediction test datav /% ε1/%σ1/σ3 15 -4 -2 0 -4 -2 0 2 4ε3/%Toyoura sand prediction test datav /% ε1/%σ1/σ3 Extension σ1=196 kPa 135 -4 -2 0 412 (e) Compression when σ1=p0 (f) Extension when σ1=p0 413 Fig. 11 Comparison between predicted and test results under conventional triaxial conditions . 414 4.1.2 Toyoura sand under true triaxial conditions 415 True triaxial drained tests under constant mean principal stress p=196kPa [29] can be used to 416 demonstrate the capability of the 3D elastoplastic models for sand [29, 41]. The values of θ are 15°, 417 30° and 45°, respectively . Triaxial tests under p=196 w hen θ=0° (triaxial compression ) and 60° 418 (triaxial extension ) have been used to verify the proposed non -orthogonal model, as shown in Fig. 419 11(c) and (d) . The capability of the proposed model in capturi ng the behaviour of sand under 3D 420 stress conditions can be verified by the predicting test results of sand, as shown in Fig. 12(a) -( c). 421 The comparisons are arrange in terms of the relation among the principal stress ratio σ1/σ3, 422 volumetric strain εv, and principal strains εi (i=1,2,3). In addition to well captur e the volume change 423 from negative to positive under 3D stress conditions, it is also critical for a 3D model for sand to 424 reflect three -dimensional strength behaviour and the three- dimensional plastic flow direction. For 425 the strength property of sand under 3D stress conditions ( θ=0°, 15°, 35°, 45° and 60 °), predictions 426 are in good agreement with the experiment results . Mor eover, an appropriate description of the 427 plastic flow direction in the deviatoric plane , as illustrated in Section 3.1, pr ovides the proposed 428 model with the ability to capture the relative relationship between different principal strains , as 429 shown in Fig. 12(a) -(c). 430 21 3 -4 -2 0 2 4Toyoura sand prediction test data εi /%p=196 kPa θ=15° 2-6 -4 -2 0v /%7σ1/σ3 15ε3ε2 ε1 431 (a) θ=15° 432 3 -4 -2 0 2 4Toyoura sand prediction test datav /% εi /%p=196 kPa θ=30°7σ1/σ3 15 2-6 -4 -2 0ε1ε2 ε3 433 (b) θ=30° 434 -4 -2 0 2 43ε3ε2 ε1 Toyoura sand prediction test data7 εi/%σ1/σ3 p=196 kPa θ=45° 1 2-6 -4 -2 05 v /% 435 (c) θ=45° 436 Fig. 12 C omparison between predicted and test results under true traixial conditions when p=p0. 437 5 Conclusions 438 In this paper, a non-orthogonal elastoplastic model was proposed within the framework of the 439 non-orthogonal elastoplastic model . Dilatancy behaviours of sand from initial negative to 440 subsequent positive along different stress paths were well described . Only seven material parameters 441 with clear physical meanings are needed and can be easily determined by conventional triaxial tests. 442 22 Model performance and comparison with test results of sand show that the proposed model is 443 capable of capturing the salient behaviour of sand under 3D stress conditions. The key 444 characteristics are listed as follows: 445 1. The framework of the 3D non-orthogonal elastoplastic model was established by introducing 446 the characteristic stress, which makes the developed model be born with the potential to 447 capture 3D mechanical behaviour of soil . 448 2. The plastic flow direction for describing the dilatancy of sand was determined by the non - 449 orthogonal plastic flow rule. Moreover, the complementarity between the non -orthogonal 450 plastic flow rule and the characteristic stress in capturing the 3D plastic flo w dir ection of sand 451 was illustrated. 452 3. A hardening parameter was newly constructed for capturing the negative and positive dilatancy 453 of sand during the shearing process, which is coordinated with the non -orthogonal plastic flow 454 direction . In addition, its independency from the stress path was demonstrated . 455 Acknowledgements 456 This study was supported by the National Key Research and Development Program of China 457 (Grant No. 2018YFC1504301 ), National Natural Science Foundation of China (Grant Nos. 458 51778026 , 5142 1005 and U1839201 ), an d the China Scho lar Council . The European Commission 459 for the support to the projects FRAMED “Fracture Across Scales and Materials, Processes and 460 Disciplines” (Grant No. 734485), and HERCULES “Towards geohazards resilient infrastructure 461 under changing climates” (Grant No. 778360) are gratefully acknowledged . 462 Appendix A Characteristic stress 463 The characteristic stress concept was proposed to establish a unified strength criterion 464 (simplified as β criterion) in a simple way [40], and has been applied to develop 3D elastoplastic 465 model for soil [19]. The expression of the characteristic stress t ensor can be written as: 466             1 1 2 2 3 3 3 12 ij r i j i j i j r r rc p n n n n n np p p                      (38) 467 where cij is the characteristic stress tensor whose principal directions coincide with those of the 468 ordinary stress tensor σij. A nd the expression of the characteristic principal stress rr()kkc p p 469 can be obtained the based on Eq. (38). σk (k =1, 2, 3) is the principal stress. pr is the reference stress 470 23 for dimensionless transformation, and pr=1kPa is commonly used. β is the characteristic stress index, 471 which is a non -dimensional material parameter that ranges from 0 to 1.  k in in Eq. (38) is the 472 corresponding directional cosine of the direction vector e(k) of the principal stress , that is, the 473 principal stress vector can be expressed as cke(k). Therefore, the spectral decomposition of the 474 characteristic stress tensor can be illustrated by the schematic diagram shown in Fig. 13. 475 476 Fig. 13 Spectral decomposition of the characteristic stress tensor . 477 Based on the friction rule on the octahedral plane under the action of the cij, a nonlinear unified 478 β strength criterion [40] was developed, which is used as the failure criterion in this paper : 479 s f ncFc (39) 480 where Ff is the characteristic failure stress ratio. The specific expressions of the characteristic stress 481 variables and parameters are listed in Tab. A -1, which are also com pared with the ordinary stress. 482 Tab. A -1 Comparison between the ordinary stress and the characteristic stress. 483 Ordinary stress Characteristic stress Stress tensor σij cij Principal stress , 1,2,3ii , 1,2,3i ir rc p ip  Mean stress p=σ ijδij /3 cn=cijδij /3 General shear stress   3 2ij ij ij ij q p p         s n n3 2ij ij ij ij c c c c c     Stress Lode angle  23 1 2 33arctan2    23 1 2 33arctan2cc c c c O 1x 2x 3x 1e 2e 3e (1)e (1) 1ec 1e 2e 3e (3)e (3) 3ec 1e 2e 3e (2)e (2) 2ec (1) 3n (1) 1n (1) 2n (3) 3n (3) 1n (3) 2n (2) 3n (2) 1n (2) 2n (1)e (3)e (2)e 24 Stress ratio η=q/p χ=c s/cn Phase t ransition stress ratio 6sin 3 sinM         1 sin 1 sin3 1 sin 2 1 sinF        Failure stress ratio f f f6sin 3 sinM         ff f ff1 sin 1 sin3 1 sin 2 1 sinF        Appendix B Specific expressions of fractional derivative 484 The specific expressions of the Caputo fractional derivative of a power function xn is given and 485 expressed as: 486     01 () 1n n n xnx xDxxn       (40) 487 where 0 xD means the Caputo derivative of order μ. 0 is chosen as the starting point and x is the 488 terminal point in the operator. Γ(·) indicates the Gamma function and its expression is 489  1 0dzze , which satisfies the recursive relation, i.e., Γ( x+1)= x·Γ(x). In addition, it is 490 worth to note that  0 0xDC , where C is a constant. The linear relationship of the Caputo 491 fractional derivatives is also necessary and can be expressed as: 492        0 1 1 2 2 1 0 1 2 0 2 x x x D A f x A f x A D f x A D f x       (41) 493 494