考虑剪胀性的面层基层混合铣刨料模型及其参数敏感性分析

Modeling of Mixed Milling Materials for Surface Base Layers Considering Shear Expansion and Sensitivity Analysis of Its Parameters

摘要

为了研究面基混合铣刨料的力学性能及其尚不明确的本构模型，在三轴试验的基础上，基于Duncan-Chang模型和Pastor-Zienkiewicz剪胀方程构建了非线性剪胀本构模型，通过应用次弹性理论推导出该本构模型的刚度矩阵，并借助ABAQUS有限元软件平台的UMAT子程序对该模型进行了二次开发，验证了模型的准确性和可靠性.此外，还对该模型的参数进行了全体敏感性分析及关键参数响应面分析.结果表明，在轴向应变为10%左右时，不同围压下两种材料的强度均逐渐达到残余强度，体积应变与轴向应变呈现出先减小后增大的趋势，且在一定范围内表现为明显的剪胀特性.利用二次开发出的UMAT子程序对三轴试验进行模拟计算，其结果与试验结果相比展现出高度的吻合度.研究指标偏差应力（σ₁-σ₃）受到参数φ₀、Δφ、K的影响尤为突出，受参数n、φ_cr、α、R_f的影响则相对较弱，轴向应变ε₁受参数R_f、K、φ₀的影响变得显著，而n、Δφ、φ_cr、α的影响则显得较为有限.对（σ₁-σ₃）的敏感性程度由高到低次序为φ₀、φ₀²、Δφ、φ₀Δφ、KΔφ、K、K²、Δφ²、Kφ₀，对ε₁的敏感性程度由高到低次序为R_f、φ₀、R_f²、K、R_fφ₀、KR_f、φ₀²、Kφ₀、K².通过响应面法分析比较关键参数之间交互作用，发现它们对（σ₁-σ₃）的具体影响程度为：φ₀>Δφ>K；对ε₁的具体影响程度为：R_f>φ₀>K.研究成果可为面基混合铣刨料的工程应用提供理论依据及为参数选取提供参考.

Abstract

In order to investigate the mechanical properties of face-based mixed milling planes and their ontological model that is not unspecified yet. Based on the triaxial test， a nonlinear shear expansion constitutive model was constructed based on the Duncan-Chang model and the Pastor-Zienkiewicz shear expansion equation， and the stiffness matrix of this constitutive model was derived by applying the theory of infra-elasticity. Moreover， the model was developed twice with the help of the UMAT subroutine of the ABAQUS finite element software platform， which verified the model accuracy and reliability of the model. In addition， the overall sensitivity analysis and the response surface analysis of the key parameters of the model were also carried out. The results show that at about 10% axial strain， the strength of the two materials gradually reaches the residual strength under different circumferential pressures， and the volumetric strain and axial strain show a tendency to decrease and then increase， and exhibit obvious shear expansion characteristics within a certain range. The triaxial tests were simulated using a second developed UMAT subroutine， and the results showed a high degree of agreement when compared to the test results. The study index deviatoric stresses （σ₁-σ₃） are particularly affected by the parameters φ₀， \mathrm{\Delta }φ， and K， while they are relatively weakly affected by the parameters n， φ_cr， α， and R_f， and the axial strain ε₁ becomes significantly affected by the parameters R_f， K， and φ₀， while the effects of n， Δφ， φ_cr， and α appear to be more limited. The degree of sensitivity to （σ₁-σ₃） is φ₀， {\phi }_{\mathrm{0}}^{\mathrm{2}}， \mathrm{\Delta }\phi， φ₀\mathrm{\Delta }\phi， K\mathrm{\Delta }\phi， K， K²， \mathrm{\Delta }{\phi }^{\mathrm{2}}， K{\phi }_{\mathrm{0}} in descending order， and the degree of sensitivity to ε₁ is R_f， φ₀， R_{\mathrm{f}}^{\mathrm{2}}， K， R_fφ₀， KR_f， {\phi }_{\mathrm{0}}^{\mathrm{2}}， Kφ₀， K² in descending order. The interactions of the key parameters were analyzed through the response surface methodology for comparing the interaction effects between them， and it is found that they are not significant for （σ₁-σ₃） is： φ₀ >\mathrm{\Delta }\phi >K； the specific influence on ε₁ is： R_f >φ₀ >K. The research results can provide a theoretical basis for the engineering application of surface mixed milling planer and a reference for parameter selection.

Graphical abstract

关键词

本构模型 / 剪胀性 / 次弹性理论 / 敏感性分析 / 响应面

随着道路交通流量的激增与车辆荷载的加重，造成道路面层和基层损坏问题日益突出，给道路的安全、畅通带来了严峻挑战^［1-4］.在道路的后续翻新修补中，机械铣刨作业会产生面层沥青与基层碎石混合的铣刨料，这种面基混合铣刨料因含有基层碎石料，在荷载作用下具有非线性、压硬性、摩擦性及应力路径相关性等特性^［5-7］，此外，碎石料在低密度和高围压条件下具有剪缩性，而在高密度和较低围压条件下还具有剪胀性^［8-9］.因此，在研究面基混合铣刨料的力学响应及构建其本构模型时，必须既考虑其非线性，也反映剪缩与剪胀的复杂特性.

目前在岩土工程领域中应用广泛的本构模型有：Duncan等^［10-11］提出的E-ν型和用体积模量代替泊松比的E-B型增量弹性的Duncan-Chang模型；K-G模型^［12］通过将应力和应变分解为球张量和偏张量，建立了平均应力dp与体积应变\mathrm{d}{\epsilon }_{\mathrm{v}}、广义剪应力dq与广义剪应变\mathrm{d}\overline{\epsilon }之间的关联；在K-G模型的基础上，研究者们改进和扩展出了Domaschuk-Valliappan模型^［13］、Naylor模型^［14］、Izumi-Verruijt耦合模型（三参数模型）^［15］、四川大学K-G模型^［16］和清华非线性解耦K-G模型^［17］等；Roscoe等^［18］基于深入的试验观察、临界状态概念的引入以及加工硬化原理和能量方程等提出剑桥模型；Burland^［19］对剑桥模型引入椭圆屈服面、相适应的流动准则，以塑性体应变为硬化参数以及基于试验数据的分析和比较提出修正剑桥模型；沈珠江^［20］在吸收了邓肯模型和剑桥模型的优点后提出了沈珠江双屈服面弹塑性模型（南水模型）；殷宗泽^［21］在深入研究土体复杂应力应变关系的基础上，对剪胀性与剪缩性、硬化与软化特性以及工程实践的需求等多方面因素的综合考虑提出殷宗泽双屈服面弹塑性模型；姚仰平等^［22-23］依据土的临界状态理论、现代塑性力学理论和试验数据分析提出了统一硬化（UH）本构模型；程展林等^［24-25］提出了一种新的“三参量”非线性剪胀模型，该模型可视为邓肯-张模型的改进型模型等.现阶段对于土体本构模型的研究已取得了丰硕的成果，但每种模型均存在一定的适用范围和限制条件，而针对适用于面基混合铣刨料的本构模型还不明确，找到合适的本构模型对基混合铣刨料的工程再应用具有重要的意义.

本文在碎石料和面基混合铣刨料大型三轴试验的基础上，基于Duncan-Chang模型和Pastor- Zienkiewicz剪胀方程构建出非线性剪胀本构模型，通过应用次弹性理论推导出该本构模型的刚度矩阵，借助ABAQUS有限元软件平台的UMAT子程序对该非线性剪胀模型进行了二次开发，并将模拟结果与试验结果进行对比，以验证该模型的准确性和可靠性.此外，还对非线性剪胀模型的参数进行了全体敏感性分析及关键参数响应面分析.研究成果可为面基混合铣刨料的工程应用提供理论依据及为参数选取提供参考.

1 材料力学性能试验

1.1 试验材料

