不使用状态变量和塑性势的土壤状态相关非关联行为的分数阶建模

摘要

研究发现，颗粒土的本构行为取决于其密度和压力（即材料状态）。为了捕捉这种状态依赖性，根据经验提出了各种状态变量，并将其引入到现有的塑性势函数中，这不可避免地导致了一些模型参数的复杂性和无意义性。本研究的目的是在不使用预先定义的塑性势和状态变量的情况下，从理论上研究颗粒土的状态相关非关联行为。通过将应力分数算子引入边界面塑性，建立了一种新的颗粒土状态相关非关联模型。与以往使用经验状态变量的研究不同，本研究中的土壤状态和非结合性是通过解析解考虑的，在不使用附加塑性势的情况下，可以获得状态相关的塑性流动规律和相应的硬化模量。还讨论了与一个众所周知的经验状态变量的可能数学联系。发现塑性流动与加载方向以及材料硬化之间的非关联性受分数阶导数控制。为了验证所提方法，对不同颗粒土的一系列排水和不排水三轴试验结果进行了模拟和比较，结果表明，模型预测与相应的试验结果吻合良好。

1介绍

众所周知，砂石等粒状土的强度和变形特性在很大程度上取决于其密度和压力（材料状态）[1]. 在本构建模过程中适当考虑状态相关性之前，通常需要不同的模型参数来模拟具有不同初始密度或承受不同围压的粒状土的应力应变行为[2,三,4].

为了更好地理解和统一粒状土的本构模型，提出了各种不同的经验状态变量，例如阈值和当前孔隙比之间的差异与阈值和临界孔隙比之间差异的状态比[5]、电流与临界孔隙比之比[6]扰动状态下的扰动概念[7]，电流与临界平均有效应力的应力比[8]和使用最广泛的状态变量(ψ)由当前和临界状态孔隙比之间的差异定义[9]. 可以发现，近年来，制定合理的状态相关非关联塑性流动规则（或应力-微分关系）至关重要。一种流行的方法是修改现有的应力-微分方程，例如Cam-clay（CC）应力-积分方程[10]和Rowe的应力-微分方程[11]，通过经验整合ψ[12,13,14,15,16]. 不可否认，这种方法总是可以大大提高模型性能；然而，状态变量与现有本构参数之间的经验相关性必然会导致更多的模型参数。研究发现，分数力学是捕捉弛豫的有效方法[17]，扩散[18,19,20]，和应力应变[14,21]材料性能。为了在不损失建模能力的情况下减少模型参数的数量，Sun和Shen[21]通过简单地对屈服面进行分数阶导数，提出了颗粒土的非关联塑性流动规律，其中获得的矢量（塑性流动方向）不再垂直于屈服面，即使不使用额外的塑性势。这种非正态性随着分数阶的增加而增加(α)减少[14,22,23]. 为了考虑状态依赖性，提出了状态依赖分数塑性模型[14]通过经验整合ψ然而，这使得获得的应力-微分方程的参数缺乏物理意义。

本研究试图从理论上研究颗粒土的状态相关非相关应力应变行为。利用严格的数学方法，建立了一个不使用任何经验状态变量和塑性势的状态相关非关联本构模型。本文提出了状态相关塑性流动规则和相关硬化规则的分析推导，而不是通过经验合并状态变量来模拟土壤状态的相关性。由于分数导数是以积分形式定义的，土壤状态是通过从下限（当前应力状态）到上限（临界应力状态”）的积分范围来捕获的。本文主要分为四个部分：第四部分。 2定义了本研究中使用的基本本构关系和相关分数导数；第节。 三在不使用塑性势的情况下，发展了一种新的状态相关本构模型，其中状态相关分数塑性流动规则是解析推导出来的；第节。 4给出了模型参数的辨识和灵敏度分析；第节。 5根据不同颗粒土的一系列实验室试验结果进行模型验证；第节。 6研究结束。为了简单起见，本研究中的所有推导和讨论均限于均质和各向同性材料。

2符号和定义

2.1组织关系

