0.   引言

粗粒土包括堆石料、砾石土、砂砾石料等粒状土，广泛应用于岩土工程的各个领域。在过去较长一段时间内，岩土数值分析中通常假定土体为均质材料；近年来，随着精细化仿真的发展，人们逐渐认识到岩土工程中存在大量材料空间变异性问题^[1-4]，其中，对粗粒土而言，材料的空间变异性主要体现在级配和密度的非均匀性^[2-3]。以土石坝工程为例，传统数值分析中将坝体相同分区的堆石料视为均质材料，而实际上由于施工的不确定性等原因，堆石料的级配和密度（孔隙比）存在显著的非均匀性。刘永涛等^[3]通过施工期坝料检测发现坝体主堆石区的孔隙率约在0.14~0.21，服从正态分布；类似地，朱晟等^[2]发现两河口大坝堆石料级配分形维数约在2.347~2.681，服从正态分布，且进一步通过数值分析指出：若忽略堆石料的空间变异性，会在一定程度上低估坝体的应力和变形，对坝体安全控制是不利的。因此，从精细化模拟的角度出发，需要考虑土体材料的空间变异性。

事实上，在数值分析中考虑土体材料空间变异性的方法主要有两种：一种方法是直接基于经典本构模型，考虑模型参数的空间不确定性，将材料的空间变异性问题转化为模型参数的空间不确定性问题。该方法适用于模型参数较少的工况，例如，在边坡稳定性分析中主要涉及到强度参数c和φ的空间变异性，因此该方法得以广泛应用^[1]。但工程的应力变形分析涉及到的模型参数更多，即便是参数较少的邓肯-张E-B模型，参数也有7个，且各个参数之间并非独立，而是存在耦合关系，采用模型参数不确定性的方法分别生成各个参数的空间随机场^[4]，易忽略参数之间的耦合关系。因此，另外一种方法则需追本溯源，建立土体状态相关本构模型^[2]。粗粒土的空间变异性问题，在宏观层面表现为强度、变形等多个力学参数的空间变异性，但本源是级配和孔隙比这两个物理状态量的空间变异性，直接建立能够考虑粗粒土级配、孔隙比影响的状态相关本构模型是更为合理的思路。

研究表明，粗粒土的力学性质与其状态密切相关，级配、密度和应力水平是影响粗粒土性质的3个主要状态量^[5-8]。从包含或可描述的状态量个数的角度分类，粗粒土的本构模型经历了如下4个阶段的发展。①无状态量：线弹性模型是一种最基本和最简单的本构模型，常用于变形的简单估算，其特点是应力应变在加卸载时呈线性关系，无法反映状态量的影响。②一个状态量：以邓肯-张非线性模型、“南水”弹塑性模型为代表，在国内土石坝工程中应用最为广泛，能够较好地描述粗粒土在不同应力状态下的非线性变形规律，因此，模型可视为包含了应力水平这一个状态量。目前大多数粗粒土本构模型可归为此类^[9]，其特点是对于同一种土料，若初始级配或密度改变，模型的参数都将变化，即这类模型本质上是将相同粗粒土在不同级配、不同密度时视为不同材料。③两个状态量：剑桥模型是基于黏土提出，其重要意义是首次提出了土体临界状态的概念，在此理论基础上，人们针对堆石料、砂土等粒状土的应力变形特性发展出了一系列状态相关本构模型，其中以Li等^[10]提出的砂土二阶状态相关模型为代表，能够同时反映密度和应力水平这两个状态量对粒状土性质的影响，现已广泛应用于堆石料等粗粒土的本构模型之中^[11-14]。④3个状态量：在二阶状态相关本构模型的基础上，进一步将级配引入到粒状土本构模型中，已有较为有益的尝试^{[5, 7, 15]}。其中，典型的方法是将级配以不均匀系数C_u或分形维数D来表征，但验证结果表明当C_u或围压较低时，模型的预测会出现较大偏差^[7]。笔者也曾在模型中以分形维数D来表征级配，但其局限性在于粗粒土的级配需满足分形函数^[6]。