试验材料选用碎石料（作参照）和面基混合铣刨料（面层沥青铣刨料和基层碎石料按照质量1∶1比例混合），通过筛分比试验测定，碎石料粒径主要集中在5~40 mm，面基混合铣刨料粒径主要集中在0.075~40 mm，其颗粒相对密度分别为2.75、2.62，限制粒径d₆₀分别为15.5 mm和10.0 mm，不均匀系数分别为1.61、19.23，曲率系数分别为1.05、1.73.根据《土工试验方法标准》（GB/T 50123—2019）的规定，测得碎石料和面基混合铣刨料的最大干重度分别为17.39 kN/m³和18.31 kN/m³，最小干重度分别为 13.3 kN/m³和13.2 kN/m³，含水率分别为12%和10%.试验材料样品及其颗粒级配曲线如图1所示.

1.2 试验仪器及方案

试验使用的仪器为大型三轴剪切试验系统，主要有三轴压力室、轴压系统、侧压系统和孔隙水压力测读系统等.其主要技术指标为：围压（{\sigma }_{\mathrm{c}}）0~ 3.0 MPa，轴向应力（{\sigma }_{\mathrm{a}}）0~21 MPa，轴向行程（U_{\mathrm{a}}）0~300 mm.试验采用固结排水剪切试验（CD），试样为圆柱形，其尺寸采用直径300 mm×高度600 mm.

针对碎石料、面基混合铣刨料，分别开展围压为100 kPa、200 kPa、300 kPa、400 kPa的大型三轴剪切试验.试样制备、试验前准备和加载试验工作严格按照《土工试验方法标准》（GB/T 50123—2019）的规定执行.将级配均匀的试验材料装填入预先准备好的筒内，分5层均匀摊铺并振动夯实，确保压实度均达到96%；扎紧橡皮膜，安装加压帽并将其抽至真空；施加围压开始试验，剪切过程中以0.4 mm/min的速率进行加载，直到试样剪切至轴向应变的25%时或试样完全破坏，停止试验.其部分试验过程如 图2所示.

1.3 试验结果

图3为碎石料和面基混合铣刨料的偏差应力-轴向应变曲线及体积应变-轴向应变曲线.由图可知，不同围压下两材料偏差应力-轴向应变关系曲线的变化趋势较为相似，均具有明显偏差应力峰值特征，随着轴向加载不断增加，表现出弹性变形阶段、塑性变形阶段和试样破坏阶段；而体积应变-轴向应变关系曲线在加载过程中呈现出体积应变随轴向应变增加而先减小后增大的特征，即在一定范围内先出现剪胀效应，随后受到外部作用力影响逐渐恢复并增大体积应变，整体呈现出非线性的曲线特征.

从图3（a）、（c）可知，从试验开始加载到试样轴向应变{\epsilon }_{\mathrm{1}}达到0.02左右时，两材料的偏差应力与轴向应变成正比，其关系曲线基本呈直线，在此阶段可近似为线弹性材料；当两材料的变形超过弹性轴向应变临界点后，偏差应力的增长速率要小于轴向应变增加速率，表现为{\epsilon }_{\mathrm{1}}~\left({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}\right)关系曲线开始不断向轴向应变轴偏移，两材料强度接近极限强度时，偏差应力的增长速率迅速变小，相反轴向应变增加速率急剧增大，塑性特性更显著；当两材料抗剪强度达到极限时，试样将不能承受更大的加载，此时试样发生剪切破坏，达到的强度值即为峰值抗压强度，与之相对应的应变值则为峰值应变.此外，当轴向应变达到0.01左右时，两材料强度趋于一个定值，即残余强度，且碎石料的残余强度略大于面基混合铣刨料（约5%），低围压的残余强度小于高围压的残余强度（25%~50%）.

从图3（b）、3（d）可观察到，随着围压从100 kPa逐渐增大至400 kPa，碎石料及面基混合铣刨料的体积应变随轴向应变的增大而呈现不同程度的增长，且两材料均展现出了典型的先剪缩后剪胀行为.围压的增加导致剪切初期体积应变值显著上升，同时剪缩段的长度也相应增加.值得注意的是，当围压增大时，碎石料及面基混合铣刨料在后期的剪胀段呈现减小趋势，其体积应变与轴向应变的关系曲线逐渐向剪缩方向转换，特别在围压达到400 kPa时，这两种材料在剪切后期表现出显著的剪缩特性，这与传统粗粒土剪切试验所揭示的规律相吻合.同时，表明面基混合铣刨料与碎石料的力学特性基本一致.

2 非线性剪胀本构模型

2.1 Duncan-Chang模型

为了本文理论模型的完整性，仍对Duncan- Chang模型作简单介绍^［10］.Duncan-Chang模型是一种建立在增量广义胡克定律基础上的非线性弹性模型，由邓肯等人根据双曲线应力应变关系提出.它能够反映土体变形的非线性特征，模型参数物理意义明确，且可通过静三轴试验确定.其切线变形模量E_{\mathrm{t}}的表达式为：

E_{\mathrm{t}}=K{p}_{\mathrm{a}}{\left(\frac{{\sigma }_{\mathrm{3}}}{p_{\mathrm{a}}}\right)}^n{\left[\mathrm{1}-\frac{R_{\mathrm{f}}\left({\sigma }_{\mathrm{1}}-{\sigma }_{\mathrm{3}}\right)(\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\phi )}{\mathrm{2}c\mathrm{c}\mathrm{o}\mathrm{s}\phi +\mathrm{2}{\sigma }_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\phi }\right]}^{\mathrm{2}}（1）

式中：K为弹性模量基数；n为弹性模量指数；R_{\mathrm{f}}为破坏比；p_{\mathrm{a}}为大气压，p_{\mathrm{a}}=101.4 kPa；c和\phi分别为土体的黏聚力与内摩擦角，对于粗粒土，内摩擦角对限制压力的依赖关系可用式（2）表示.

\phi ={\phi }_{\mathrm{0}}-∆\phi \mathrm{l}\mathrm{g}\left(\frac{{\sigma }_{\mathrm{3}}}{p_{\mathrm{a}}}\right)（2）

式中：{\phi }_{\mathrm{0}}为参考围限压力下的内摩擦角；∆\phi为内摩擦角随围压的减小量.

2.2 Pastor-Zienkiewicz剪胀方程

在岩土工程领域，尽管Duncan-Chang模型在描述土体强度变形特性方面具有诸多优点，但其无法充分反映土体的剪胀性质.通过对碎石料和面基混合铣刨料进行的力学性能试验研究表明，这2种材料在低围压条件下均表现出明显的剪胀性.为了更准确地描述和预测两材料的剪胀行为，有必要引入专门的剪胀方程，从而更精确地预测其在受剪过程中的力学响应.

关于土体剪胀模型及理论研究成果较多.其中几种在工程实践中尤为常用，Nova等^［26］基于土体应力状态和颗粒结构，引入剪胀系数等参数，建立了剪切与体积应变间的定量关系.Rowe^［27］研究了剪胀比与应力比、初始孔隙比的关系，提出剪胀方程.徐舜华等^［28］在Rowe基础上进行了改进和方程拓展.Manzari等^［29］和Gajo等^［30］基于临界状态下土体特性提出了剪胀方程.这些方程在描述特定条件下的土体剪胀行为时具有优势，但也存在局限性.Nova剪胀理论未充分考虑岩土材料达到临界状态时的特性，在模拟接近或达到临界状态时会出现偏差.Rowe剪胀方程及徐舜华修正方程在高围压下可能低估材料压缩性，而在剪胀阶段可能高估其剪胀性，导致计算结果精确度偏低.Manzari和Gajo的剪胀方程包含较多参数，确定这些参数需依赖大量试验数据和复杂拟合过程，不利于实际工程应用.

为了方便模型参数计算能精确反映模型剪胀性，本文选用Pastor等学者^［31-32］提出的Pastor- Zienkiewicz剪胀方程，该方程基于弹塑性理论并引入相关状态参数，使其具有广泛的适用性和更高的预测精度，且该方程中仅有2个参数，均可通过常规三轴试验得到.Ravanbakhsh等^［33］、Cola等^［34］和Song等^［35］应用该方程精准预测了砂土的三轴试验结果.Heidarzadeh等^［36］和Mira等^［37］基于该方程的理论计算结果与试验数据的对比验证，认为该剪胀方程能够准确反映体积应变和广义剪应变增量之间的动态关系.该剪胀方程的一般表达式为：

\frac{\mathrm{d}{\epsilon }_{\mathrm{v}}}{\mathrm{d}{\epsilon }_{\mathrm{s}}}=\left(\mathrm{1}+\alpha \right)\left(M-\frac{q}{p}\right)（3）