根据土力学中的传统假设，压缩应力和应变被认为是正的，而广义应力和变形被认为是负的。本研究中提到的所有应力都是利用Terzhgi有效应力原理得到的有效应力。在弹塑性模型中，有四个主要部分，即弹性刚度张量(电子)，塑性载荷张量(米)，塑性流动张量(n个)和硬化模量(H（H）)。因此，总应变(\（瓦雷普西隆{ij}）)可以分解为以下弹性(\（\varepsilon_{ij}^{e}\）)和塑料(\（\varepsilon_{ij}^{p}\）)部件：

$$\Delta\varepsilon\｛ij｝=\Delta\varepsilon\｛ij｝^｛e｝+\Delta\varepsilon\｛ij｝^｛p｝$$

(1)

哪里\（i，j=1，2，3）Δ表示增量，上标e（电子）和第页分别表示弹性和塑性组件。基于胡克定律，增量弹性应变张量(\（Delta\varepsilon_{ij}^{e}\）)可以与增量有效应力张量相关(\（Delta\sigma'_{ij}\）)由

$$\Delta\sigma'_{ij}=E_{ijkl}\Delta\varepsilon_{kl}^{E}$$

(2)

哪里\（E_{\mathit{ijkl}}\）表示四阶弹性刚度张量，其定义如下[24]:

$$E_{ijkl}=（K-2G/3）\δ_{ij}\δ_}+G（δ_{ik}\δ_{jl}+\δ__{il}\△_{jk}）$$

(3)

在哪儿\（δ{ij}）是Kronecker三角洲。K（K）和G公司分别是体积模量和剪切模量，可以用[25]

$$\begin{aligned}&K=\frac{1+e}{\kappa}p'，\end{aligned}$$

(4)

$$\begin｛aligned｝&G=\frac｛3（1-2）｝｛2+2｝K，\end｛aligned｝$$

（5）

哪里κ是膨胀线的梯度\（e-\ln p'）平面；e（电子）是样品的当前孔隙比；ν是泊松比。\（p'=\sigma'{ij}\delta{ij{/3）是平均有效主应力。此外，广义剪切应力\（q=平方{3/2s_{ij}秒_｛ij｝｝\），其中\（s{ij}=\sigma'{ij{-\sigma'{kk}\delta{ijneneneep/3），是偏应力张量。相应的体积应变(\（瓦雷普西隆{v}）)可以定义为\（\varepsilon{v}=\varepsilon{ij}\delta{ij{），而广义剪切应变\（\varepsilon_{s}=\sqrt{2/3e_{ij}电子_｛ij｝｝\），其中\（e_{ij}=\varepsilon_{ij{-\varepsilon_{v}\delta_{ij}\），是偏应变张量。

在塑性方法中，增量塑性应变张量(\（Delta\varepsilon_{ij}^{p}\）)可以与增量有效应力张量相关，从而[26]

$$\Delta\varepsilon_{ij}^{p}=\frac{1}{H} n个_{ij}米_{kl}\Delta\sigma'_{klneneneep$$

(6)

其中塑性荷载张量(\（m _｛ij｝\）)垂直于屈服面，因此可以通过对屈服面进行一阶导数来确定；塑性流动张量(\（n _｛ij｝\）)与非关联岩土材料的屈服面不垂直，可以通过对屈服面进行分数阶导数来确定，如Sumelka和Nowak所示[23]和Sun等人[27]; 硬化模量(H（H）)对于粒状土，涉及尺寸硬化和屈服面相对于边界面的相对位置，稍后将进行定义。组合公式(2)–(6)，可给出以下弹塑性本构关系：

$$\Delta\sigma'_{ij}=\biggl[E_{ijkl}-\frac{E_{ijct}n个_{ct}电子_{klrs}米{rs}}{H+m_{ct}电子_{ctab}个_{ab}}\biggr]\Delta\varepsilon_{kl}$$

(7)

2.2分数导数与屈服面

分数导数有许多不同的定义[17,18,19,20,28,29]，其中每一个都是积分形式，并且有点复杂。因此，只有一些函数，例如幂律函数，有一定的解析解，而其余的函数需要数值近似[30,31]. 遵循分数塑性[14,22,32]以下是著名的卡普托左侧（等式(8))和右侧（等式(9))分数导数[33,34]在本研究中使用：