综上所述，构建能综合考虑粗粒土级配、密度和应力水平的三阶状态相关本构模型，既是理论发展的前沿趋势，也是工程实践的迫切需求，目前的研究还远不充分。基于此，本文引入归一化的级配参数来定量表征粗粒土级配的变化，提出了考虑级配、孔隙比和应力水平影响的临界状态统一方程；并结合广义塑性理论，建立了粗粒土三阶状态相关本构模型；最后，利用堆石料在不同级配、不同密实度、不同围压下的三轴固结排水剪切试验，验证了该模型的适用性，为粗粒土工程中考虑材料空间变异性的精细化数值模拟奠定了基础。

1.   临界状态线与等向固结线

1.1   试验方案

蔡正银等^[6]以某土石坝的堆石料作为试验土料，开展了4个级配、4个相对密度以及4个围压条件下共64组大型三轴固结排水试验，试样尺寸为Φ300mm×700 mm，最大粒径为60 mm，试验方案如表 1所示。

表  1  试验方案

Table  1.  Test schemes

试验变量	变量值
级配参数IG	0.163，0.207，0.225，0.305
相对密度Dr	0.60，0.75，0.9，1.0
围压σ3/kPa	300，600，1000，1500

下载: 导出CSV 

| 显示表格

表 1中，土料级配用归一化的级配参数I_G来表征，其定义如图 1所示：I_G=B_I/(B_I+B_P)，其中B_I为级配曲线与d=d_max，d=0.075 mm以及横坐标轴所围成的面积，(B_I+B_P)为d=d_max和d=0.075 mm所围成的面积。

图  1  级配参数I_G的定义

Figure  1.  Definition of gradation index I_G

下载: 全尺寸图片 幻灯片

表 1中的4条初始级配曲线如图 2所示，对应的级配参数I_G分别为0.163，0.207，0.225，0.305，可见，I_G越大，表示土料细颗粒含量越多，粗颗粒含量越少。

图  2  试验初始级配曲线

Figure  2.  Gradation curves in tests

下载: 全尺寸图片 幻灯片

对4个不同级配的土料，设计的试验初始相对密度分别为0.60，0.75，0.90，1.00，换算得到对应的初始孔隙比e₀，如图 3所示。在此基础上，对每种密度的试样分别在300，600，1000，1500 kPa这4种围压下开展三轴固结排水剪切试验。

图  3  不同级配土料初始孔隙比与相对密度的对应关系

Figure  3.  Corresponding relationships between initial void ratio and relative density of soil materials at various gradations

下载: 全尺寸图片 幻灯片

1.2   等向固结点与临界状态点

等向固结点表示的是土体在等向压缩条件下孔隙比与应力的对应关系，对于三轴试验而言，等向固结点可以视为试样固结稳定后、剪切开始前的状态。

以级配I_G=0.207为例，当初始孔隙比e₀=0.287（D_r=0.9）时，试样分别在围压300，600，1000，1500 kPa下达到固结稳定，对应的固结孔隙比e_i分别为0.267，0.254，0.245，0.238；试样固结稳定后继续进行排水剪切，得到的应力应变曲线如图 4所示，当轴向应变达到20%时，剪应力和体变基本稳定，即达到了临界状态点，此时对应的孔隙比为临界状态孔隙比e_c。

图  4  I_G=0.207，e₀=0.287（D_r=0.9）应力应变曲线

Figure  4.  Stress-strain curves of specimens at I_G=0.207 and e₀=0.287

下载: 全尺寸图片 幻灯片