式中：\mathrm{d}{\epsilon }_{\mathrm{v}}和\mathrm{d}{\epsilon }_{\mathrm{s}}分别是剪胀引起增量体积应变和广义剪应变；\alpha为试验确定的模型参数；M为p-q平面上临界状态线的斜率，在三轴压缩应力状态下可通过土体的残余内摩擦角估算：

M=\frac{\mathrm{6}\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{c}\mathrm{r}}}{\mathrm{3}-\mathrm{s}\mathrm{i}\mathrm{n}{\phi }_{\mathrm{c}\mathrm{r}}}（4）

式中：{\phi }_{\mathrm{c}\mathrm{r}}为土体的残余内摩擦角.

在三轴应力状态下，平均应力p和广义剪应力q的表达式分别为：

p=\frac{{\sigma }_{\mathrm{11}}+{\sigma }_{\mathrm{22}}+{\sigma }_{\mathrm{33}}}{\mathrm{3}}（5a）

q=\frac{\mathrm{1}}{\sqrt[]{\mathrm{2}}}\left[{\left({\sigma }_{\mathrm{11}}-{\sigma }_{\mathrm{22}}\right)}^{\mathrm{2}}+{\left({\sigma }_{\mathrm{22}}-{\sigma }_{\mathrm{33}}\right)}^{\mathrm{2}}+\right.

{\left.{\left({\sigma }_{\mathrm{11}}-{\sigma }_{\mathrm{33}}\right)}^{\mathrm{2}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}（5b）

2.3 非线性剪胀本构模型

程展林等^［24］学者为揭示土体在剪切过程中的剪胀行为特性，将土的压硬性描述为p对偏应变{\epsilon }_{\mathrm{s}}的影响，同时，将剪胀性阐释为q对体应变{\epsilon }_{\mathrm{v}}的影响.鉴于这种交叉影响，传统的广义胡克定律无法准确描述这种复杂的应力-应变关系，因此，提出了粗粒土非线性剪胀模型，其应力-应变关系增量形式为：

\mathrm{d}{\epsilon }_{\mathrm{v}}=\frac{\mathrm{d}p}{K_{\mathrm{p}}}+\frac{\mathrm{d}q}{K_{\mathrm{q}}}（6a）

\mathrm{d}{\epsilon }_{\mathrm{s}}=\frac{\mathrm{d}q}{G}（6b）

式中：\mathrm{d}{\epsilon }_{\mathrm{v}}和\mathrm{d}{\epsilon }_{\mathrm{s}}分别为体积应变和广义剪应变的增量；\mathrm{d}p、\mathrm{d}q分别为平均应力增量和广义剪应力增量；K_{\mathrm{p}}、K_{\mathrm{q}}、G分别称为体变模量、剪胀模量、剪切模量，以下对参数K_{\mathrm{p}}、K_{\mathrm{q}}、G进行具体推导.

2.3.1 体变模量K_{\mathrm{p}}

基于胡克定律的基本原理，可以得出体变模量K_{\mathrm{p}}的具体表达式为：

K_{\mathrm{p}}=\frac{E_{\mathrm{u}\mathrm{r}}}{\mathrm{3}(\mathrm{1}-\mathrm{2}\mu )}（7）

式中：\mu为材料弹性阶段的泊松比，假定其为常数^［38］；E_{\mathrm{u}\mathrm{r}}为土体的切线弹性模量，可由式（8）确定.

E_{\mathrm{u}\mathrm{r}}=K_{\mathrm{u}\mathrm{r}}{p}_{\mathrm{a}}{\left(\frac{{\sigma }_{\mathrm{3}}}{p_{\mathrm{a}}}\right)}^n（8）

式中：K_{\mathrm{u}\mathrm{r}}为试验确定的模型参数.

2.3.2 剪胀模量K_{\mathrm{q}}

通过Pastor-Zienkiewicz剪胀方程并结合三轴应变关系，由剪胀作用引起的体积应变增量\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{q}}可表示为：

\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{q}}=\frac{\mathrm{2}\left(\mathrm{1}+\alpha \right)\left(M-\frac{q}{p}\right)}{\frac{\mathrm{2}}{\mathrm{3}}\left(\mathrm{1}+\alpha \right)\left(M-\frac{q}{p}\right)+\mathrm{2}}\mathrm{d}{\epsilon }_{\mathrm{1}}^{\mathrm{q}}（9）

在三轴试验应力条件下，有：

\mathrm{d}{\epsilon }_{\mathrm{1}}^{\mathrm{e}}=\frac{\mathrm{1}}{E_{\mathrm{u}\mathrm{r}}}\mathrm{d}{\sigma }_{\mathrm{1}}（10）

\mathrm{d}{\epsilon }_{\mathrm{1}}=\frac{\mathrm{1}}{E_{\mathrm{t}}}\mathrm{d}{\sigma }_{\mathrm{1}}（11）

联立式（10）和式（11），可得剪胀作用引起的轴向应变\mathrm{d}{\epsilon }_{\mathrm{1}}^{\mathrm{q}}为：

\mathrm{d}{\epsilon }_{\mathrm{1}}^{\mathrm{q}}=\mathrm{d}{\epsilon }_{\mathrm{1}}-\mathrm{d}{\epsilon }_{\mathrm{1}}^{\mathrm{e}}=\frac{E_{\mathrm{u}\mathrm{r}}-E_{\mathrm{t}}}{E_{\mathrm{u}\mathrm{r}}{E}_{\mathrm{t}}}\mathrm{d}{\sigma }_{\mathrm{1}}（12）

\mathrm{d}q=\mathrm{d}{\sigma }_{\mathrm{1}}（13）

联立式（9）、式（12）和式（13），可得剪胀作用引起产生的剪胀模量K_{\mathrm{q}}：

K_{\mathrm{q}}=\frac{\mathrm{d}q}{\mathrm{d}{\epsilon }_{\mathrm{v}}^{\mathrm{q}}}=\frac{\frac{\mathrm{2}}{\mathrm{3}}\left(\mathrm{1}+\alpha \right)\left(M-\frac{q}{p}\right)+\mathrm{2}}{\mathrm{2}\left(\mathrm{1}+\alpha \right)\left(M-\frac{q}{p}\right)}\frac{E_{\mathrm{u}\mathrm{r}}{E}_{\mathrm{t}}}{E_{\mathrm{u}\mathrm{r}}-E_{\mathrm{t}}}（14）

2.3.3 剪切模量G

将式（9）代入式（3），可以得到剪胀作用引起的广义剪应变增量\mathrm{d}{\epsilon }_{\mathrm{s}}^{\mathrm{q}}的表达式为：

\mathrm{d}{\epsilon }_{\mathrm{s}}^{\mathrm{q}}=\frac{\mathrm{2}}{\frac{\mathrm{2}}{\mathrm{3}}\left(\mathrm{1}+\alpha \right)\left(M-\frac{q}{p}\right)+\mathrm{2}}\mathrm{d}{\epsilon }_{\mathrm{1}}^{\mathrm{q}}（15）

在三轴应力条件下，可以推导出由弹性作用引起的广义剪应变增量\mathrm{d}{\epsilon }_{\mathrm{s}}^{\mathrm{e}}的表达式为：

\mathrm{d}{\epsilon }_{\mathrm{s}}^{\mathrm{e}}=\frac{\mathrm{2}}{\mathrm{3}}\left(\mathrm{d}{\epsilon }_{\mathrm{1}}^{\mathrm{e}}-\mathrm{d}{\epsilon }_{\mathrm{3}}^{\mathrm{e}}\right)=\frac{\mathrm{2}}{\mathrm{3}}\left(\mathrm{1}+\mu \right)\mathrm{d}{\epsilon }_{\mathrm{1}}^{\mathrm{e}}（16）

广义剪应变增量\mathrm{d}{\epsilon }_{\mathrm{s}}的表达式为：

\mathrm{d}{\epsilon }_{\mathrm{s}}=\mathrm{d}{\epsilon }_{\mathrm{s}}^{\mathrm{e}}+\mathrm{d}{\epsilon }_{\mathrm{s}}^{\mathrm{q}}（17）

将式（15）~式（17）代入式（6b），可以推导出剪切模量G的具体表达式为：