$$开始{aligned}&{}{{sigma'{c}}D_{\sigma'}^{\alpha}f\bigl（\sigma'\bigr）=\frac{1}{\varGamma（n-\alpha）}\int_{\sigma'_{c}}^{\sigma'}\frac}f^{（n）}（\chi），D\chi}{（\sigma'\chi>\sigma'_{c}，\end{aligned}$$

(8)

$$开始{对齐}&{}_{\sigma'}D_{\sigma'_{c}}^{\alpha}f\bigl（\sigma'\bigr）=\frac{（-1）^{n}}{\varGamma（n-\alpha）}\int_{\σa'}^{\ sigma'{c}{\frac}f^{（n）}（\chi），D\chi}{（\higma'）^{\阿尔pha+1-n}}，\quad\sigma'{c}>\sigma'，\end{aligned}$$

(9)

哪里D类(\（=\partial^{\alpha}/\ partial\sigma^{\prime\，\ alpha}\）)表示用于获得屈服函数上的分数应力梯度的偏导数（f）值得注意的是，卡普托分数导数具有奇异核，这可能不适用于更具体的情况，例如捕捉不同尺度的材料异质性和结构。这个极限的一个可能的解决方案是提出一个新的没有奇异核的分数阶导数，例如Caputo–Fabrizio导数[35]和Yang–Srivastava–Machado衍生物[28]. 值得指出的是，Yang–Srivastava–Machado分数导数[28]是Riemann–Liouville分数阶导数的一个著名推广，具有奇异核，在求解局部热学和力学问题方面具有很大潜力。然而，在本研究中，只考虑了与塑性相关的力学现象；因此，原始的卡普托分数导数[33,34]足够且有效[35]. Following Sun等人[27,36]，使用修改的Cam-clay（MCC）屈服函数