分别将等向固结点和临界状态点绘制于e-${\text{(}}p{\text{/}}{p_{\text{a}}}{\text{)}}$^ξ平面(ξ=0.7），如图 5所示，可以初步得到：当级配I_G与初始孔隙比e₀一定时，随着围压${\sigma _3}$的增加，在e-${\text{(}}p{\text{/}}{p_{\text{a}}}{\text{)}}$^ξ平面，等向固结点和临界状态点的轨迹都能用直线较好地描述，轨迹线分别为等向固结线和临界状态线。在此基础上，进一步探讨当初始级配I_G、初始孔隙比e₀和围压${\sigma _3}$都变化时，如何统一描述e-${\text{(}}p{\text{/}}{p_{\text{a}}}{\text{)}}$^ξ平面内粗粒土的等向固结线和临界状态线。

图  5  等向固结点与临界状态点

Figure  5.  Isotropic consolidation points and critical state points

下载: 全尺寸图片 幻灯片

1.3   临界状态线

研究表明^{[6, 11]}，粗粒土的临界状态线在e-${\text{(}}p{\text{/}}{p_{\text{a}}}{\text{)}}$^ξ平面内可用直线描述，其表达式为

式中：e_c为临界状态孔隙比；e_Γ，${\lambda _{\text{c}}}$和ξ为材料参数；${p_{\text{a}}}$为标准大气压。对于粗粒土，参数ξ一般取为0.7^[6]。

以I_G=0.163（细粒含量最少）和0.305（细粒含量最多）这两个级配土料为例，将堆石料在临界状态时试验值绘制在e-${\text{(}}p{\text{/}}{p_{\text{a}}}{\text{)}}$^ξ平面，如图 6所示。

图  6  堆石料在e-(p/p_a)^ξ的临界状态线

Figure  6.  Critical state lines of rockfill in e-(p/p_a)^ξ space

下载: 全尺寸图片 幻灯片

图 6显示，当级配I_G相同时，不同初始孔隙比e₀的试样各自对应一条临界状态线，且各条线之间基本平行，即斜率${\lambda _{\text{c}}}$与e₀无关，截距e_Γ与e₀相关。当级配I_G不同时，I_G=0.163的斜率${\lambda _{\text{c}}}$为0.0165；I_G=0.305的斜率${\lambda _{\text{c}}}$为0.0123，说明斜率${\lambda _{\text{c}}}$与I_G相关。

一方面，斜率${\lambda _{\text{c}}}$与e₀无关、与I_G相关，将不同初始级配I_G试样的斜率${\lambda _{\text{c}}}$绘制于${\lambda _{\text{c}}}$-I_G平面，如图 7所示，随着I_G的增大，即粗颗粒含量降低、细颗粒含量增加，临界状态线的斜率${\lambda _{\text{c}}}$降低，可表示为

图  7  参数${\lambda _{\text{c}}}$与I_G的关系

Figure  7.  Relationship between parameter ${\lambda _{\text{c}}}$ and I_G

下载: 全尺寸图片 幻灯片

式中：${\lambda _{{\text{c0}}}}$，${\alpha _{{\lambda _{\text{c}}}}}$为材料参数，对于本文堆石料，${\lambda _{{\text{c0}}}}$= 0.0213，${\alpha _{{\lambda _{\text{c}}}}}$=0.0295。

另一方面，截距e_Γ与I_G和e₀都相关。将不同初始级配I_G、不同初始孔隙比e₀试样的截距e_Γ绘制于e_Γ-e₀平面，如图 8所示，可表示为

图  8  参数e_Γ与I_G和e₀的关系

Figure  8.  Relationship among parameters e_Γ, I_G, e₀

下载: 全尺寸图片 幻灯片

式中：${e_{\Gamma 0}}$，${\alpha _\Gamma }$和${\chi _\Gamma }$为材料参数，对于本文堆石料，${e_{\Gamma 0}}$=0.269，${\alpha _\Gamma }$=0.260，${\chi _\Gamma }$=0.602。

将式（3），（2）代入式（1）可得，考虑初始级配I_G、初始孔隙比e₀的临界状态线可统一表示为

1.4   等向固结线

类似地，将等向固结线在e-${\text{(}}p{\text{/}}{p_{\text{a}}}{\text{)}}$^ξ平面内表示为直线，其表达式为

式中：e₀为初始孔隙比；斜率${\lambda _{\text{i}}}$为材料参数。

继续以I_G=0.163，0.305的试验结果为例，将堆石料的等向固结点绘制在e-(p/p_a)^ξ平面，如图 9所示。当级配I_G相同时，不同初始孔隙比e₀的试样对应的等向固结线也基本平行，即斜率${\lambda _{\text{i}}}$与e₀无关。当级配不同I_G时，I_G=0.163的斜率${\lambda _{\text{i}}}$为0.00687；I_G=0.305的斜率${\lambda _{\text{i}}}$为0.0123，说明斜率${\lambda _{\text{i}}}$与I_G相关。

图  9  堆石料在e-(p/p_a)^ξ的等向固结状态线

Figure  9.  Isotropic consolidation lines of rockfill in e-(p/p_a)^ξ space

下载: 全尺寸图片 幻灯片

进一步地，将不同初始级配I_G试样的斜率${\lambda _{\text{i}}}$绘制于${\lambda _{\text{i}}}$-I_G平面，如图 10所示，随着级配参数I_G的增大，等向固结线的斜率${\lambda _{\text{i}}}$降低，可表示为

图  10  参数λ_i与I_G的关系

Figure  10.  Relationship between parameter λ_i and I_G

下载: 全尺寸图片 幻灯片

式中：${\lambda _{{\text{i0}}}}$，${\alpha _{{\lambda _{\text{i}}}}}$为材料参数，对于本文堆石料，${\lambda _{{\text{i0}}}}$= 0.00867，${\alpha _{{\lambda _{\text{i}}}}}$=0.0111。

将式（6）代入式（5）可得，考虑初始级配I_G、初始孔隙比e₀的等向固结线可统一表示为

2.   本构模型构建

2.1   广义塑性模型基本理论

在广义塑性理论中没有屈服函数f、塑性势函数g以及硬化参数h等概念，而是直接定义了流动方向n_g、加载方向n_f和塑性模量H，其刚度矩阵为

式中：${\boldsymbol{D}^{{\text{ep}}}}$为弹塑性刚度矩阵，${\boldsymbol{D}^{\text{e}}}$为弹性刚度矩阵。

流动方向n_g为

由于岩土材料采用不相关联流动法则更为合适，因此广义塑性模型中一般选择的加载方向n_f与塑性流动方向n_g不等，即d_f与d_g不相等，n_f可表示为

2.2   塑性流动方向

式（8）中n_g用来确定塑性流动方向，相当于经典弹塑性理论中的∂g/∂σ，其中g为塑性势函数。对于${n_{\text{g}}}$的确定，广义塑性模型的灵活性体现出了极大的优势，可以通过三轴试验数据直接拟合得到剪胀方程d_g，代入式（9）即可确定流动方向${n_{\text{g}}}$。笔者^[16]此前从考虑颗粒破碎耗能的角度提出了堆石料的剪胀方程，本文在此基础上进一步引入状态参量，得到

式中：d_g=${\text{d}}\varepsilon _{\text{v}}^{\text{p}}/{\text{d}}\varepsilon _{\text{s}}^{\text{p}}$为剪胀比；η=q/p为应力比；M_c为临界状态应力比，本文堆石料M_c=1.72^[3]；n_d为材料参数，为-ln(M_d/M_c)与$\psi$的斜率；β为材料参数，通过式（11）直接拟合d_g-η曲线来确定；$\psi$为状态参量：

剪胀方程在两个边界条件应满足d_g预测值为0。以I_G=0.163，e₀=0.279 (D_r=1.0)，${\sigma _3}$=1500 kPa时的试样为例，如图 11所示，式（11）预测曲线显示，当η=M_c时，e=e_c，则$\psi$=0，d_g=0；当η=M_d时（体变由剪缩变为剪胀的相变点，设该点对应的应力比为M_d），${M_{\text{c}}}\exp ( - {n_{\text{d}}}\psi )$=M_d，则d_g=0。可见，式（11）同时满足上述两个边界条件。

图  11  堆石料的剪胀特性

Figure  11.  Dilatancy behavior of rockfills

下载: 全尺寸图片 幻灯片

2.3   加载方向

式（8）中n_f用来确定加载方向，相当于经典塑性理论中的∂f/∂σ，其中f为屈服函数。在广义塑性模型中，通常直接定义d_f，并基于式（10）来确定n_f。若已知屈服函数f的表达式，式（10）中的d_f也可通过∂f/∂p求取。笔者^[16]提出了一个能够同时描述出粗粒土椭圆、子弹型和水滴形的屈服面函数：

以广泛使用的剑桥模型、修正剑桥模型的屈服函数与式（13）进行对比，同时对文献[17]中堆石料的屈服轨迹试验值开展预测，如图 12所示，可见式（13）对于粗粒土的屈服轨迹描述效果较好。

图  12  典型堆石料屈服轨迹

Figure  12.  Typical yield loci of rockfill

下载: 全尺寸图片 幻灯片

将式（13）中的屈服函数对p求偏导数，即可得到d_f的表达式为

2.4   塑性模量

参考剑桥模型，首先从等向压缩条件出发，即q=dq=0，η=0时，e-lnp平面内，等向压缩塑性体积模量H₀可表示为

式中：$\lambda$，κ分别为e-lnp平面内等向压缩和卸载曲线的斜率。

研究表明^[18]，粗粒土的等向压缩曲线在e-lnp平面不为直线，即加载斜率$\lambda$在e-lnp平面不为定值；卸载曲线在e-lnp平面为直线，卸载斜率κ则在e-lnp平面为定值，如图 13所示。

图  13  e-lnp平面内粗粒土的等向压缩与卸载曲线

Figure  13.  Isotropic compression and unloading curves of rockfill in e-lnp plane

下载: 全尺寸图片 幻灯片

由于本文已基于e-(p/p_a)^ξ平面提出了等向压缩曲线的统一表达式，即式（7）。因此，首先在式（15）的基础上进行坐标变换，将e-lnp坐标系变换为e-(p/p_a)^ξ坐标系，得到等向压缩条件下塑性模量H₀的表达式为

式中：${\lambda _i}$和κ_i分别为等向压缩曲线在e-($p/{p_{\text{a}}}$)^ξ平面的加载斜率和卸载斜率。κ_i在e-($p/{p_{\text{a}}}$)^ξ平面的值可通过e-lnp平面的定值κ进行坐标换算：

由式（7）可知，e_i为固结完成、剪切加载起点的孔隙比，与级配参数I_G和初始孔隙比e₀相关，因此在式（16）的基础上首先引入h₀(1-h_ee_i)这一因式综合考虑级配和孔隙比对塑性模量H的影响。然后，塑性模量H还应满足的边界条件为：当应力比η达到峰值${M_{\text{f}}}$时，塑性模量H为0。基于上述考虑，本文最终确定的塑性体积模量H为

式中：${h_{\text{0}}}$和${h_{\text{e}}}$为材料参数，通过式（16）中对不同固结孔隙比e_i对应的等向压缩塑性体积模量H₀进行最优化拟合确定；${n_{\text{f}}}$为状态相关材料参数，为ln(${M_{\text{f}}}$/${M_{\text{c}}}$)与$\psi$的斜率，${M_{\text{f}}}$为峰值应力比。

2.5   弹性模量

弹塑性模型假定土体应变由弹性应变和塑性应变两部分组成，其中弹性应变可由广义虎克定律确定，需要2个参数，即泊松比$\nu$和弹性模量E。其中，对于弹性部分，粗粒土的泊松比$\nu$一般取0.3左右的常数，根据式（15），弹性体积模量B可表示为

由式（19）可得弹性模量E可表示为

3.   本构模型验证

3.1   模型参数确定

本文模型涉及到的模型参数主要为3类，共14个，如表 2所示：第一类为临界状态参数，共6个；第二类为等向固结/卸载参数，共3个：第三类为状态相关参数，共5个。各参数的物理意义和确定方法如前文所述。

表  2  本文堆石料的模型参数

Table  2.  Model parameters of rockfills in this study

临界状态参数	λc0	αλc	eΓ0	αΓ	χΓ	Mc
参数值	0.0213	0.0295	0.269	0.260	0.602	1.72
等向固结参数	λi0	αλi	κ/10-3
参数值	0.00867	0.0111	6.1
状态相关参数	nd	β	h0	he	nf
参数值	0.748	0.51	1.35	0.98	4.92

下载: 导出CSV 

| 显示表格

3.2   相同I_G相同D_r不同σ₃

对于固定级配、固定相对密度工况，选择三组试验试样，初始级配I_G和相对密度D_r分别为：(a) I_G=0.207，D_r=0.6；(b) I_G=0.305，D_r=0.9；(c) I_G=0.163，D_r=1.0，试验围压${\sigma _3}$均为300，600，900，1500 kPa，如图 14所示。

图  14  相同I_G相同D_r不同${\sigma _3}$三轴试验验证

Figure  14.  Triaxial test verification of specimens at constant I_G constant D_r and various ${\sigma _3}$

下载: 全尺寸图片 幻灯片

图 14表明，在任意初始级配I_G和相对密度D_r条件下，围压${\sigma _3}$越低，应力软化和体变剪胀特征越显著；当围压持续增大时，会逐渐过渡到应力硬化和体变剪缩型，说明粗粒土的性质具有显著的应力相关性。

一般而言，现有模型对于应变硬化和体积剪缩型的应力应变规律都能较好描述，部分模型则无法反映应变软化和体积剪胀现象。图 14中所示的3试样包括了应变硬化和体积剪缩特征（图 14（a）），以及显著的应变软化和体积剪胀特征（图 14（b），（c）），而本文模型都能较好描述。

3.3   相同I_G不同D_r相同${\sigma _3}$

对于固定级配、固定围压工况，选择3试样进行验证，初始级配I_G和围压${\sigma _3}$分别为：(a) I_G=0.163，${\sigma _3}$=600 kPa；(b) I_G=0.255，${\sigma _3}$=900 kPa；(c) I_G=0.305，${\sigma _3}$=1500 kPa，相对密度D_r都为0.60，0.75，0.90和1.0，如图 15所示。

图  15  相同I_G不同D_r相同σ₃三轴试验验证

Figure  15.  Triaxial test verification of specimens at constant I_G various D_r and constant σ₃

下载: 全尺寸图片 幻灯片

当初始级配I_G和围压${\sigma _3}$固定时，相对密度D_r的变化对于应力应变规律的影响甚为显著，说明粗粒土的性质具有显著的密度相关性：随着相对密度D_r的增大，应力由应变硬化型逐渐过渡到应变软化型，体变也由剪缩型过渡到剪胀型。本文模型能较好反映这一演化规律，如图 15所示。

3.4   不同I_G相同D_r相同σ₃

对于固定相对密度、固定围压工况，选择3组试样进行验证，相对密度D_r和围压${\sigma _3}$分别为：(a) D_r=0.75，${\sigma _3}$=600 kPa；(b) D_r=0.90，${\sigma _3}$=900 kPa；(c) D_r=1.0，${\sigma _3}$=1500 kPa，初始级配I_G都为0.163，0.207，0.255和0.305，如图 16所示。

图  16  不同I_G相同D_r相同${\sigma _3}$三轴试验验证

Figure  16.  Triaxial test verification of specimens at various I_G constant D_r and constant ${\sigma _3}$

下载: 全尺寸图片 幻灯片

当相对密度D_r和围压${\sigma _3}$固定时，级配的变化对于应力应变规律的影响也较为显著，说明粗粒土的性质具有显著的级配相关性：随着级配参数I_G的增大（细颗粒含量增多），应力的应变软化特征、体变的剪胀特征越显著，模型能较好反映这一规律，如图 16所示。

综上所述，粗粒土的应力应变具有显著的状态相关性，当围压${\sigma _3}$越低、相对密度D_r越大、级配参数I_G越大（细颗粒含量越多）时，粗粒土会表现出应力软化和体变剪胀特征，反之则表现为应力硬化和体变剪缩。本文模型对于粗粒土在任意级配、任意相对密度和任意围压下的应力应变规律都能较好地预测。

4.   结论

本文基于4种级配、4种相对密度和4种围压组合下的64组堆石料三轴固结排水剪切试验，研究了考虑粗粒土级配、密度和应力水平影响的三阶状态相关本构模型，得到以下3点结论。

（1）粗粒土的临界状态线在e-(p/p_a)^ξ平面内可用直线表示，其截距与初始孔隙比成正相关；与级配参数成负相关，即细颗粒含量越多，截距越小。斜率只与级配参数成负相关，细颗粒含量越多，斜率越小。

（2）粗粒土的等向固结线在e-(p/p_a)^ξ平面内可用直线表示，其截距为初始孔隙比；斜率与级配参数成负相关，即细颗粒含量越多，斜率越小。

（3）粗粒土的应力应变特性具有显著的状态相关性，经典状态相关理论包含密度和应力水平这两个状态量，本文提出了综合考虑级配、密度和应力水平的临界状态线与等向固结线方程，并在此基础上构建了三阶状态相关本构模型。模型验证结果表明，本文模型采用同一套模型参数，就能较好地预测粗粒土在不同级配、不同相对密度和不同围压下的应力应变规律，为粗粒土工程中考虑材料空间变异性的精细化数值模拟奠定了基础。