G=\frac{\mathrm{d}q}{\mathrm{d}{\epsilon }_{\mathrm{s}}}=\frac{E_{\mathrm{u}\mathrm{r}}{E}_{\mathrm{t}}}{\frac{\left(E_{\mathrm{u}\mathrm{r}}-E_{\mathrm{t}}\right)}{\frac{\mathrm{1}}{\mathrm{3}}\left(\mathrm{1}+\alpha \right)\left(M-\frac{q}{p}\right)+\mathrm{1}}+\frac{\mathrm{2}}{\mathrm{3}}\left(\mathrm{1}+\mu \right)E_{\mathrm{t}}}（18）

3 模型三轴应力空间推广

3.1 刚度矩阵推导

为了推广该模型在三轴应力空间下的更广泛应用，充分利用非线性剪胀模型具有应力增量张量和应变增量张量之间一一对应的弹性关系特性，引入了次弹性理论^［12］（也称作最小弹性模型）.假设研究材料为各向同性，且变形是瞬时的，不随时间累积.则次弹性理论的三维关系张量形式为：

\mathrm{d}{\mathit{\sigma }}_{ij}={\mathit{D}}_{ijkl}\left({\sigma }_{mn}\right)\mathrm{d}{\mathit{\epsilon }}_{kl}（19）

式中：\mathrm{d}{\mathit{\sigma }}_{ij}\mathrm{、}\mathrm{d}{\mathit{\epsilon }}_{kl}分别为增量应力张量和增量应变张量；{\mathit{D}}_{ijkl}\left({\sigma }_{mn}\right)为次弹性理论的刚度矩阵.次弹性理论转换成矩阵形式可写为：

\left\{\begin{array}{l}\mathrm{d}{\sigma }_{\mathrm{11}}\\ \mathrm{d}{\sigma }_{\mathrm{22}}\\ \mathrm{d}{\sigma }_{\mathrm{33}}\\ \mathrm{d}{\sigma }_{\mathrm{12}}\\ \mathrm{d}{\sigma }_{\mathrm{23}}\\ \mathrm{d}{\sigma }_{\mathrm{31}}\end{array}\right\}=

\left[\begin{array}{llllll}{D}_{\mathrm{1111}}& D_{\mathrm{1122}}& D_{\mathrm{1133}}& D_{\mathrm{1112}}& D_{\mathrm{1123}}& D_{\mathrm{1131}}\\ D_{\mathrm{2211}}& D_{\mathrm{2222}}& D_{\mathrm{2233}}& D_{\mathrm{2212}}& D_{\mathrm{2223}}& D_{\mathrm{2231}}\\ D_{\mathrm{3311}}& D_{\mathrm{3322}}& D_{\mathrm{3333}}& D_{\mathrm{3312}}& D_{\mathrm{3323}}& D_{\mathrm{3331}}\\ D_{\mathrm{1211}}& D_{\mathrm{1222}}& D_{\mathrm{1233}}& D_{\mathrm{1212}}& D_{\mathrm{1223}}& D_{\mathrm{1231}}\\ D_{\mathrm{2311}}& D_{\mathrm{2322}}& D_{\mathrm{2333}}& D_{\mathrm{2312}}& D_{\mathrm{2323}}& D_{\mathrm{2331}}\\ D_{\mathrm{3111}}& D_{\mathrm{3122}}& D_{\mathrm{3133}}& D_{\mathrm{3112}}& D_{\mathrm{3123}}& D_{\mathrm{3131}}\end{array}\right]\left\{\begin{array}{l}\mathrm{d}{\epsilon }_{\mathrm{11}}\\ \mathrm{d}{\epsilon }_{\mathrm{22}}\\ \mathrm{d}{\epsilon }_{\mathrm{33}}\\ \mathrm{d}{\epsilon }_{\mathrm{12}}\\ \mathrm{d}{\epsilon }_{\mathrm{23}}\\ \mathrm{d}{\epsilon }_{\mathrm{31}}\end{array}\right\}（20）

假设所研究的材料具备各向同性特性，则{\mathit{D}}_{ijkl}\left({\sigma }_{mn}\right)的一般形式为：

\begin{array}{l}{\mathit{D}}_{ijkl}\left({\sigma }_{mn}\right)={\beta }_{\mathrm{1}}{\delta }_{ij}{\delta }_{kl}+{\beta }_{\mathrm{2}}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{jk}{\delta }_{il}\right)+{\beta }_{\mathrm{3}}{\sigma }_{ij}{\delta }_{kl}+{\beta }_{\mathrm{4}}{\delta }_{ij}{\sigma }_{kl}+{\beta }_{\mathrm{5}}\left({\delta }_{ik}{\sigma }_{jl}+{\delta }_{il}{\sigma }_{jk}+{\delta }_{jk}{\sigma }_{il}+\left.{\delta }_{jl}{\sigma }_{ik}\right)\right.+{\beta }_{\mathrm{6}}{\delta }_{ij}{\sigma }_{km}{\sigma }_{ml}+{\beta }_{\mathrm{7}}{\delta }_{kl}{\sigma }_{im}{\sigma }_{mj}+{\beta }_{\mathrm{8}}\left({\delta }_{ik}{\sigma }_{jm}{\sigma }_{ml}+{\delta }_{li}{\sigma }_{jm}{\sigma }_{mk}+{\delta }_{jk}{\sigma }_{im}{\sigma }_{ml}+\right.\left.{\delta }_{jl}{\sigma }_{im}{\sigma }_{mk}\right)+\\ {\beta }_{\mathrm{9}}{\sigma }_{ij}{\sigma }_{kl}+{\beta }_{\mathrm{10}}{\sigma }_{ij}{\sigma }_{km}{\sigma }_{mi}+{\beta }_{\mathrm{11}}{\sigma }_{im}{\sigma }_{mj}{\sigma }_{kl}+{\beta }_{\mathrm{12}}{\sigma }_{im}{\sigma }_{mj}{\sigma }_{kn}{\sigma }_{ml}\end{array}（21）

式中：{\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}},\cdots ,{\beta }_{\mathrm{12}}为12个刚度矩阵参数，它们仅受应力张量的不变量所影响，其他因素对其并无直接作用.\delta为kronecker delta函数，当i=j时，{\delta }_{ij}=\mathrm{1}；当i\ne j时，{\delta }_{ij}=\mathrm{0}.

由于常规试验无法全面确定所有12阶刚度矩阵参数，潘家军等^［38］提出了假设：认为参数{\beta }_{\mathrm{5}}\sim {\beta }_{\mathrm{12}}与高阶应力张量不变量有关，其对于刚度矩阵的影响相对较小，因此可以将参数{\beta }_{\mathrm{5}}\sim {\beta }_{\mathrm{12}}设为0.基于这一假设，式（21）得以简化为：

\begin{array}{l}{\mathit{D}}_{ijkl}\left({\sigma }_{mn}\right)={\beta }_{\mathrm{1}}{\delta }_{ij}{\delta }_{kl}+{\beta }_{\mathrm{2}}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{jk}{\delta }_{il}\right)+\\ {\beta }_{\mathrm{3}}{\sigma }_{ij}{\delta }_{kl}+{\beta }_{\mathrm{4}}{\delta }_{ij}{\sigma }_{kl}\end{array}（22）

将式（22）代入式（19），增量应力张量的表达式可写为：

\mathrm{d}{\mathit{\sigma }}_{ij}={\beta }_{\mathrm{1}}{\delta }_{ij}\mathrm{d}{\epsilon }_{kk}+\mathrm{2}{\beta }_{\mathrm{2}}\mathrm{d}{\epsilon }_{ij}+{\beta }_{\mathrm{3}}{\sigma }_{ij}\mathrm{d}{\epsilon }_{kk}+{\beta }_{\mathrm{4}}{\delta }_{ij}{\sigma }_{ij}\mathrm{d}{\epsilon }_{kl}（23）

对比式（20）和式（23），可以推导出刚度矩阵的各个元素具体表达式为：

D_{\mathrm{1111}}={\beta }_{\mathrm{1}}+\mathrm{2}{\beta }_{\mathrm{2}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{11}}

D_{\mathrm{1122}}=D_{\mathrm{1133}}={\beta }_{\mathrm{1}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{11}}