$$f=\bigl（2p'-p'_{0}\bigr）^{2}+\biggl（\frac{2q}{M}\biggr）^{2]-p_{0{^{\prime 2}=0$$

(10)

其中临界状态应力比M（M）可以表示为[37]

$$M=M_{c}\biggl[\frac{2c^{4}}{1+c^{4]-（1-c^{4}）\sin（3\theta）}\bigr]^{1/4}$$

(11)

在哪儿\（M _｛c｝\）表示三轴压缩的临界状态应力比；c（c）表示三轴拉伸的临界状态应力比(\（M _｛e｝\）)和\（M_{c}\）.Lode角(θ)表示为\（θ=\arccos\{9s_{ij}秒_{jk}秒{ki}/[2（3s_{lr}秒_{lr}/2）^{3/2}]\}/3\）.\（[-\pi/6，\pi/6]\中的θ），其中\（θ=pi/6）对应于三轴压缩，而\（θ=-\pi/6）对应于三轴拉伸。\（p'_{0}\）控制屈服面的大小（硬化）。此外，在方程式(9)，伽马函数定义如下：\（\varGamma（x）=\int_{0}^{\infty}e^{-\tau}\tau^{x-1}\，d\tau\）.α是分数阶，由于正塑性耗散的热力学限制，应至少小于2[27]. 正如孙和申所讨论的那样[21],α表示塑性流动与加载方向之间的非关联程度以及材料变形的状态依赖程度。\（西格玛）是当前有效应力，而\（西格玛{c}）是相应的临界状态应力。使用\（西格玛因为积分极限是由临界状态土力学的基本假设[10]：剪切时，土壤将持续变形，最终接近由临界状态线（CSL）表示的临界状态\（e-\ln p'）和\（p'-q\）平面分别如下：

$$\开始{aligned}&p'_{c}=p_{r}\exp\biggl（\frac{e_{varGamma}-e}{lambda}\biggr），结束{aligned}$$

(12)

$$\begin{aligned}&q_{c}=q+M\bigl（p'-p'{c}\bigr），\end{aligned}$$

(13)

哪里\（p_{r}\）是单位压力；λ是CSL在\（e-\ln p'）平面。\（e_{\varGamma}\）表示CSL的截距\（p'=1\mbox{kPa}）.\（p'\）和\（p'{c}\）分别是电流和相应临界状态下的平均有效主应力，而q个和\（q{c}\）分别是当前和相应临界状态下的偏应力。值得注意的是，在本研究中，我们连接了当前的应力点\（（p'，q）\）通过使用临界状态应力比与相应的CSLM（M）如图所示1在土壤达到临界状态之前，与CSL相关的当前应力点有三个可能的位置\（p'-q\）平面，即在“干”侧（上方）、“湿”侧（下方）和CSL上。如果当前状态位于CSL上，则\（p'=p'{c}\）和\（q=q{c}\），其中只指示了一个临界状态应力点，并且材料当前正在经历临界状态流动。然而，在其余两种情况下，可以观察到\（p'{c}\）和\（q{c}\）应代表CSL上的两个独立临界状态压力点（A和B）。如图所示1，应力点A表示有膨胀趋势的土壤，而应力点B表示有收缩趋势的土壤。\（p’_｛c｝\）通过从当前应力点（a或B）延伸水平平行线来获得\（p’\）-轴并与CSL相交，而\（q{c}\）通过将从当前应力点（a或B）延伸的垂直平行线与q个-轴。当前应力点与临界状态线的相对位置定义了哪个方程(8)或(9)将使用。更具体地说，如果\（p'>p'{c}\）或\（q>q{c}\），等式(8)应该使用，反之亦然。然而，如中所示[36]即使使用了当前应力点（或分数导数）的不同位置（或定义），最终也获得了唯一的状态相关塑性流动规则，而不使用任何塑性势和经验状态参数。

图1

当前应力点可能位置的示意图

三状态相关分数模型

3.1状态相关塑性流动

在本研究中，状态相关的分数塑性流动张量(\（n{ij}\）)定义为

$$n_{ij}=\frac{1}{\sqrt{1+d_{g}^{2}}\biggl[\frac{d_{g}}{3}\delta_{ij}+\frac{3s{ij{}}{2q}\bigr]$$

(14)

哪里\（d{g}\）是应力-延展比，不同研究人员对其有几种不同的定义，例如CC表达式[26,38]，椭圆表达式[12,39]，Rowe的表达式[6,40]然而，无论经典塑性模型中使用哪种剪胀方程，都必须加入经验状态指数，以统一建模岩土材料在广泛密度和压力下的状态相关行为[三]. 在本研究中，不使用经验状态指数和塑性势的三维状态相关剪胀比可以根据Sun等人[27,36]，其中统一的应力-直径比(\（d{g}\）)可以通过使用等式获得(8)和(9)在MCC表面进行分数阶导数，以便

$$开始{对齐}[b]d_{g}&=-\frac{{}}{p'}d{p'{c}}^{\alpha}f（p'）}{{}_{q{c}{d{q}^{alpha}f（q）}=-\frac{{{}{p'{c}d{p'}^{阿尔法}f（p’）}{}_{q} D类{q{c}}^{α}f（q）}\\&=M^{1+\alpha}\frac{（p'-p'{c}）+（2-\alpha）（p'{c{-p'{0}/2）}{（q-q{c{）+$$

(15)

其中，材料流动被发现受到几个因素的影响，包括洛德角(θ)通过M（M），电流应力\（（p'，q）\）和临界状态应力\（（p′{c}，q′{c{）\）; 最重要的是，它还取决于从当前状态到相应临界状态的应力距离(\（p'-p'_｛c｝\）和\（q-q{c}\）)。很容易找到\（d_{g}=0\）处于临界状态，其中\（p'=p'{c}\）和\（q=q{c}\），表示无塑性体积应变。MCC上的分数梯度，使用公式(8)和(9)，本质上考虑了剪切过程中土壤的状态信息。因此，与经典的各向同性和各向异性塑性模型不同[41,42,43,44,45]其中，材料状态对塑性流动的影响需要经验纳入状态参数（例如。ψ)在本研究中，对状态相关塑性流动进行了数学推导。通过使用公式(15)，影响α可以得到颗粒土在剪切过程中的应力-扩散特性。如Sun等人所示[36]，应力剪胀曲线表现出顺时针旋转和向下移动α增加。它简化为经典的MCC应力-微分模型[10]带有\（阿尔法=1）.

还应注意，应力-延展比只有两种可能性\（d_{g}=0\）一个处于临界状态，另一个处于相变状态。处于临界状态，其中\（p'=p'{c}\）和\（q=q{c}\），应力-直径比自动等于零，表示无塑性体积应变。在相变状态下，\（d{g}=0\）可以通过适当的值来确保α，可通过重新排列等式来定义(16)如下：

$$\阿尔法=\压裂{2M^{2} 第页'{ct}-2\eta{t}^{2} 第页'_{t}}{2M^{2} 第页'_{ct}-M^{2} 第页'{t}-\eta{t}^{2} 第页'_{t}}$$

(16)

在哪儿\（p'{ct}=p_{r}\exp（（e_{\varGamma}-e{t}）/\lambda）\）；\（e_{t}\）,\（p'{t}\）和\（\ta{t}\）分别为相变状态下的孔隙比、有效平均主应力和应力比。由于使用了垂直(\（q-q{c}\）)和水平(\（p'-p'_｛c｝\）)距离通过α相变线不是固定的，而是随着物质的初始状态而改变的。为了进一步澄清，对α带有状态参数(米)在Li和Dafalias[1]可以制作，其中米用于捕捉状态相关的塑性流动，并确保相变线随材料状态的变化。因此，对状态相关的塑性流动进行了分析。

3.2可能的数学关系\（I_{p}\）

关于与经验状态变量的数学联系的另一个有趣发现\（I_{p}\）[8]可以通过进一步代入公式(13)到等式(15)。然后，可以将上述应力-微分方程重排如下：

$$开始{对齐}[b]d_{g}&=M^{1+\alpha}\frac{（I{p}-1）+（2-\alpha）[1-I{p}/2-I{p{\eta^{2}/（2M^{2{）]}{（1-I{p{）M+（2-\ alphaα/2-（1-\alpha/2）\eta^{2}/M^{2{]I_{p}+（1-\alpha）}{（1-I_{p}）+（2-\alpha）[q/Mp'_{c}+（I{p}-1）]}，\结束{对齐}$$

(17)

哪里\（I_{p}=p'/p_{c}\）是Wang等人定义的状态压力指数[8]. 导出的应力-直径比与经验状态压力指数之间的数学关系\（I_{p}\）观察到。与其他研究不同[46]根据经验合并的\（I_{p}\），目前的研究从自然数学推导中反映了状态压力的依赖性。然而，应该指出\（d{g}\）和\（I_{p}\）特别是将此模型与其他模型联系在一起[1,13,14,15,41,47]使用的\（I_{p}\）或ψ[9]通过对特定屈服面（例如MCC）进行分数导数来推导这种关系的附加优点。当使用其他屈服面时，这种连接通常也存在吗？需要进行进一步的研究。然而，应该强调的是，由于分数导数的积分定义，本研究中提出的分数方法内在地考虑了从当前状态到临界状态的状态信息，如等式中所示分数导数的初始定义所示(8)和(9)。通过对其他屈服面进行分数导数，也可以从理论上导出状态相关的应力剪胀方程。

3.3边界面和加载方向

为了简单起见，边界曲面(f̄)假设与屈服面形状相同，即。

$$\bar{f}=\bigl（2\bar{p}'-\bar{p{0}\biger）^{2}+\biggl（\frac{2\bar}{q}}{M}\biggr）^{2]-\bar{p}_{0}^{\prime 2}=0$$

（18）

哪里\（\bar{p}'_{0}\）是边界曲面与横坐标的截距，控制边界曲面的大小。图像应力点\（（\bar{p}'，\bar{q}）\）边界曲面上可以用标量表示ρ作为

$$\开始{aligned}&\bar{p}'=\rho\bar{p2}'_{0}\end{aligned}$$

(19)

$$\begin{aligned}&\bar{q}=\rho\eta\bar{p}'_{0}，\end{aligned}$$

(20)

其中应力比η可以使用径向映射规则定义[39]如下：

$$\eta=\frac{q}{p'}=\frac}\bar{q}}{\bar{p}'}$$

(21)

边界内表面塑性[48]，加载张量(\（m{ij}\）)垂直于边界曲面，因此可以通过对边界曲面函数进行一阶导数获得，如下所示：

$$m_{ij}=\压裂{1}{\sqrt{1+d_{f}^{2}}\biggl[\压裂{d_{f}}{3}\delta_{ij}+\压裂{3s{ij{}}{2q}\bigr]$$

(22)

其中荷载比\（d{f}\）公式化为

$$d_{f}=\frac{M^{2}-\eta^{2{2}{2\eta}$$

(23)

进一步代入等式(19)–(21)到等式(18)，标量(ρ)其确定图像应力点，可以如下获得：

$$\rho=\frac{1}{1+（\eta/M）^{2}}$$

(24)

此外，初始边界曲面的位置(\（\bar{p}'_{0i}\）)可通过将正常压缩线和膨胀线相交于\（e-p’\）平面如下：

$$\bar{p}'_{0i}=2p_{r}\exp\biggl[\frac{e_{varGamma}-e_{0}-\kappa\ln p'_{ic}}{\lambda-\kappa}\biggr]$$

(25)

哪里\（e{0}\）是剪切前的初始孔隙比。\（p'{ic}\）是初始围压。边界曲面的演化(\（\bar｛p｝'_｛0｝\）)可进一步获得如下：

$$\bar{p}'_{0}=\bar{p2}'_0i}\exp\biggl（\frac{1+e_{0}}{\lambda-\kappa}\varepsilon_{v}^{p}\biggr）$$

(26)

请注意，公式的详细推导(25)和(26)可以在孙和沈找到[21]和Sun等人[49]，因此为了简单起见，这里不再重复。

3.4硬化模量

根据Dafalias的说法[48]，硬化模量H（H）由边界表面的尺寸（硬化）以及加载表面和边界表面之间的距离来确定。H（H）还观察到取决于物质状态，其中状态变量通常是通过经验合并的，例如Wang等人[8]. 它通常分解为两个组件[12,39,47]:

$$H=H_{b}+H_{delta}$$

（27）

哪里\（H_{b}\）通过在经历各向同性硬化的边界表面上应用一致性条件来确定：

$$H_{b}=-\frac{\partial\bar{f}}{\parial\bar{p'}_{0}}\frac}\partial/bar{p'}_{0}}{\ partial\ varepsilon_{v}^{p}}\frac{n_{v{}{\Vert\frac{\ parial\bar{f}{\protial\bar{\boldsymbol{\sigma}}}\Vert}=\frac[1+e_{0}{\lambda-\kappa}\frac{\bar{p'}M_{c}^{2} d日_{g} \\sqrt{d_{g}^{2}+1}}{\sqrt{（2\rho-1）^{2} M（M）_{c} ^{4}+4\rho^{2}\eta^{2{}}$$

(28)

很容易找到\（H_{b}\）由于依赖于\（d{g}\）物质状态。\（H_{\增量}\）与距离之间的比率有关(δ)从当前应力点到图像应力点的距离(\（δ_｛\max｝\）)从应力原点到图像应力点[21,39]:

$$H_{\delta}=H_{0}页'\frac{1+e_{0}}{\lambda-\kappa}\frac}\delta}{\delta_{\max}-\delta}$$

(29)

哪里\（h{0}=h_{1} e（电子）-h{2}\）；\（h{1}\）和\（h _｛2｝\）是材料常数。如预期，硬化模量\（H=+\infty\）在荷载开始时，其中\（\delta_{\max}-\delta\到0\）,\（eta=0）和\（\dot{\eta}=0\），表示无塑性应变状态。\（H=0）处于临界状态，其中\（p'=p'{c}\）,\（q=q{c}\）,\（\bar{p}'{0}=p'{0{）.

4参数识别和敏感性分析

共有九个参数(\（M_{c}\）,λ,\（e_{\varGamma}\）,c（c）,α,\（h{1}\）,\（h{2}\）,κ,ν)在所提出的状态相关模型中，所有这些都可以从传统的三轴试验中确定。下面详细阐述了参数识别和灵敏度分析。

与经典塑性模型类似，有四个临界状态参数，即。\（M _｛c｝\）,λ,\（e_{\varGamma}\）和c（c）.临界状态应力比(\（M_{c}\）)通过测量临界状态线的梯度来确定\（p'-q\）平面。随着\（M_{c}\），颗粒土的预测抗剪强度增加，而体积膨胀减小，如图所示2注意，图中用于模拟的其他模型参数2是：\（λ=0.057）,\（e_{varGamma}=0.64\）,\（c=1）,\（α=1.12）,\（h _｛1｝=h _｛2｝=2.5 \）,\（卡帕=0.0013）和\（nu=0.25\）.

图2

的影响\（M_{c}\）颗粒土的应力应变特性

λ和\（e_{\varGamma}\）是中的关键状态参数\（e-\ln p'）平面。λ可以通过测量临界状态线的梯度来获得\（e-\ln p'）平面，而\（e_{\varGamma}\）可以通过临界状态线的截距来确定\（p'=1）.c（c）可以通过进一步进行三轴拉伸试验来确定临界状态应力比(\（M_{e}\）)进行三轴拉伸，然后\（c=M{e}/M{c}\）如图所示三，随着λ预测的峰值应力和体积膨胀减小。然而，如图所示4，随着\（e_{\varGamma}\），预测的峰值应力和体积膨胀增加。请注意，图中用于模拟的其他模型参数三和4是：\（M{c}=1.7\）,\（c=1）,\（α=1.12）,\（h{1}=h{2}=2.5\）,\（卡帕=0.0013）和\（nu=0.25\）.

图3

的影响λ颗粒土的应力应变特性

图4

的影响\（e_{\varGamma}\）颗粒土的应力应变特性

有一个新参数，分数阶α这与传统塑性模型中的不同。它控制颗粒土的状态相关塑性流动。因此，它可以通过使用相变状态下的应力剪胀比来获得，即方程(16)。如图所示5，随着α预测的峰值应力随着体积膨胀的增大而减小。由更高级别的α更快地达到临界状态；当α从1.0增加到1.4。注意，图中用于模拟的其他模型参数5是：\（M{c}=1.7\）,\（c=1）,\（λ=0.057）,\（e_{varGamma}=0.64\）,\（h{1}=h{2}=2.5\）,\（卡帕=0.0013）和\（nu=0.25\）.

图5

的影响α颗粒土的应力应变特性

硬化参数\（h{0}\）与孔隙比相关\（h{1}\）和\（h{2}\）如Li和Dafalias所示[1]还有孙和申[21]，可以通过拟合\（瓦雷普西隆{1}-q\）不同孔隙比样品之间的关系。影响\（h{0}\）关于粒状土的应力应变特性，见图6其中，剪切强度增加，体积膨胀增加。注意，图中用于模拟的其他模型参数6是：\（M{c}=1.7\）,\（c=1\）,\（λ=0.057）,\（e_{varGamma}=0.64\）,\（α=1.12）,\（卡帕=0.0013）和\（nu=0.25\）.

图6

的影响\（h{1}\）和\（h{2}\）颗粒土的应力应变特性

弹性常数κ可以通过测量膨胀线的梯度来确定\（e-\ln p'）平面。如图所示7，随着κ峰值偏应力增加，而初始偏应力随着κ此外，还观察到\（瓦雷普西隆{v}）和\（瓦雷普西隆{1}）.在初始加载阶段，其中\（瓦雷普西隆{1}）是小的，对于具有高κ然而κ在\（瓦雷普西隆{v}）当\（瓦雷普西隆{1}）增加，其中塑料部件\（瓦雷普西隆{1}）和\（瓦雷普西隆{v}）成为主导。

图7

的影响κ颗粒土的应力应变特性

泊松比ν定义了材料的横向变形能力，可通过在初始加载阶段使用以下方程获得[12]:

$$\nu\approx-\frac｛\varepsilon\｛3｝｝｛\varepsilon\｛1｝｝$$

(30)

影响ν被发现是有限的，因此此处不作简单介绍。然而，根据弹性理论，可以观察到更大的侧向变形ν表中列出了用于预测不同颗粒土应力应变行为的每个模型参数的详细值1值得注意的是，对于模拟三轴试验\（θ=pi/6）,c（c）将不会使用。

表1模型参数

5模型验证

本节展示了模型捕捉不同颗粒土（包括不同初始状态的堆石和砂）状态相关本构行为的能力。塔城堆石料的一系列三轴试验结果[12]萨克拉门托河砂[51]和渥太华沙[50]在图中进行了模拟8——14.塔城堆石料排水应力应变特性[12]图中模拟了四种不同的初始孔隙比8——11，而图12——14提出萨克拉门托河砂排水和不排水三轴特性的模型预测[51]和渥太华沙[50]. 值得注意的是，所有测试结果均由离散数据点表示，而连续线用于模型预测。

图8

塔城堆石体排水特性的模型预测\（e{0}=0.189）（数据来自Xiao等人[12])

图9

塔城堆石体排水特性的模型预测\（e{0}=0.244）（数据来自Xiao等人[12])

图10

塔城堆石排水特性的模型预测\（e{0}=0.285）（数据来自Xiao等人[12])

图11

塔城堆石体排水特性的模型预测\（e{0}=0.317）（数据来自Xiao等人[12])

图12

萨克拉门托河砂排水特性的模型预测[51]

Xiao等人[12]报道了塔城堆石料不同初始孔隙比的排水三轴试验结果。该材料主要由中等直径的亚棱角至圆形颗粒组成(\（d{50}）)23 mm和均匀系数(\（C_{u}\）)第5.4节。样品通过分层压实制备，直径约为300 mm，高度约为600 mm。初始孔隙比和相应的围压见图8——11因此这里不再重复。从图中可以观察到8——11即使不使用状态变量和塑性势函数，该模型也能很好地模拟塔城堆石料在不同初始状态（围压和孔隙比）下的状态相关本构行为。塔城堆石料的应变硬化和软化特性以及相应的体积收缩和膨胀特性都可以得到合理的描述。

Lee和Seed[51]对萨克拉门托河砂在不同初始条件下的排水和不排水三轴特性进行了综合研究。萨克拉门托河砂主要由亚圆形石英骨料和偶尔的贝壳碎片组成\（d{50}）0.22 mm和a\（C_{u}\）等于1.45。经测试，最小和最大孔隙比分别为0.61和1.03。样品的初始直径约为42.67 mm，高度约为103.63 mm。初始孔隙比0.87用于模拟排水试验结果，而初始孔隙比为0.86、0.86、086和0.85用于模拟分别在98、294、490和1069 kPa围压下进行的不排水试验。图12显示了萨克拉门托河砂排水本构行为的模型模拟，可以观察到模型模拟与相应试验结果之间的良好一致性。图13显示了萨克拉门托河砂在不同初始条件下不排水行为的模型模拟。随着轴向应变的增加，模拟的偏应力增加，直到达到临界（稳定）状态流，这与相应的试验结果非常一致，这进一步验证了新开发的状态相关应力-直径关系。

图13

萨克拉门托河砂不排水特性的模型预测[51]

图14给出了均匀渥太华砂不排水本构行为的模拟结果[50]主要由圆形/亚圆形石英骨料组成。模拟了初始孔隙比分别为0.793、0.793和0.805的渥太华砂在348、475和550kPa相应围压下的试验结果。结果表明，模型预测与试验结果吻合良好。模拟的偏应力最初增加，直到达到峰值，然后在临界状态下显著降低，表明材料处于静态液化状态。具有相同初始孔隙比的样品在进一步剪切时接近相同的偏应力，所有这些都可以通过所提出的分数塑性模型合理地描述。

图14

渥太华砂不排水特性的模型预测[50]

6结论

研究发现，颗粒土的应力应变特性与状态有关。为了捕捉这种状态依赖性，通过将经验状态变量纳入附加塑性势函数，提出了许多复杂的本构模型。为了简化建模方法，通过对屈服函数进行分数阶导数，提出了一种新的状态相关分数塑性模型，该模型不使用预定义的状态变量和塑性势。然后对模型参数进行了详细的识别和敏感性分析。为了进一步验证所建立的模型，对不同颗粒土的一系列排水和不排水试验结果进行了模拟和讨论。主要调查结果总结如下：

• (1) 在不使用预先定义的状态变量和塑性势函数的情况下，利用分数阶塑性理论解析地导出了一个新的状态相关应力-微分方程。通过分数应力梯度的严格数学定义，模拟了非关联流动对材料状态的依赖性。

• (2) Wang等人提出的状态相关扩容方程和状态压力指数之间可能存在的数学联系[8]文中还讨论了状态压力对颗粒土应力-扩散现象的依赖性。

• (3) 材料状态通过从当前应力状态到相应临界应力状态的垂直和水平距离影响非结合程度和硬化模量\（p'-q\）平面。

• (4) 随着分数阶数的增加，预测的峰值应力减小，而体积膨胀增加。分数阶越高的样本越快达到临界状态；随着分数阶的增加，从应变软化行为到应变硬化行为的转变也增加。

• （5） 模型参数均可根据传统三轴试验结果确定。结果表明，所提出的状态相关分数塑性模型能够很好地描述不同颗粒土在各种荷载条件下的应力应变行为。