5 A non-orthogonal elastoplastic model for clay is proposed by combining the non- orthogonal plastic 6 flow rule with the critical state concept and presented from the perspective of the magnitude and 7 direction of the plastic strain increment . The magnitude is obtained based on the improved elliptical 8 yield function and the plastic volumetric strain dependent hardening parameter. The direction is 9 determined by applying the non -orthogonal plastic flow rule with the R iemann- Liouville fractional 10 derivative to the yield function without the necessity of additional plastic potential function . The 11 presented approach gives rise to a simple model with five parameters. All parameters have clear 12 physical meaning and can be easily identified by triaxial test s. The model performance is shown by 13 analyzing the evolution process of the yield surface, the hardening rule and the plastic flow direction . 14 The capability of the proposed model to capture the mechanical behavio urs of cl ay with different 15 stiffness is also confirmed by predicting test results from the l iterature. 16 Keywords : Clay ; Constitutive model; Elastoplasticity; Non-orthogonal plastic flow rule; Stress 17 dilatancy; Critical state 18 19 1 Ph.D. Candidate. Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, 20 Beijing University of Technology, Beijing 100124, China. E -mail: liangjingyu@emials.bjut.edu.cn . 21 2* Professor. Key Laboratory of Urban Security and Disaster Engineering of Ministry of Education, 22 Beijing University of Technology, Beijing 100124, China (C orresponding author). E -mail: 23 dechun@bjut.edu.cn . 24 3 Ph.D. Candidate. Key La boratory of Urban Security and Disaster Engineering of Ministry of Education, 25 Beijing University of Technology, Beijing 100124, China. E -mail: zhouxin@emials.bjut.edu.cn . 26 4 Professor. Key Laboratory of Urb an Security and Disaster Engineering of Ministry of Education, 27 Beijing University of Technology, Beijing 100124, China. E -mail: duxiuli@bjut.edu.cn . 28 5 Professor. Institute of Geotechnical Engineering , University of Natural Resources and Life Sciences, 29 Vienna , 1180, Austria. E -mail: wei.wu@boku.ac.at . 30 This is a pre-print of the paper published in Computers and Geotechnics, Volume 116, December 2019, 103200 DOI: https://doi.org/10.1016/j.compgeo.2019.103200 2 1 In troduction 31 The constitutive model can be used to determine the two basic elements of the strain increment 32 produced by the current stress increment, namely (i ) magnitude and (ii) direction . Within the 33 framework of plastic ity theory, the magnitude of the plastic strain increment can be determined by 34 the consistency condition together with the yield function and the hardening parameter . The 35 direction of the plastic strain increment, i.e., the plastic flow direction, is usually determined by the 36 orthogonal gradient of the plastic potential surface . If the plastic potential function is chosen to be 37 the same with the yield function, the associated flow rule is adopted [1-4], in which the plastic flow 38 direction is orthogonal to the yield surface. But the associated flow rule has been reported to be 39 unsuitable for describing the mechanical behaviours of geomaterials [5-7]. I f the plastic potential 40 function is not identical to the yield function, the non -associated flow rule is adopt ed [8-11], and the 41 plastic flow direction perpendicular to the plastic potential surface is non -orthogonal to the yield 42 surface. It is also called the non -orthogonal conditions, but the non -orthogonali ty with respect to the 43 yield surface is accomplished by constructing a new plastic potential function [12-15]. However, it is 44 difficult to construct a required plastic potential function , especially for geomaterials with pressure 45 sensitivity and dilatancy [16-19]. In addition, some additional parameters without clear physical 46 meanings might be needed. 47 The direction of plastic strain increment can also be directly determined by the non -orthogonal 48 gradient of the yield function. Therefore, the plastic potential function is no longer necessary for 49 determining the plastic flow direction . The establishing process of the elastoplastic model becomes 50 simpler and few material parameters are required . The fractional partial derivative possesses the 51 ability to adjust the non-orthogonal gradient direction of the differentiable function [20, 21]. In order 52 to determine the general tensor form of the plastic strain increment based on the commonly used 53 yield function for soil, a general form of non -orthogonal plastic flow rule [22, 23] is proposed based 54 on the fractional partial derivative a nd the covariant transformation. The proposed non -orthogonal 55 plastic flow rule has been applied to establish elastoplasic models for soi l [22] and concrete [23]. 56 Nevertheless, the non- orthogonal plastic flow rule has not been directly introduced into the critical 57 state theory . Their com bination is promising for reproducing the mechanical behaviour of soil. 58 During establishing non -orthogonal elastoplastic models for soil [21, 22], the plastic flow direction is 59 expressed by the C aputo fractional derivative. However, t he Caputo fractional derivative of constant 60 3 is equal to zero, which makes it unsuitable for the ellipitical yield function of the modified Cam - 61 clay model when the fractional order is greater than one. 62 In this paper, a framework of non-or thogonal elastoplastic critical state model for clay is 63 presented . A general non -orthogonal plastic flow rule with the Riemann- Liouville (RL) fractional 64 derivative is developed to determine the plastic flow direction without recourse t o t he plastic 65 potential function. The proposed model is capable of reasonably capturing the mechanical behaviour 66 of clay with different stiffness. 67 2 S tress -strain relationship of the non -orthogonal model for clay 68 In what follows, the main ingredients of the non-or thogonal elastoplastic model are 69 recapitulated. The strain increment is decomposed into two components: the reversible elastic strain 70 increment edij and the irreversible plastic strain increment pdij : 71 epd d dij ij ij   (1) 72 where edij can be determined by the Hooke’s law and the pdij can be determined based on the 73 non-orthogonal plastic flow rule in this paper . 74 2.1 Elastic strain increment 75 The elastic strain increment can be obtained from the Hooke's law , in which t he principal axis 76 direction of the elastic strain increment is coincide with that of stress increment and the magnitude 77 of the elastic strain increment can be determined by the elastic parameters as follows: 78 e1d d dij ij kk ijEE    (2) 79 where ν is the Poisson’s ratio. E is the Young’s modulus , which can be determined from the linear 80 relationship between void ratio e and the natural logarithm of mean stress p, i.e., e-lnp relationship, 81 as shown in Fig. 1. The expression of E in Hooke’s law is : 82  013 1 2eEp (3) 83 where e0 is the initial void ratio. κ is the swelling index of the swelling line (SWL: e=ΠS–κ· lnp, i n 84 which ΠS is the void ratio at p=1 k Pa on the SWL) , as shown in Fig. 1. The symbol λ in Fig. 1 is the 85 compression index of the normally consolidated line (NCL: e=Π N–λ· lnp, where Π N i s the void ratio 86 at p=1 kPa on the NCL). The determined method of E is the same as that in the critical state theory 87 [1]. Apparently, the modulus E in Eq. (3) i s variable and depends on p, which reflects t he pressure 88 4 dependence of clay . 89 90 Fig. 1 e–lnp relationship under isotropic consolidation conditions. 91 2.2 Plastic strain increment 92 The plastic strain increment can be expressed as: 93 pdij ijm (4) 94 where Λ is the plastic multiplier that represent s the magnitude of d𝜀𝑖𝑗p. Λ can be determined by 95 submitting the yield function and the hardening parameter into the consistency condition. mij 96 represents the direction of d𝜀𝑖𝑗p, which can be determined by the non -or thogonal plastic flow rule. 97 2.2.1 Magnitude of plastic strain increment 98 Withi n the framework of plastic ity theory, yield function and hardening parameter are the main 99 parts to determine the magnitude of plastic strain increment. Yield surface is the equivalent surface 100 of the hardening parameter , and its general form is: 101  p,ij kl f f H (5) 102 where H is the hardening parameter , which is the function of specific plastic strain. 103 The consistency condition ensures the consistency of the stress state and the hardening 104 parameter, which makes the yield function with the hardening parameter value always pass through 105 the current stress state point. It is used to close the set of constitu tive equations and determine the 106 plastic multiplier Λ. The consistency condition states that d f=0 at yield, which is the total dif ferential 107 of the yield function f in Eq. (5). Therefore, for the isotropic hardening rule, the consistency 108 condition can be expressed as: 109 p pd d d 0ij kl ij klf f HfH       (6) 110 0A xA 0lnp lnp e 0e xe e pe ee   lnp 0A NCL SWL 5 By submitting Eq. (4) i nto Eq. (6), the plastic multiplier Λ can be obtained: 111 pdij ij kl klf fHmH    (7) 112 Once the yield function and the hardening parameter is decided, the expression of Λ can be 113 accordingly obtained. 114 In this paper, an improved elliptical yield surface , which has the same form and convexity with 115 the MCC yield function [1], is adopted and expressed as : 116  2 2 20x f q N p p p    (8) 117 where px is the equivalent consolidation stress , which is the intersection coordinate of the yield curve 118 and the p-axis and can be expressed by any stress state ( p, q) on the yield surface, i.e., 119 22 xp q N p p . N is the axial ratio of the vertical and horizontal axe s of the elliptic yield curve . 120 If N is forced to be the same with the critical state stress ratio M, Eq. (8) becom es the yield function 121 of the MCC model . If N does not equal M, there will be two possible scenarios for the improved 122 yield surface based on the relative relationship between N and M: (a) N<M; (b) N’>M, where N’ is 123 used as the symbol for the condition N’>M. The corresponding shapes of yield curve a re illustrated 124 in Fig. 2. 125 126 (a) (b) 127 Fig. 2 Possible shape s of the yield curve when : (a) N<M; (b) N’>M. 128 129 An experimentally validated form for the isotropic hardening parameter , which is the same 130 with the MCC model [1], is used in this paper to reflect the strain hardening feature of clay . B ased 131 q p O xp N f 2xp M 1 q p O xp 'f 2xp N M 1 6 on the e –lnp relationship under isotropic compression conditions as shown in Fig. 1, the hardening 132 parameter H can be expressed as: 133 p vp 0lnxpHlp   (9) 134 where lp=(λ-κ)/(1+ e0). px connects the hardening parameter 𝜀vp with the yield function, and its 135 specific expression associated with 𝜀vp can be derived from Eq. (9): 136 p v 0 pexpxppl  (10) 137 With the expression s of the yield function in Eq. (8) a nd the volumetric strain dependent 138 hardening parameter in Eq. (9) t ogether with its incremental expression d𝜀vp=Λ·𝑚p from Eq. 139 (4), the expression of Λ in Eq. (7) can be converted accordingly as follows: 140 p vdd pffpqpq fm  (11) 141 The partial derivatives of the yield function f in Eq. (11) can be derived from Eqs. (8) and 142 (9). Theref ore , Λ can be e xpresse d as: 143  2 p 22 d 2 dx xpN p p p q qlN pp m (12) 144 By substituting px=q2/(N2p)+p from Eq. (8) i nto Eq. (12), 145    2 2 2 p 2 2 2d 2 d pN p q p pq q l p q N p m   (13) 146 Particularly, as indicate d by Eq. (13), mp as sociated with the volumetric strain dependent hardening 147 parameter affects the calculation of Λ, which will be discussed in the below section . 148 2.2.2 Direction of plastic st rain increment 149 The direction of plastic strain increment is determined by the non -orthogonal plastic flow rule , 150 which is the non-orthogonal gradient of the yield surface . By com bining the fractional derivative 151 and the covariant transformation, Lu et al. [22, 23] proposed a non-or thogonal plastic flow rule in 152 general form, which can be expressed as: 153 pdkl ij ij kl ijSfmS   (14) 154 7 where Skl is the variable that construct s the yield function, which can be the stress invariant s of σij, 155 like the commonly used (p, q, θ). The plastic flow direction mij is actually t he non-orthogonal 156 gradient of the yield surface , which is determined by the fractional order partial der ivative 157 kl fS in Skl space firstly and then transformed into σij space by the covariant transformation 158 Skl/σij. Actually, t he gener al tensor form of mij in Eq. (14) i s the kernel for deriving the general 159 expression of d𝜀𝑖𝑗p that is essential for model applicat ion and the numerical research. 160 The geometric schematic diagram of the non -orthogonal plastic flow rule expressed based on 161 the fractional derivative is shown in Fig. 3. T he components ( m1, m2, m3) of the plastic strain 162 increment dire ction can be determined by the non -orthogonal plastic flow rule. The resultant 163 direction vector m indicated by the red arrow in Fig. 3(a) is the fractional order gradient of the yield 164 function , which is not orthogonal to the yield su rface. By contrast, the v ector n indicated by the 165 black arrow is the gradient of the yield surface , i.e., m when μ=1, which is orthogonal to the yield 166 surface, as shown in Fig. 3(b). 167 168 (a) (b) 169 Fig. 3 Comparison of the fractional gradient vector when μ≠1 and the normal vector when μ=1. 170 171 The RL operator of the fractional partial derivative is used in this paper to derive the explicit 172 solution of the fractional derivative of Eq. (8). The specific expression s of the RL fractional 173 derivative for a power function and a constant as well as the linear relation of the fractional 174 derivative are required and listed as follows: 175     01 () 1n nn xn xD x xxn        (15) 176 O 1 2 3 m 3m 1m A 2m   Yield surface: ijf O 1 2 3 A   Yield surface: ijf n 3n 1n 2n 8   0() 1xC CxDCx     (16) 177        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       (17) 178 where 0Dμ x is the operator of the fractional derivative of order μ, and 0 is chosen as the starting point 179 of the fractional derivative. Γ(𝑧)=∫𝑒−𝜏𝜏𝑧−1d𝜏∞ 0,(Re(𝑧)> 1) is the Euler Gamma function . 180 Based on Eq. (14), the non -orthogonal plastic flow direction mij can be expr essed as : 181 ij ij ijf p f qmpq        (18) 182 where the fractional partial derivative of the yield function in Eq. (8) w ith respect to p and q can 183 be expressed as (mp, mq)= ( fp , fq ), which is determined by applying the RL operator 184 to Eq. (8). B ased on Eqs. (15)-(17), px=q2/(N2p)+p from Eq. (8) an d the recursive propert y of the 185 gamma function, i.e., Γ( x+1)= x·Γ(x), (mp, mq) can be thus obtained : 186        2 2 2 2 2 223 ,,33pqq N p q mmp              (19) 187 The covariant transformation coefficient betwe en (p, q) and t he ordinary stress tensor σij can be 188 expressed as: 189  1 3 3 2ij ij ij ij ijp p q q    (20) 190 By substituting Eqs. (19) a nd (20) i nto Eq. (18), the explicit expression of the plastic strain 191 increment direction is expr essed into the tensor form : 192        2 2 2 2 2 223 3 1 3 3 3 2ij ij ij ijq N p q p mpq               (21) 193 By submitting Eq. (13) and Eq. (21) into Eq. (4), t he plastic strain increment is obtained: 194         2 2 2 1 p 2 2 2 2 2 2 0d 2 d 33 1d1 3 2 2 2ij ij ij ijN p q p pq q p q p e N p q p q N p            (22) 195 Specifically, when μ=1 and N=M, Eq. (22) can also degenerate into the explicit expression of the 196 MCC model : 197 9     2 2 2 p 2 2 2 2 2 2 0d 2 d 3 d13ij ij ij ijM p q p pq q pp e M p q p q M p        (23) 198 Finally , the stress -strain relationship for the non-or thogonal critical state elastoplastic 199 constitutive model is obtained by submit ting Eqs. (2) and (22) i nto Eq. (1). 200 2.3 Parameter calibration 201 There are only five material parameters with clear physical meaning in the proposed non- 202 orthogonal critical state model, among which four material parameters (λ , κ, M and ν) are identical 203 to the MCC model . The critical state stress ratio M can be represented by the internal friction angel 204 φ under triaxial compression conditions : 205 6sin 3 sinM  (24) 206 Only the critical void ratio change ΔΠ= ΠN−ΠC is added, which is the vertical distance between 207 the NCL and the CSL , as shown in Fig. 4. U nder triaxial compression conditions from A to B shown 208 in Fig. 4, the change of the void ratio Δep is exactly equal to ΔΠ. I t also equals the void ratio variation 209 along the path AAxB, as shown in Fig. 4. The void ratio v ariations corresponding to the path AAx and 210 the path AxB are Δe and Δee, respectively. And Δe=λ·ln(px/p0) and Δee=κ·ln(px/p0) are easily obtained 211 based on the expressions of the NCL and the SWL. Therefore, Δep=Δe–Δee can be expressed as Δep 212 =ΔΠ=(λ-κ)·lnξ, where ξ=px/p0. Once the parameters λ and κ are determined, ξ and ΔΠ are one -to- 213 one correspond ing parameters . 214 exp (25) 215 216 Fig. 4 Vertical distance between the NCL and the CSL. 217 218 The fractional order μ and the axial ratio N of the yield curve are two model parameters used 219 ee 0lnp lnp e N  lnxp A C B  ln1  p= ln e       e pe xA  NCL CSL SWL 10 in this paper , which are not independent and can be determined by ξ. According to the consistency 220 condition, the yield curve finally passes through the critical state stress point ( pcs, qcs) = ( p0, M·p0). 221 By submitting (p0, M·p0) into Eq. (8) and together with ξ =px/p0 , N can b e determined : 222 1MN   (26) 223 In particular, ξ=2 can be determined from the expressions of the MCC yield function , and ΔΠ=( λ- 224 κ)·ln2 is settled accordingly . The fractional order μ can be determined based on the stress -dilatancy 225 ratio d equal s 0 at the critical stress state. The expression of d is presented in details in section 3.2. 226 And thus, by combining Eq. (26) w ith d=0 in Eq. (28), one can get μ as follows : 227 121 (27) 228 From Eqs. (26) an d (27), we can see that μ and N are not independent but coordinated with 229 the critical state concept . Only one parameter ξ or ΔΠ is needed . When ξ>2, μ>1 and N<M can be 230 obtained, which reflects that the clay is soft with the stiffness lower than that of the MCC model . 231 When ξ<2, μ<1 and N>M can be obtained, which reflects that the clay is hard with the stiffness 232 greater than that of the MCC model . C ompared with the MCC model, by adding only one parameter 233 ξ (or ΔΠ), a non-orthogonal critical state model can be easily developed and can be used to describ e 234 behaviours of clay with different stiff ness. 235 3 Performance of the model 236 The performance of the non -orthogonal elastoplastic model can be visualized by analyzing the 237 evolution of the yield surface, the hardening rule and plastic flow direction . By matching and 238 coordinating the yield curve shape , the hardening parameter and the non-or thogonal plastic flow 239 rule, the proposed model can be used to capture the behaviours of clay with different stiffness . I t 240 can also be used to predict existed test results in the literature. 241 3.1 Analysis on the evolution of yield curve and hardening rule 242 Parameters for the new Fujinomori clay (NF -clay) [24] as listed in Tab. 1 i n section 3.4 a re used 243 for demonstration. Two ξ-values: ξ=2.3 and ξ’=1.7, are used , which correspond to N=1.19<M=1.36 244 and N’=1.63>M for the yield curves , as listed in Fig. 5(a) and (b), respectively. U nder the same 245 stress step BiBi+1 on triaxial compression stress path with p=196 kPa, the ev olution process of the 246 yield curves and the hardening rules for two ξ-value s are presented in Fig. 5(a) and (b), r espectively . 247 11 As the stress state changes from Bi to Bi+1, as shown in Fig. 5(a), the yield curve changes from the 248 yield curve fi to the next yield curve fi+1, and the equivalent consolidated pressure point changes 249 from Ai to Ai+1 accordingly . Since the yield surface is the equivalent surface of the hardening 250 parameter p v, Ai and Bi (or Ai+1 and Bi+1) on the same yield curve correspond to the same value of 251 p v , that is, the difference of p v from Bi to Bi+1 is just the same with that from Ai to Ai+1. The plastic 252 deformation p v produced by the stress step BiBi+1 can be thus calculated by the change of the 253 equivalent consolidation pressure from Ai to Ai+1 based on the NCL and the SWL. Based on the 254 evolution of the yield surface and the hardening parameter, the magnitude of the plastic strain 255 increment under the same stress step can thus be calculated . The similar evolution process for ξ’=1.7 256 is show in Fig. 5(b). However, i t is obvious that t he difference of the equivalent consolidated 257 pressure between Ai and Ai+1 in Fig. 5(a), i.e., Δ px= px,i+1–px,i, is bigger than that between A’i and 258 A’i+1 in Fig. 5(b). According to Eq. (9), the change of the hardening parameter for N=1.19< M is 259 greater than that for N=1.63 >M, which leads to a larger magnitude of the plastic strain increment 260 for N<M. 261 1N 0 100 200 300 4000100200300400 0 100 200 300 4000.620.640.660.680.70 Hardening rule ξ =2.3 (ΠC=1.02)SWLΔep ei+1ei CSLBi+1 Ai+ 1Bi Ai p (kPa)Void ratio e NCL: e=ΠN- λ·lnpΠ'C=1.04 ΠC=1.02 Parameters used for prediction N =1.19 < M =1.36 fi+1 fiBi+1 Bi Ai+1Ai p (kPa)q (kPa)M 0 100 200 300 4000100200300400 0 100 200 300 4000.620.640.660.680.70 Hardening rule ξ' =1.7 (Π 'C=1.04)SWLΔep' ei'+1ei' NCL: e=ΠN- λ·lnp CSLBi'+1 Ai'+ 1Bi' A'i p (kPa)Void ratio eΠ'C=1.04 ΠC=1.02 Parameters used for prediction N' =1.63 > M =1.36 fi'+ 1 fi'Bi+1 Ai'+ 1Bi A'iM p (kPa)q (kPa) N'1 262 12 (a) (b) 263 Fig. 5 Evolution of yield curve and hardening rule for: (a) ξ=2.3; (b) ξ’=1.7. 264 265 3.2 Analysis on the non-ort hogonal plastic flow direction 266 The non -orthogonal plastic flow rule actually plays a key r ole in the coordination of the yield 267 surface with the critical state concept. Apparently, the plastic flow direction expressed by Eq. (21) 268 is affected by two model parameters μ and N. According to Eqs. (26) and (27), t wo configurations: 269 (a) μ=1.23 and N=1.19 when ξ=2.3, and ( b) μ’=0.57 and N’=1.63 when ξ’=1.7, are determined . The 270 triaxial stress pathes A0B at p=196 kPa are the same for these two configurations. The plastic flow 271 direction under the path A0B can be obtained with the determined parameters , as shown in Fig. 6. 272 The plastic strain increment direction, as indicated by the red arrow for configuration (a) in Fig. 6(a) 273 and the blue arrow for configuration (b) in Fig. 6(b), changes from horizontal to vertical as the stress 274 ratio η=q/p increases from 0 to M . For the condition of η=0 at point A0, t he plastic flow direction is 275 horizontal to the right, i.e., d𝜀sp= 0, and is independent of μ. With the increase of η, the plastic 276 flow direction changes from horizontal to vertical. When η reaches M at point B, the plastic strain 277 increment is vertically upward, i.e., d𝜀vp= 0, and the critical state is reached, but it is not necessary 278 to be the peak point of the yield curve . For the configuration (a), as shown in Fig. 6(a), the stress 279 state point passes through the pea k point of the yield curve and finally reaches the critical state point 280 B that locates on the left of the peak point. For the configuration (b), as shown in F ig. 6(b), the stress 281 state point does not reach the peak point of the yield curve, and t he critical state point B locates on 282 the right of the peak point . The comparison of the plastic flow direction s for these two configurations 283 is shown in Fig. 7. The orthogonal gradient of the yield curves at point B is also shown by black 284 dashed vectors in Fig. 7 t o further illustrate that the critical state point does not necessary to be the 285 peak point of the yield curves. 286 13 1 0 90 180 270 360 450 5400100200300400500 AxB A0 p (kPa)q (kPa)M Yield curves Peak point Stress pathPlastic flow direction μ=1.23 ( N=1.20< M) N 0 90 180 270 360 450 5400100200300400500 Yield curves Peak point Stress pathPlastic flow direction μ' =0.57 ( N'=1.63> M) A'xB A0M p (kPa)q (kPa) 1 N' 287 (a) (b) 288 Fig. 6 Plastic flow direction with fractional -order plastic flow rule for two configurations : (a) μ=1.23 289 and N=1.19 when ξ=2.3; (b) μ’=0.57 and N’=1.63 when ξ’=1.7. 290 Plastic flow direction μ =1.23 ( N =1.19< M) μ'=0.57 (N' =1.63< M) μ =1.00 0 100 200 300 400 5000100200300400500 BiB A0 p (kPa)q (kPa)M 291 Fig. 7 Comparison of the p lastic flow direction s for two configurations . 292 293 The variation rule of the plastic strain increment can also be reflected by the stress -dilatancy 294 relationship d-η, w hich can be determine d based on Eq. (21) and expressed as: 295    2 2 2p v p 22 s2 d=d 3N d         (28) 296 For the two possible configurations , the stress -dilatancy curves are indicated by red curve and blue 297 curve, respectively . The stress- dilatancy curves for configurations (a) and (b) are below and above 298 the MCC stress -dilatancy curve indicated by the black dashed curve. And eventually, the stress state 299 on the stress -dilatancy curves reaches the same critical state point where η reaches M, as shown in 300 Fig. 8. 301 14 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.00.00.30.60.91.21.51.8 MDilatancy relationship at p = 196kPa μ = 0.57 ( N'=1.63 > M) μ = 1.00 (MCC model, M=1.36) μ = 1.23 ( N=1.19 < M) η -dp v/dp s 302 Fig. 8 Stress -dilatancy curves. 303 304 3.3 Description of soil behaviour with different stiffness 305 For the condition s that μ=0.50, 0.7 5, 1.00, 1.25 and 1.5 0, the corresponding parameters N=1.67 , 306 1.52, 1.18, 1.17 and 0.96, ξ=1.67, 1.80 , 2.00, 2.3 3 and 3.00 ar e liste d in Fig. 9. U nder the same stress 307 path A0B at constant p, the yield curves change from the initial ones to the terminal ones that pass 308 through the critical state point B , as shown in Fig. 9(a) . It is obvious that the critical state point B is 309 not the peak point of the yield curves except for the condition of μ=1.00. For different yield curves, 310 the equival ent consolidation pressures are px1, px2, px3, px4 and px5, respectively. These equivalent 311 consolidation pressures can be used to determine the position of CSLs based on the NCL and the 312 SWL, as shown in Fig. 9(b), w hich result in a difference in the change of the void ratio for the stress 313 path A 0B. 314 0 100 200 300 400 500 6000100200300400500 px3px2 μ=0.50 ( N=1.67> M) μ=0.75 ( N=1.52> M) μ=1.00 ( N=1.36= M) μ=1.25 ( N=1.17< M) μ=1.50 ( N=0.96< M) B A0 p (kPa)q (kPa)M p0px1px4px5 4.8 5.2 5.6 6.0 6.4 6.80.520.560.600.640.680.72 p0 px5px4px3 px2px1SWLHardening rule ξ=1.67 ξ=1.80 ξ=2.00 ξ=2.33 ξ=3.00A0 ln pVoid ratio eNCL: e=ΠN- λ·lnp CSLCSL : e=ΠC- λ·lnp 315 (a) (b) 316 Fig. 9 Coordination of the critical state with: (a) the yield curve and (b) the hardening rule . 317 318 By further introducing the non -orthogonal plastic flow rule presented in section 3.2, the 319 15 behaviours of clay with different stiffness are thus captured and shown by the stress -strain curves in 320 Fig. 10. W ith the increase of μ-value, the predicted principal strains and the plastic volumetric strain 321 at the same stress increment is larger. The prediction result is identical to that of the MCC model 322 when μ is 1. That is, the greater the value of μ, the softer property of clay can be captured. With the 323 change of μ, the behaviours of clay with different stiffness can be reasonably described by the 324 proposed non -orthogonal critical state elastoplastic model. 325 Predictions ( μ≥1.00) μ=1.00 μ=1.25 μ=1.50Predictions ( μ≤1.00) μ=0.50 μ=0.75 μ=1.00 -10 -5 0 5 10 15 01234 0 42ε2=ε3 ε1 ε1 /% εvε3 /% v /%1/3 326 Fig. 10 Stre ss -strain relationship with different value of μ. 327 328 3.4 Model pred ictions 329 To assess the performance of the non -orthogonal critical state model for clay, experimental 330 data from literature are used. The test data of the saturated remoulded normally consolidated NF- 331 clay [24] and the Fujinomori cl ay (F-clay) [25] under triaxial compression condition s are used f or the 332 verification of the non -orthogonal critical state model . Material p arameters listed in Tab. 1 from [15] 333 are adopted for model prediction . 334 Tab. 1 Material p arameters used for model prediction 335 Material p arameters λ/(1+e0) κ/(1+e0) ν φ ξ NF-clay 0.0444 0.0047 0.3 33.7° 2.3 F-clay 0.0508 0.0112 0.3 35° 2.3 Drained triaxial compression tests with p=196 kPa on the NF -clay [24] and the F -clay [25] were 336 conducted to explore the strength and deformation properties of saturated remoulded normally 337 consolidated clay. For the NF -clay, test data and material parameters used for prediction are from 338 Chowdhury and Nakai [24]. In this paper, ξ =2.3 can be determined from Eq. (25) based on the void 339 16 ratio variation of the NF -clay. On the other hand, the MCC model (ξ=2.0) is also used for 340 comparison and demonstration. The fractional order μ =1.23 and 1.00 are determined in accordance 341 with the value of ξ=2.3 and ξ=2, respectively. A nd the axial ratio N=1.19 and 1.36 (= M) can also be 342 determined accordingly based on ξ and M. Model predictions and comparisons with test results of 343 the NF -clay indicated by the red open circles can be found in Fig. 11. Apparently, model predictions 344 with the fractional order μ=1.23 that represented by the red solid line capture the test results better 345 than predictions of the MCC model (μ =1.00) represented by the black dash line. For the F -clay, the 346 test dat a and material parameters are from Lu et al. [15]. Similar variation rule can be found in Fig. 347 12. As shown in Fig. 11 and Fig. 12, the gr eater the μ value, the greater the principal strains occurring 348 at the same principal stress ratio . Comparisons and analyses indicate that the non- orthogonal model 349 can be used to well describe the deformation properties of clay with different st iffness. 350 Predictions μ=1.23 ( N=1.20< M) μ=1.00 (MCC model)Fujinomori clay Test data -10 -5 0 5 10 15 01234 0 42ε2=ε3ε1 ε1 /% εvε3 /% v /%1/3 p=196 kPa 351 Fig. 11 Model p redictions for test results of the NF -clay (data from [24]) under triaxial compression s. 352 Predictions μ=1.23 (N=1.24< M) μ=1.00 (MCC model )Fujinomori clay Test data -10 -5 0 5 10 15 01234 εvε2=ε3ε1 2 40 ε3 /% v /%ε1 /%1/3 p=196 kPa 353 Fig. 12 Mo del p redictions for test results of the F -clay (data from [25]) under triaxial compression s. 354 355 4 Conclusions 356 17 In this paper, a non -orthogonal elastoplastic model with critical state is presented , and the 357 magnitude and direction of plastic strain increment are directly determined by the yield function. 358 The proposed non-orthogonal model possesses simple formulation s. Only five material parameters 359 with clear physical meanings are needed . They all can be de termined by the traditional laboratory 360 tests. The proposed non-orthogonal model can degrade into an orthogonal model, which is the MCC 361 model, when the fractional order μ equals 1. The critical state theory is improved to capture t he 362 mechanical behaviour s of clay with different stiffness. Moreover, it remains exploratory to be 363 improved to furnish modelling dilatancy of sand. 364 Acknowledgements 365 This work was supported by the National Key Research and Development Program of China 366 (Grant No. 2018YFC1504301 ), the National Natural Science Foundatio n of China (Grant Nos. 367 51778026 and 51421005). The financial support from the European Commission under the projects 368 HERCULES (Grant No . 778360) and RE -BUILT (Rehabilitation of the Built Environment in the 369 Context of Smart City and Sustainable Development Concepts for Knowledge Transfer and Lifelong 370 Learning ) is gratefully acknowledged. 371 372