D_{\mathrm{2211}}=D_{\mathrm{2233}}={\beta }_{\mathrm{1}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{22}}

D_{\mathrm{2222}}={\beta }_{\mathrm{1}}+\mathrm{2}{\beta }_{\mathrm{2}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{22}}

D_{\mathrm{3322}}=D_{\mathrm{3311}}={\beta }_{\mathrm{1}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{33}}

D_{\mathrm{3333}}={\beta }_{\mathrm{1}}+\mathrm{2}{\beta }_{\mathrm{2}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{33}}

D_{\mathrm{1112}}=D_{\mathrm{2212}}=D_{\mathrm{3312}}=-{\beta }_{\mathrm{4}}{\sigma }_{\mathrm{12}}

D_{\mathrm{1123}}=D_{\mathrm{2223}}=D_{\mathrm{3323}}=-{\beta }_{\mathrm{4}}{\sigma }_{\mathrm{23}}

D_{\mathrm{1131}}=D_{\mathrm{2231}}=D_{\mathrm{3331}}=-{\beta }_{\mathrm{4}}{\sigma }_{\mathrm{31}}

D_{\mathrm{1212}}=D_{\mathrm{2323}}=D_{\mathrm{3131}}=\mathrm{2}{\beta }_{\mathrm{2}}

在上述矩阵中，未明确写出的其余元素均默认为0.

3.2 矩阵参数

求解式（6a）、（6b）的逆表达式，可以表示为：

\mathrm{d}p=K_{\mathrm{p}}\mathrm{d}{\epsilon }_{\mathrm{v}}-\frac{K_{\mathrm{p}}G}{K_{\mathrm{q}}}\mathrm{d}{\epsilon }_{\mathrm{s}}（24a）

\mathrm{d}q=G\mathrm{d}{\epsilon }_{\mathrm{s}}（24b）

式中：\mathrm{d}p和\mathrm{d}q分别为土体的平均主应力增量和广义剪应力增量.将增量应力张量的表达式（23）展开，可写为：

\mathrm{d}{\sigma }_{\mathrm{11}}={\beta }_{\mathrm{1}}\mathrm{d}{\epsilon }_{kk}+\mathrm{2}{\beta }_{\mathrm{2}}\mathrm{d}{\epsilon }_{\mathrm{11}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{11}}\mathrm{d}{\epsilon }_{kk}+{\beta }_{\mathrm{4}}\left({\sigma }_{\mathrm{11}}\mathrm{d}{\epsilon }_{\mathrm{11}}+\mathrm{2}{\sigma }_{\mathrm{22}}\mathrm{d}{\epsilon }_{\mathrm{22}}\right)（25a）

\mathrm{d}{\sigma }_{\mathrm{22}}={\beta }_{\mathrm{1}}\mathrm{d}{\epsilon }_{kk}+\mathrm{2}{\beta }_{\mathrm{2}}\mathrm{d}{\epsilon }_{\mathrm{22}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{22}}\mathrm{d}{\epsilon }_{kk}+{\beta }_{\mathrm{4}}\left({\sigma }_{\mathrm{11}}\mathrm{d}{\epsilon }_{\mathrm{11}}+\mathrm{2}{\sigma }_{\mathrm{22}}\mathrm{d}{\epsilon }_{\mathrm{22}}\right)（25b）

\mathrm{d}{\sigma }_{\mathrm{33}}={\beta }_{\mathrm{1}}\mathrm{d}{\epsilon }_{kk}+\mathrm{2}{\beta }_{\mathrm{2}}\mathrm{d}{\epsilon }_{\mathrm{33}}+{\beta }_{\mathrm{3}}{\sigma }_{\mathrm{33}}\mathrm{d}{\epsilon }_{kk}+{\beta }_{\mathrm{4}}\left({\sigma }_{\mathrm{11}}\mathrm{d}{\epsilon }_{\mathrm{11}}+\mathrm{2}{\sigma }_{\mathrm{33}}\mathrm{d}{\epsilon }_{\mathrm{33}}\right)（25c）

在三轴应力条件下，由于\mathrm{d}{\epsilon }_{\mathrm{22}}=\mathrm{d}{\epsilon }_{\mathrm{33}}和{\sigma }_{\mathrm{22}}={\sigma }_{\mathrm{33}}，故\mathrm{d}p和\mathrm{d}q有：

\mathrm{d}p=\frac{\mathrm{d}{\sigma }_{\mathrm{11}}+\mathrm{2}\mathrm{d}{\sigma }_{\mathrm{33}}}{\mathrm{3}}=\left({\beta }_{\mathrm{1}}+\frac{\mathrm{2}}{\mathrm{3}}{\beta }_{\mathrm{2}}+{\beta }_{\mathrm{3}}p+{\beta }_{\mathrm{4}}p\right)\mathrm{d}{\epsilon }_{\mathrm{v}}+{\beta }_{\mathrm{4}}q\mathrm{d}\epsilon（26a）

\mathrm{d}q=\mathrm{d}{\sigma }_{\mathrm{11}}-\mathrm{d}{\sigma }_{\mathrm{33}}={\beta }_{\mathrm{3}}q\mathrm{d}{\epsilon }_{\mathrm{v}}+\mathrm{3}{\beta }_{\mathrm{2}}\mathrm{d}{\epsilon }_{\mathrm{s}}（26b）

联立求解式（24）、式（26），可得各项刚度矩阵表达式的参数分别为：

{\beta }_{\mathrm{1}}=K_{\mathrm{p}}-\frac{\mathrm{2}}{\mathrm{9}}G+\frac{p_{kl}{K}_{\mathrm{p}}G}{q{K}_{\mathrm{q}}}（27a）

{\beta }_{\mathrm{2}}=\frac{G}{\mathrm{3}}（27b）

{\beta }_{\mathrm{3}}=\mathrm{0}（27c）

{\beta }_{\mathrm{4}}=-\frac{K_{\mathrm{p}}G}{q{K}_{\mathrm{q}}}（27d）

针对式（27）参数{\beta }_{\mathrm{1}}中p_{kl}，当下标k=l=\mathrm{1}时，p_{kl}={\sigma }_{\mathrm{22}}+{\sigma }_{\mathrm{33}}-\mathrm{2}{\sigma }_{\mathrm{11}}；当下标k=l=\mathrm{2}时，p_{ij}={\sigma }_{\mathrm{11}}+{\sigma }_{\mathrm{33}}-\mathrm{2}{\sigma }_{\mathrm{22}}；当下标k=l=\mathrm{3}时，p_{ij}={\sigma }_{\mathrm{11}}+{\sigma }_{\mathrm{22}}-\mathrm{2}{\sigma }_{\mathrm{33}}.此时便可以得出次弹性理论的刚度矩阵，从而使非线性剪胀模型在三轴应力空间下得到进一步推广.

4 理论与试验结果的验证

4.1 理论模型参数

针对不同围压条件下两材料的大型三轴剪切试验结果，进行数据处理与分析得出Duncan-Chang模型非线性弹性参数：弹性模量基数K、卸载再加载时的弹性模量基数K_ur、弹性模量指数n、弹性泊松比\mu及破坏比R_f.强度参数：黏聚力c、参考围限压力下的内摩擦角φ₀、内摩擦角随围压的减小量Δφ及残余内摩擦角φ_cr.为了获取Pastor-Zienkiewicz剪胀方程中与剪胀相关的模型参数α，采用式（9）进行反推计算获得，见式（28）.具体的理论模型参数，如表1所示.

\alpha =\frac{\left(\mathrm{3}+M-\frac{q}{p}\right)\mathrm{d}{\epsilon }_{\mathrm{v}}-\mathrm{3}\left(M-\frac{q}{p}\right)\mathrm{d}{\epsilon }_{\mathrm{1}}}{\left(\mathrm{3}\mathrm{d}{\epsilon }_{\mathrm{1}}-\mathrm{d}{\epsilon }_{\mathrm{v}}\right)\left(M-\frac{q}{p}\right)}（28）

其中，\mathrm{d}{\epsilon }_{\mathrm{v}}\mathrm{和}\mathrm{d}{\epsilon }_{\mathrm{1}}分别为大型三轴剪切试验中测量的体积应变与轴向应变.

4.2 模型验证

为了验证所推导非线性剪胀模型的有效性和三轴应力空间下次弹性理论刚度矩阵的适用性，基于ABAQUS有限元软件平台的UMAT子程序，针对具有剪胀特性的材料（如碎石料和面基混合铣刨料），开发出了能在三轴应力状态下运行的非线性剪胀模型计算子程序.在验证该模型和UMAT子程序准确性的过程中，严格模拟了大型三轴剪切试验中的 加载条件、模型尺寸和加载路径，确保模拟结果的精准度和真实性.具体细节包括：计算模型为直径 300 mm、高度600 mm的圆柱体，分析流程遵循先固结后剪切的步骤；随后，运用该UMAT子程序进行模拟计算，并将模拟结果与试验结果进行细致比对，如图4所示.

值得注意的是，该模型在针对具有剪胀特性的粗粒土材料时，展现出广泛的适用性，能够精确预测其强度特征及变形特性.然而，该模型也存在一定的局限性，即该模型主要用于预测粗粒土峰值强度及其之前阶段的应力-应变关系，未能涵盖从峰值强度衰减至残余强度的材料软化过程.具体而言，当粗粒土的强度条件符合Mohr-Coulomb准则时，即标志着其达到峰值强度，此时该模型的UMAT子程序将自动停止计算过程，这意味着该模型对于峰值强度之后的状态无法进行预测.然而，这一限制对实际工程应用的影响相对有限，因为出于结构安全性的首要考虑，工程分析往往仅关注峰值强度及其之前阶段的应力-应变关系.此阶段内土体的力学特性是评估土体结构变形、稳定性以及承载能力的关键依据.

图4的对比结果表明，非线性剪胀模型在模拟不同围压条件下，碎石料和面基混合铣刨料的偏应力-轴应变及体积应变-轴应变关系曲线表现出色.该模型精确地预测了两材料在受力过程中的非线性剪胀特性.与试验数据相比，对于面基混合铣刨料，模型在100 kPa、200 kPa、300 kPa和400 kPa围压下，峰值偏应力的相对误差分别为8.8%、4.6%、3.7%和4.4%，而峰值体积应变的相对误差则在4.7%~16.5%之间波动.对于碎石料，模型在不同围压下的峰值偏应力和峰值体积应变的相对误差均值分别为6.3%和19.1%.这些数据充分证明了模拟计算结果与实验结果的良好一致性，进一步验证了非线性剪胀模型的准确性、刚度矩阵的合理性以及子程序编制的正确性.

此外，鉴于该模型预测试验结果所存在误差的原因，主要归结于以下几个方面：1）土体的各向异性和非均质性.土体本身具有各向异性和非均质性，而该模型假设土体是均质和各向同性的，因此会引发一定的误差；2）试验本身和模型参数确定误差.三轴试验过程中可能存在操作误差，如试样制备、饱和度控制、试验操作等，且模型参数如K、n、\mu、K_ur、R_f、c、φ₀、Δφ、α和φ_cr等，均通过三轴试验确定，但这些参数的测定可能受到试验方法、试验条件等因素的影响，这些因素都可能引入误差，影响模型参数的准确性；3）土体应力路径的多样性.实际土体在受力过程中经历复杂的应力路径，而该模型主要基于轴向应力增加的情况，可能无法准确模拟侧向应力变化下的土体行为.

5 模型参数敏感性分析

选用面基混合铣刨料为例作为参数敏感性分析对象.在进行非线性剪胀模型参数敏感性分析时，考虑到涉及的模型参数较多（共10个参数），难以全面考量所有因素间的交互作用^［39］，为此，首先将全体参数运用方差分析法对其进行初步筛选，以识别出对模型具有显著敏感性的关键参数；随后，针对这些筛选出的关键参数，进一步实施响应面分析，从而评估出模型的参数敏感度.

5.1 全体参数敏感性分析

本节以面基混合铣刨料的强度为例，作参数敏感性分析.在该模型的全体参数的选择上，Pastor- Zienkiewicz模型包含2个参数α和φ_cr，而Duncan-Chang模型涉及8个参数，分别是K、n、\mu、K_ur、R_f、c、φ₀、Δφ.需要注意的是，由于面基混合铣刨料的c为 0 kPa，而卸载再加载时的弹性模量基数K_ur和弹性泊松比\mu均为常数，且经大量计算发现K_ur和\mu对模型强度影响微小，因此不把参数c、K_ur和\mu作为全体参数敏感性分析的讨论对象，则在全局参数敏感性分析中，共需考虑7个参数.

针对K、n、R_f、φ₀、Δφ、φ_cr、α共7个模型参数进行研究因子设计，研究指标为围压100 kPa时所对应的（σ₁-σ₃）和ε₁.在研究因子设计水平中，7个模型参数被视为7个研究因子，以考察其在不同取值下对研究指标的影响.设置低（-1）和高（1）两个水平因子进行研究，低水平设定为模型参数的90%，而高水平则设定为模型参数的1.1倍.各研究因子及水平设计 如表2所示.

为了有效地分析出7个参数对研究指标的敏感性程度，设计出了一个包含14种工况的部分因子数值计算方案.详细方案及计算结果见表3.

在模型全体参数的敏感性分析中采用了方差分析法，这是一种能够系统而全面评估各因子对研究指标影响程度的有效方法，通过该方法能够精确地确定模型参数的显著性，以及因子间交互作用的显著性，方差分析结果见表4.在分析过程中，将P值作为显著性水平的度量，提供了一个概率框架来评估研究因子的交互效应.具体来说，当P值小于0.01时，认为该研究因子具有高度的显著性，并用**进行标记；若P值在0.01~0.05之间，则判定该研究因子显著，用*进行标注；当P值介于0.05和0.1之间时，称该研究因子为边缘显著；而P值大于或等于0.1时，则意味着该研究因子不显著.

从表4中可以看出，全体模型参数对（σ₁-σ₃）的敏感性排序为φ₀、Δφ、K、n、φ_cr、α、R_f，其中φ₀具有高度的显著性，Δφ具有显著性.对ε₁的敏感性排序为R_f、K、φ₀、n、Δφ、φ_cr、α，其中对ε₁影响高度显著的是R_f.综合来看，各参数对（σ₁-σ₃）和ε₁的敏感性排序有较大差异，对于研究指标（σ₁-σ₃），受到参数φ₀、Δφ、K的影响尤为突出，相比之下，参数n、φ_cr、α、R_f的影响则相对较弱；然而，当关注研究指标ε₁时，参数R_f、K、φ₀的影响变得显著，而n、Δφ、φ_cr、α的影响则显得较为有限.

5.2 关键参数响应面分析

通过全体参数敏感性分析，可知各参数对（σ₁-σ₃）和ε₁的敏感性排序存在显著差异，鉴于这种差异性，分别针对（σ₁-σ₃）和ε₁的关键参数开展响应面分析，以更精确地揭示它们之间的相互作用和影响机制.

5.2.1 偏差应力

针对φ₀、Δφ、K这3个关键参数进行响应面分析，以探究它们之间两两交互作用对（σ₁-σ₃）的影响，在此不研究三阶及更高阶的交互效应（影响极小）.每个关键参数设计三个水平：高水平（1）为模型参数的1.1倍，而低水平（-1）则为模型参数的90%.研究因子与水平取值见表5.

根据所设定的水平值制定了研究因子设计方案，并将相应的计算结果汇总在表6中.

根据各试验方案的计算结果计算K、φ₀、Δφ这3个参数以及参数之间二阶交互效应的方差大小.对表6数据结果进行方差分析，偏差应力的方差分析结果如表7所示.可以看出，对偏差应力的敏感性程度由高到低次序为φ₀、φ₀²、Δφ、φ₀Δφ、KΔφ、K、K²、Δφ²、Kφ₀，其中φ₀、φ₀²、Δφ为高度显著，φ₀Δφ为显著.

为了深入剖析参数φ₀、Δφ和K之间交互作用对（σ₁-σ₃）的具体影响程度，通过绘制各参数间响应面图及等值线图（如图5所示）来分析其内在关系.由图5（a）、（b）可知，当保持Δφ=8.35°时，无论参数K处于高水平还是低水平时，（σ₁-σ₃）随着φ₀的增大而增大，并且在这一过程中，等值线的分布也变得更加密集；而在φ₀处于任意水平时，（σ₁-σ₃）几乎不会随着K的增大而发生显著变化，说明φ₀的变化对（σ₁-σ₃）的变化影响更加显著，显著性随着φ₀的减小而减弱，相比之下，K的变化对（σ₁-σ₃）的影响几乎可以忽略不计.

由图5（c）、（d）可知，在保持K=365.59时，当φ₀处于任意水平时，（σ₁-σ₃）几乎不受Δφ变化的影响；而当Δφ处于任意水平时，φ₀的变化对（σ₁-σ₃）的变化影响更加显著，显著性随着φ₀的减小而减弱.由 图5（e）、（f）可知，当保持φ₀=46.47°时，K、Δφ的水平的变化对偏差应力（σ₁-σ₃）的变化影响均不显著，随着K、Δφ的增大或减小，（σ₁-σ₃）基本无变化.

5.2.2 轴向应变

采用响应面分析法探究参数R_f、K和φ₀之间两两交互作用对ε₁的具体影响.同样，设计三个水平，高水平（1）取模型参数的1.1倍，低水平（-1）取模型参数的90%.研究因子与水平取值见表8，研究因子设计方案见表9.

根据研究因子设计方案对各关键参数以及参数之间二阶交互效应的方差进行计算，方差分析结果详见表10.可以看出，对轴应变的敏感性程度由高到低次序为R_f、φ₀、R_f²、K、R_fφ₀、KR_f、φ₀²、Kφ₀、K²，其中R_f、φ₀、R_f²、K、R_fφ₀高度显著，KR_f、φ₀²显著.

鉴于参数R_f、K和φ₀之间对ε₁有显著的交互作用，通过绘制响应面图及等值线图（如图6所示）来深入剖析其内在关系.由图6（a）、（b）可知，当φ₀保持在46.47°时，若K处于低水平，观察到ε₁会随着R_f的增大而增大，且此时等值线分布密集，说明此时R_f的变化对ε₁的影响更为显著，而这种显著性会随着K的减小而逐渐减弱；若R_f处于高水平时，随着K的增大，ε₁呈现减小的趋势，且此时等值线的分布变得相对稀疏，说明此时K的变化对ε₁的影响不显著，显著性随着R_f的减小而减弱.此外，值得注意的是，云图呈现出弧度过渡的形态，而非简单的线性过渡，证明参数R_f和K间存在显著的交互作用，并共同影响着ε₁的变化.

由图6（c）、（d）可知，当K保持在365.59时，当φ₀处于高水平时，R_f的变化对ε₁的影响更为显著，这种显著性随着φ₀的减小而逐渐减弱.然而，在R_f处于低水平的情况下，φ₀的变化对ε₁的影响则相对不显著，并且随着R_f的增大，其显著性反而增强.

由图6（e）、（f）可知，当R_f保持在0.81时，在K处于低水平时，φ₀的变化对ε₁的影响更加显著，显著性随着K的增大而减弱；在φ₀处于高水平时，K的变化对ε₁的影响更加显著，显著性随着φ₀的减小而减弱.K与φ₀之间存在一定的交互作用，且共同影响ε₁的变化.

6 结 论

1）在轴向应变为10%左右时，不同围压下两种材料的强度都逐渐达到残余强度，其中碎石料的残余强度略高出面基混合铣刨料约5%.此外，低围压条件下的残余强度相较于高围压条件要低，差距在25%~50%之间.体积应变与轴向应变之间的关系呈现出一种先减小后增大的趋势，且在一定范围内表现为明显的剪胀特性.

2）在三轴应力状态下，利用ABAQUS有限元软件平台的UMAT子程序对文中非线性剪胀模型进行二次开发，该模型计算结果与试验结果相比展现出高度的吻合度，验证了模型的准确性和有效性.

3）通过全体参数敏感性分析发现，对于研究指标（σ₁-σ₃），受到参数φ₀、Δφ、K的影响尤为突出，相比之下，参数n、φ_cr、α、R_f的影响则相对较弱；然而，对于研究指标ε₁，参数R_f、K、φ₀的影响变得显著，而n、Δφ、φ_cr、α的影响则显得较为有限.

4）通过关键参数及其二阶效应敏感性分析，发现对（σ₁-σ₃）的敏感性程度由高到低次序为φ₀、φ₀²、Δφ、φ₀Δφ、KΔφ、K、K²、Δφ²、Kφ₀，其中φ₀、φ₀²、Δφ为高度显著，φ₀Δφ为显著.对轴应变的敏感性程度由高到低次序为R_f、φ₀、R_f²、K、R_fφ₀、KR_f、φ₀²、Kφ₀、K²，其中R_f、φ₀、R_f²、K、R_fφ₀高度显著，KR_f、φ₀²显著.

5）无论K、Δφ的增大或减小，φ₀的变化对（σ₁-σ₃）的影响更加显著.当K处于低水平时，R_f的变化对ε₁的影响更为显著；当φ₀处于高水平时，R_f的变化和K的变化对ε₁的影响更加显著.

6）文中开发出的面基混合铣刨料的非线性剪胀模型，不仅适用于具有剪胀性岩土材料的数值分析，还能为处理复杂的实际工程问题提供高效、实用的解决方案.另外，该模型有望通过等效刚度和等效强度的方法进一步拓展至土工合成材料加筋粗粒土结构，从而推动加筋土结构数值计算领域的发展.

参考文献

[1] 赵晓华，杨海益，姚莹，等 ．高速公路桥梁路段交通安全风险评价及影响要素挖掘［J］．同济大学学报（自然科学版），2022，50（11）：1637-1646．

[2] ZHAO X H， YANG H Y， YAO Y，et al ．Traffic safety evaluation and feature analysis of freeway bridge sections［J］．Journal of Tongji University （Natural Science），2022，50（11）：1637-1646．（in Chinese）

[3] 杨晓光， 胡仕， 张梦雅 ．智能高速公路交通应用技术发展综述［J］．中国公路学报， 2023， 36（10）： 142-164．

[4] YANG X G， HU S， ZHANG M Y.Development of intelligent motorway traffic application technologies：a review［J］. China Journal of Highway and Transport， 2023， 36（10）： 142-164．（in Chinese）

[5] 裴玉龙， 何庆龄， 侯琳， 等. 基于文献计量分析的交通量预测研究综述［J］. 公路交通科技， 2023， 40（2）： 171-181．

[6] PEI Y L， HE Q L， HOU L，et al ．Review of studies on traffic volume prediction based on bibliometric analysis［J］．Journal of Highway and Transportation Research and Development， 2023， 40（2）： 171-181．（in Chinese）

[7] 陈昌富， 李伟， 张嘉睿， 等 ．山区公路边坡工程智能分析与设计研究进展［J］．湖南大学学报（自然科学版）， 2022， 49（7）：15-31．

[8] CHEN C F， LI W， ZHANG J R，et al ．State-of-the-art of intelligent analysis and design in slope engineering of highways in mountainous areas［J］. Journal of Hunan University （Natural Sciences），2022， 49（7）：15-31．（in Chinese）

[9] YU L T， XIE J， LI R，et al. Study on the performance of emulsified asphalt recycled subgrade based on the evaluation of semi-rigid milling material［J］. Construction and Building Materials，2022，324：126614．

[10] 童海涛，郭知涛，秦凯强 ．高速公路改扩建铣刨料再生利用技术案例分析［J］．公路，2020，65（3）：243-247．

[11] TONG H T， GUO Z T， QIN K Q.Case study on recycling technology of milling materials in expressway reconstruction and expansion［J］．Highway， 2020， 65（3）： 243-247．（in Chinese）

[12] 田源， 牛冬瑜. 铣刨料掺量及掺配方式对再生水稳碎石路用性能的影响［J］. 长安大学学报（自然科学版）， 2020， 40（4）：39-49．

[13] TIAN Y， NIU D Y. Effect of milling material content and blending method on road performance of regeneration cement stabilized gravel［J］. Journal of Chang’an University （Natural Science Edition）， 2020， 40（4）： 39-49．（in Chinese）

[14] VERMA G， KUMAR B ．Prediction of compaction parameters for fine-grained and coarse-grained soils：a review［J］．International Journal of Geotechnical Engineering， 2020， 14（8）： 970-977．

[15] 蒋明杰，吉恩跃，王天成，等 ．粗粒土抗剪强度的缩尺效应规律试验研究［J］．岩土工程学报， 2023， 45（4）： 855-861．

[16] JIANG M J， JI E Y， WANG T C，et al ．Experimental study on laws of scale effects of shear strength of coarse-grained soils［J］．Chinese Journal of Geotechnical Engineering， 2023， 45（4）： 855-861．（in Chinese）

[17] DUNCAN J M， CHANG C Y ．Nonlinear analysis of stress and strain in soils［J］．Journal of the Soil Mechanics and Foundations Division，1970，96（5）：1629-1653．

[18] DUNCAN J M， BYRNE P， WONG K S， et al. Strength， stress-strain and bulk modulus parameters for finite element analysis of stress and movement in soil masses［R］. In： Report No. UCB/GT/80-01， Berkerley： University of California， 1980.

[19] 李广信 ．高等土力学［M］．北京：清华大学出版社，2004．

[20] LI G X ．Advanced soil mechanics［M］. Beijing：Tsinghua University Press，2004．（in Chinese）

[21] DOMASCHUK L， VALLIAPPAN P. Nonlinear settlement analysis by finite element［J］. Journal of the Geotechnical Engineering Division，1975，101（7）：601-614．

[22] NAYLOR D J ．Stress-strain laws for soils［M］. Essex：Applied Science Publishers Ltd.， 1978．

[23] IZUMI N， VERRUIJT A. The dynamic interaction between water flow and soil deformation［J］. Soils and Foundations，1982，22（1）：29-44．

[24] 曾以宁，屈智炯，刘开明，等.土的非线性K-G模型的试验研究［J］. 成都科技大学学报，1985，17（4）：143-149.

[25] ZENG Y N， QU Z J， LIU K M，et al ．Experimental study on nonlinear K-G model of soil ［J］．Journal of Chengdu University of Science and Technology，1985，17（4）：143-149．（in Chinese）

[26] 张丙印，于玉贞，张建民 ．高土石坝的若干关键技术问题［C］//第八届土力学及基础工程学术研讨会论文集. 北京： 中国建筑工业出版社， 2003： 163-186．

[27] ZHANG B Y， YU Y Z， ZHANG J M ．The key questions of high earth-rock dam［C］//Proceedings of the 8th Symposium on Soil Mechanics and Foundation．Beijing：China Architecture and Building Press， 2003： 163-186．（in Chinese）

[28] ROSCOE K H， SCHOFIELD A N， THURAIRAJAH A. Yielding of clays in states wetter than critical［J］. Géotechnique， 1963， 13（3）：211-240．

[29] BURLAND J B. The yielding and dilation of clay［J］.Geotechinique，1965，15（2）： 211-214．

[30] 沈珠江. 理论土力学［M］. 北京： 中国水利水电出版社，2000．

[31] SHEN Z J. Theoretical soil mechanics［M］. Beijing：China Water & Power Press， 2000．（in Chinese）

[32] 殷宗泽 ．一个土体的双屈服面应力-应变模型［J］．岩土工程学报，1988，10（4）：64-71．

[33] YIN Z Z ．A double yield surface stress-strain model of soil［J］. Chinese Journal of Geotechnical Engineering，1988， 10（4）： 64-71．（in Chinese）

[34] 姚仰平 ．UH模型系列研究［J］. 岩土工程学报，2015，37（2）：193-217．

[35] YAO Y P ．Advanced UH models for soils［J］. Chinese Journal of Geotechnical Engineering，2015，37（2）：193-217．（in Chinese）

[36] ZHU E Y， YAO Y P ．Structured UH model for clays［J］.Transportation Geotechnics，2015，3： 68-79．

[37] 程展林，姜景山，丁红顺，等 ．粗粒土非线性剪胀模型研究［J］．岩土工程学报，2010，32（3）：460-467．

[38] CHENG Z L， JIANG J S， DING H S，et al. Nonlinear dilatancy model for coarse-grained soils［J］. Chinese Journal of Geotechnical Engineering，2010，32（3）：460-467．（in Chinese）

[39] 程展林，陈鸥，左永振，等 ．再论粗粒土剪胀性模型［J］．长江科学院院报，2011，28（6）：39-44．

[40] CHENG Z L， CHEN O， ZUO Y Z，et al. Further discussion on dilatancy model for coarse-grained soils［J］. Journal of Yangtze River Scientific Research Institute，2011，28（6）：39-44．（in Chinese）

[41] NOVA R， WOOD D M ．A constitutive model for sand in triaxial compression［J］. International Journal for Numerical and Analytical Methods in Geomechanics，1979，3（3）：255-278．

[42] ROWE P W. The stress-dilatancy relation for static equilibrium of an assembly of particles in contact［J］. Proceedings of the Royal Society of London Series A Mathematical and Physical Sciences，1962， 269（1339）： 500-527．

[43] 徐舜华，郑刚，徐光黎 ．考虑剪切硬化的砂土临界状态本构模型［J］．岩土工程学报，2009， 31（6）： 953-958．

[44] XU S H， ZHENG G， XU G L. Critical state constitutive model of sand with shear hardening［J］. Chinese Journal of Geotechnical Engineering， 2009， 31（6）： 953-958．（in Chinese）

[45] MANZARI M T， DAFALIAS Y F. A critical state two-surface plasticity model for sands［J］. Géotechnique，1997，47（2）：255-272．

[46] GAJO A， WOOD M ．Severn-Trent sand：a kinematic-hardening constitutive model：the q-p formulation［J］．Géotechnique，1999，49（5）： 595-614．

[47] PASTOR M， ZIENKIEWICZ O C， LEUNG K H. Simple model for transient soil loading in earthquake analysis Ⅱ. Non-associative models for sands［J］. International Journal for Numerical and Analytical Methods in Geomechanics，1985，9（5）：477-498．

[48] PASTOR M， ZIENKIEWICZ O C， CHAN A H C. Generalized plasticity and the modelling of soil behaviour［J］. International Journal for Numerical and Analytical Methods in Geomechanics，1990， 14（3）： 151-190．

[49] RAVANBAKHSH E， HAMIDI A. Development of a generalized plasticity constitutive model for cemented sands using critical state concepts［J］. International Journal of Geotechnical Engineering，2013， 7（4）： 364-373．

[50] COLA S， TONNI L. Adapting a generalized plasticity model to reproduce the stress-strain response of silty soils forming the Venice lagoon basin［M］//Solid Mechanics and Its Applications．Dordrecht： Springer Netherlands， 2007： 743-758．

[51] SONG F， JIN Y T， LIU H B，et al ．Analyzing the deformation and failure of geosynthetic-encased granular soil in the triaxial stress condition［J］. Geotextiles and Geomembranes，2020，48（6）：886-896．

[52] HEIDARZADEH H， OLIAEI M ．Development of a generalized model using a new plastic modulus based on bounding surface plasticity［J］．Acta Geotechnica，2018，13（4）：925-941．

[53] MIRA P， TONNI L， PASTOR M，et al ．A generalized midpoint algorithm for the integration of a generalized plasticity model for sands［J］. International Journal for Numerical Methods in Engineering，2009，77（9）：1201-1223．

[54] 潘家军，程展林， 饶锡保，等 ．一种粗粒土非线性剪胀模型的扩展及其验证［J］．岩石力学与工程学报，2014，33（）：4321-4325．

[55] PAN J J， CHENG Z L， RAO X B，et al ．A nonlinear dilatancy model’s extention and validation for coarsed-grained soils［J］．Chinese Journal of Rock Mechanics and Engineering，2014， 33（Sup.2）： 4321-4325．（in Chinese）

[56] 冯亚新，江兆强，孙一清，等 ．基于改进响应面的邓肯-张E-B模型参数全局敏感性分析［J］．水利水电科技进展， 2023， 43（2）：44-50．

[57] FENG Y X， JIANG Z Q， SUN Y Q，et al. Global sensitivity analysis of parameters in Duncan-Chang E-B model based on improved response surface［J］. Advances in Science and Technology of Water Resources，2023，43（2）：44-50．（in Chinese）

基金资助

国家自然科学基金资助项目(52278328)