黄土是第四纪以来外营力作用下形成的沉积土，具有疏松多孔的结构特性^{[1 − 3]}。随着西部城市经济的发展，铁路成为沟通城市的关键^{[4 − 6]}。为解决西北地区多沟壑、地表起伏的问题，多采用压实黄土作为铁路路基填方材料^{[7 − 8]}。然而不同部位黄土压实度不尽相同，在降雨入渗等外界因素的影响下^{[9 − 10]}，压实黄土含水率不断增加，黏聚力及颗粒间摩擦力逐渐减小^{[11 − 12]}，抗剪强度降低^{[13 − 15]}，导致压实黄土在短时间内出现不均匀沉降，极易产生铁路工程隐患^{[16 − 17]}。因此,研究饱和压实黄土的力学特性，对工程建设具有重要意义。

为了描述压实黄土的力学特性，诸多学者先后使用了Duncan-Chang本构模型^{[18 − 19]}、幂函数修正模型^{[20 ]}，虽然这些模型种类繁多，但多适用于应变硬化的土体，而饱和填方黄土在固结不排水试验中表现出明显的应变软化特性，此时传统的双曲线模型或幂函数修正模型已不再适用。而造成应变软化的主要原因是土骨架、土颗粒的破坏与重组，为此，一些学者提出损伤因子的概念^{[21 − 23]}，但是目前多数损伤理论仍以解析理论为主，没有相应的数值算法，进而无法应用于有限元软件的迭代分析。

鉴于此，本文提出修正Duncan-Chang损伤本构模型，采取西安—延安（西延）高铁填方区Qp₃黄土，进行固结不排水室内三轴试验，对试样的初始变形模量与峰值强度进行简要分析。根据修正Duncan-Chang损伤本构模型的张量形式编写出vumat子程序，与三轴试验结果进行验证。最后将vumat子程序应用于填方路基模型中，研究不同压实度与不同填方坡度下的路基瞬时沉降规律。

1.   取样与基本物理参数分析

试样取自陕西关中盆地洛川县西延高铁后子头乡段（图1），该场地地势平缓，线路右侧局部发育小型冲沟。用于填土的${\mathrm{Q}\mathrm{p}}_3^{eol}$黄土呈淡黄色，以粉质颗粒为主、具有大空隙及垂直节理发育。最大填方高度约11.4 m（图2a）。

图  1  研究区域

Figure  1.  Study area

下载: 全尺寸图片 幻灯片

图  2  取样与制样流程

Figure  2.  Sampling and preparation process

下载: 全尺寸图片 幻灯片

取样时首先用直径15 cm、高20 cm的PVC管套入原状土柱，铲刀切削，保鲜膜包裹，防止水分流散，之后使用气泡膜小心缠绕，以防止运输过程中受到干扰（图2b—e）。根据《土工试验方法标准》（GB/T 50123—2019）^[24]进行土的物理性质指标测定，得到试样的基本物理性质参数（表1）。从试验结果看，试样孔隙率高，干密度低，含水率较高，具有典型的黄土特征；使用Battersize2000激光粒度仪对黄土试样进行颗分试验（图2f），结果表明粒径级配曲线有2个峰值，分别位于0.4 μm和40 μm附近。中值粒径$D_{50}$=29.88 μm，平均粒径为40.05 μm，黏粒含量（粒径小于5 μm的颗粒）较少，占比7.5%，粉粒含量（粒径介于2～50 μm的颗粒）较多，占比67.5%，属于粉质黄土（图3）。

表  1  土样的基本物理参数

Table  1.  Basic physical parameters of soil samples

参数	天然含水率/%	比重	天然密度/（g∙cm−3）	干密度/（g∙cm−3）	孔隙比	塑限/%	液限/%
取值	13.6	2.65	1.58	1.45	0.78	17	28

下载: 导出CSV 

| 显示表格

图  3  填方黄土粒径级配

Figure  3.  Particle size gradation of fill loess

下载: 全尺寸图片 幻灯片

2.   三轴剪切试验

2.1   试验仪器

采用英国GDS仪器公司生产的三轴试验系统进行三轴试验。依据《铁路路基设计规范》（TB 10001—2016）^[25]与设计资料，压实度$K$≥0.90，因此，为探究不同压实度填方路基的物理力学特征，三轴试验采用$K$=0.90、0.95和1.00的试验方案，试样尺寸为直径39.8 mm、高80 mm。通过击实试验，确定洛川黄土最优含水率为14%，对应的最大干密度为1.76 g/cm³。定义$K$=（${\rho }_{\mathrm{d}}$/${\rho }_{\mathrm{dmax}}$）×100%，${\rho }_{\mathrm{d}}$为试验黄土的实际干密度，${\rho }_{\mathrm{dmax}}$为最大干密度。不同压实度黄土的制作流程为：（a）将干燥后的黄土彻底粉碎，过筛；（b）将含水率配至14%后闷土72 h待水分均匀扩散；（c）将配置好的土样分4层击实，压实度分别为$K$=0.90（对应${\rho }_{\mathrm{d}}$=1.58 g/cm³）、$K$=0.95（对应${\rho }_{\mathrm{d}}$=1.67 g/cm³）、$K$=1.00（对应${\rho }_{\mathrm{d}}$=1.76 g/cm³）；（d）放入真空饱和缸中饱和24 h（图2g—j）。

由于路堤瞬时沉降过程时间较短，内部水分不能及时排出，因此对试样进行固结不排水试验（consolidated undrain，CU）。结合实际工况，试验中三轴室通过气压控制器施加40，60，100 kPa的围压（${\sigma }_3$），竖向加载柱塞施加轴向压力（${\sigma }_1$）。试件在三轴室固结，直至95%的超孔隙水压力消散；在不排水的情况下，以0.05 mm/min的位移速率剪切试件，直至试件破坏。根据《土工试验方法标准》（GB/T 50123—2019），选取应力-应变曲线上峰值偏差应力作为抗剪强度。

2.2   试验结果分析

2.2.1   应力-应变特征

图4为不同K、${\sigma }_3$下填方黄土的应力-应变关系，反映了黄土的抗剪强度的发展情况。以K、${\sigma }_3$为自变量，填方黄土的应力-应变主要表现出以下特征：（1）相同K下填方黄土的峰值强度${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$随${\sigma }_3$的增加而增加，相同${\sigma }_3$条件下黄土的${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$随着K的增加而增加；（2）$K$=0.90，${\sigma }_3$=40 kPa时，填方黄土应力-应变曲线呈非线性增加，后保持平稳且略有下降；$K$=1.00，${\sigma }_3$=60，100 kPa时，填方黄土在剪切过程中表现出明显的应变软化特性；（3）当${\sigma }_3$从40 kPa增至100 kPa，K=0.90，0.95时，填方黄土对应的${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$增幅分别为125.09%、82.02%；$K$=1.00时，填方黄土对应的${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$增幅为27.70%，说明随K的增加，填方黄土试样${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$的增幅逐渐减小；（4）当${\sigma }_3$从40 kPa增至100 kPa，$K$=0.90时，填方黄土$E_{\mathrm{i}}$的增幅为64.5%；$K$=0.95时，$E_{\mathrm{i}}$增幅为58.84%；当K达到1.00>时，填方黄土$E_{\mathrm{i}}$的增幅达到144.56%，说明随K的增加，黄土试样$E_{\mathrm{i}}$的增幅整体呈递增趋势。

图  4  不同压实度下填方黄土应力-应变曲线

Figure  4.  Stress-strain curves of fill loess under different compactness

下载: 全尺寸图片 幻灯片

2.2.2   初始变形模量与峰值强度

Duncan-Chang本构模型广泛应用于预测及处理路基沉降问题^{[26 − 28]}，而$E_{\mathrm{i}}$与${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$是应用于该本构模型的重要物理力学参数，因此分析参数$E_{\mathrm{i}}$、${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$与$K$、${\sigma }_3$的关系很有必要。基于固结不排水三轴剪切试验的数据，分别以$E_{\mathrm{i}}$与${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$为自变量，$K$与${\sigma }_3$为因变量绘制出三维数值云图如图5（a）、（b）所示。

图  5  压实黄土变形与强度参数三维云图

Figure  5.  Deformation and strength parameters of compacted loess

下载: 全尺寸图片 幻灯片

由图可以看出，$E_{\mathrm{i}}$与${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$随$K$、${\sigma }_3$的增加而增加，但由三维云图投影在XY面等值线的外法线指向不同，可以判断出$E_{\mathrm{i}}$与${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$在某一点（${\lambda }_{\mathrm{c}}$, ${\sigma }_3$）的梯度$({\displaystyle \frac{\partial E_{\mathrm{i}}}{\partial K}},{\displaystyle \frac{\partial E{}_{\mathrm{i}}}{\partial {\sigma }_3}})$与$({\displaystyle \frac{\partial {({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}}{\partial K}},{\displaystyle \frac{\partial {({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}}{\partial {\sigma }_3}})$方向相反，$E_{\mathrm{i}}$的梯度指向为图5（a）中XY平面左向紫色区域，即压实度和围压减小方向；${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$的梯度指向为图5（b）中XY平面红色区域，即压实度和围压增加的方向。因此，随着$K$与${\sigma }_3$的增加，$E_{\mathrm{i}}$的变化率逐渐减小，而${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$的变化率逐渐增加。

为了进一步量化不同土样的$E_{\mathrm{i}}$与${({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}$，根据图5（a）（b）进行非线性三维表面拟合（图5cd），得到拟合式（1）拟合优度$R^2$分别为0.94和0.98，因此式（1）可以较为准确地表述初始变形模量、峰值强度与压实度、围压之间关系。

3.   模型构建及数值实现

3.1   统计损伤模型建立

Duncan-Chang本构模型多用于描述路基土的应力-应变关系，但是该本构模型仅适用于应变硬化土样，而路基瞬时沉降时，饱和填方黄土土样表现为应变软化特征，因此需要对Duncan-Chang模型进行修正。由于压实黄土受压过程中，土骨架破坏重组，抗剪能力减弱，实质上反映出土体的损伤过程^{[29, 30]}。故本文引入损伤变量对Duncan-Chang模型进行修正，便于预测不同围压下饱和压实黄土的应力应变状态。假设峰值破坏前，土体试样的应力应变满足Duncan-Chang本构模型：

式中：${\epsilon }_1$——轴向应变/%；

a——参数，$a={\displaystyle \frac{1}{E{}_{\mathrm{i}}}}$；

b——参数，$b={{\displaystyle \frac{1}{({\sigma }_1-{\sigma }_3)}}}_{\mathrm{ult}}$；

$({\sigma }_1-{\sigma }_3{)}_{\mathrm{ult}}$——极限偏差应力/kPa。

在固结不排水三轴剪切试验中，由于$\mathrm{d}{\sigma }_2$= $\mathrm{d}{\sigma }_3$=0，所以切线模量（$E_{\mathrm{t}}$）可以表示为：

由于$\lg \left(E_{\mathrm{i}}/p_{\mathrm{a}}\right)$与$\lg \left({\sigma }_3/p_{\mathrm{a}}\right)$两者近似呈线性关系^[31]，所以$E_i$也可表示为$E_{\mathrm{i}}=T{p}_{\mathrm{a}}{\left({\displaystyle \frac{{\sigma }_3}{p_{\mathrm{a}}}}\right)}^n$，$p_{\mathrm{a}}$为大气压（101.4 kPa），$T、C$为试验常数^[32]。定义破坏比$R_{\mathrm{f}}={\displaystyle \frac{{({\sigma }_1-{\sigma }_3)}_{\mathrm{f}}}{{({\sigma }_1-{\sigma }_3)}_{\mathrm{ult}}}}$，可得：

由于采用Mohr-Coulomb屈服准则，即：

式中：$c、\phi$——土体的黏聚力/kPa、内摩擦角/（°）。

由式（4）（5）得：

当土体满足Mohr-Coulomb准则进入屈服阶段后采用损伤理论进行修正，根据等效应变假设：

式中：$\sigma 、\epsilon$——名义应力、名义应变；

${\sigma }^{\ast }、{\epsilon }^{\ast }$——等效应力、等效应变；

$D$——损伤变量。

根据连续介质力学，损伤变量$D$可定义为某一荷载作用下，土体内已损伤破坏的微元个数$N_{\mathrm{t}}$与未损伤破坏微元个数$N$的比值：

假设连续随机破坏的微元体服从Weibull分布，其破坏概率密度函数为：

式中：$F$——屈服准则；

$m$、$F_0$——Weibull分布参数。

根据破坏的概率密度函数和损伤变量的定义，由式（7）可得到土的统计损伤变量：

由于室内三轴剪切试验测得的应力为名义应力，故依据广义胡克定律和式（7）（8）可得：

式中：ν——泊松比。

Mohr-Coulomb准则可由双剪统一强度理论退化而来，故根据李杭州等^[33]的研究，土的统计损伤变量可表示为：

其中：

将式（7）（8）带入式（2）中可得：

进而式（3）可变为：

为方便后续的二次开发验证，假设式（14）中$E_{\mathrm{t}}$为固结不排水（consolidation undrained，CU）试验峰值点（$E_{\mathrm{f}}$）的左斜率，即$E_{\mathrm{t}}=E_{{\mathrm{f}}^{-}}$。根据广义胡克定律，及式（16）得到饱和压实黄土损伤本构模型的最大主应力形式：

3.2   参数求解

根据上文所述，服从Weibull分布损伤的Duncan-Chang本构模型需$T$、$n$、$R_{\mathrm{f}}$、$c$、$\phi$、$m$、$F_0$共7个参数。其中$m$、$F_0$可以通过对三轴试验数据进行线性化处理获取：

其中：

对试验数据进行拟合可获得参数$m$和$A$，进而根据式（21）求取参数$F_0$, $T$、$n$、$R_{\mathrm{f}}$、$c$、$\phi$等材料参数可结合文献[34]将CU三轴试验数据进行处理后确定。

4.   模型验证

4.1   常规三轴压缩算例

ABAQUS在处理未内置的本构模型时，提供了隐式和显示两种二次开发的子程序接口，由于路堤的瞬时沉降时间较短，符合vumat显示计算的工况条件，故本文采用vumat子程序进行验证。

vumat子程序是ABAQUS提供给用户进行ABAQUS/Explicit材料本构模型二次开发的用户子程序接口，其核心任务是对单元的每个积分点进行以本构模型增量形式为基础的算法迭代，进而实现某一时刻应力${\sigma }_{ij}^n$向下一分析步${\sigma }_{ij}^{n+1}$的更新。

根据广义胡克定律，填方黄土试样破坏前，其增量形式为：

式中：$G、\lambda$——剪切模量、拉梅常数；

${\delta }_{ij}$——克罗内克函数。

当试样应力关系满足Mohr-Coulomb准则时，破坏产生，且损伤也随之进行，此时模型增量形式变为：

式中：$D_{\mathrm{e}}^{\mathrm{d}}$——带有损伤修正的刚度矩阵。

根据式（22）（23）以及3.1节所述内容，使用Fortran语言编写vumat，具体流程如图6所示。首先在ABAQUS材料模块完成用于vumat的材料参数定义，之后依据ABAQUS操作流程，来到作业模块，完成双精度设置，调用vumat子程序文件。在初始分析步时，假设本构满足初始变形模量为$E{}_{\mathrm{i}}$的线弹性迭代，之后在下一分析步中通过Duncan-Chang本构模型由式（6）对变形模量$E$进行修正，在这一修正阶段，模型的应力更新根据式（23）仍然服从弹性迭代，在进行第$n$+1步的应力${\sigma }_{n+1}$更新时，若主应力状态满足Mohr-Coulomb破坏准则，模型会进入损伤阶段，由式（16）通过损伤变量$D$更新变形模量$E_{n+1}$，进而完成应力${\sigma }_{n+1}$的更新。以此流程对每个积分点进行计算，并将每一分析步更新的应力状态返回至ABAQUS主程序中，最终生成odb结果文件。

图  6  vumat运算流程图

Figure  6.  Flow chart of vumat operation

下载: 全尺寸图片 幻灯片

经过对西延高铁填方区黄土试样的CU试验数据进行处理，得到不同围压下压实度K=0.90、0.95、1.00时的模型参数（表2）。

表  2  模型参数

Table  2.  Model parameters

参数	K	σ3/kPa	$T$	$n$	$R_{\mathrm{f}}$	Weibull参数		c/kPa	φ/（°）
						m	$F_0$
取值	0.90	40	235.962	0.558	0.952	0.262	9.211	3.70	15
		60			0.879	0.266	12.279
		100			0.848	0.400	103.073
	0.95	40	302.924	0.522	0.830	0.220	4.307	12.90	16
		60			0.853	0.339	44.517
		100			0.865	0.405	250.158
	1.00	40	407.532	0.978	0.798	0.369	89.871	17.30	18
		60			0.819	0.341	61.321
		100			0.853	0.408	135.915

下载: 导出CSV 

| 显示表格

本文建立单位长度有限元模型，边界条件与试验一致，网格采用C3D8网格，依据CU试验分析步设置为两个阶段：①施加围压的预加载阶段；②施加轴向偏差应力的加载阶段。基于建立的三维有限元模型和表2 的填方黄土物理力学参数，本文模拟了围压下不同压实度饱和路基黄土的CU试验，并将试验结果、Duncan-Chang（图中简称“D-C”）模型结果及本文模型结果进行对比，如图7所示。数值模拟的结果可以较为准确地反映西延高铁路堤填方黄土破坏前应力-应变轨迹，破坏时的峰值强度及其所对应的应变，而破坏后根据损伤修正的Duncan-Chang本构模型，虽然模拟结果与试验有些差异，但揭示了剪切过程中填方黄土应变软化的力学特性，在路基沉降分析中具有重要意义。

图  7  不同围压、不同压实度模型结果对比

Figure  7.  Comparison of results of different compacting models under different confining pressures

下载: 全尺寸图片 幻灯片

4.2   填方路基分析

考虑到填方坡度与黄土压实度都会对路基的瞬时沉降量产生影响，本文建立起如图8所示的填方路基模型，模型厚度为1，并且由路堤、原始地基和填方黄土三部分组成，其中路堤长15.0 m、高4.6 m，原始地基长50 m、深度25 m。根据《建筑地基基础设计规范》（GB 5007—2011）^[35]填方坡度一般为1∶0.3～1∶1.25，故计算坡脚（α）设置为45°、60°及73°，填方黄土深度为10 m，上底长25 m，下底与上述坡脚相对应，分别为15.0，19.2，21.9 m；压实度分别为0.90，0.95，1.00。为便于计算，路堤整体采用ABAQUS内置的弹性本构；原始地基采用Mohr-Coulomb弹塑性本构；填方黄土采用本文提出的修正Duncan-Chang损伤本构模型，各部分材料参数见表3。模型边界条件为：底面固定，左右两边在X方向进行约束，对路基整体施加填方路基竖直方向上的瞬时沉降。

图  8  填方路基模型图

Figure  8.  Fill subgrade model drawing

下载: 全尺寸图片 幻灯片

表  3  路基材料参数

Table  3.  Subgrade material parameter

材料	$E_{\mathrm{i}}$/MPa	E/MPa	ν	K	T	n	Rf	Weibull参数		c/kPa	φ/（°）
								m	$F_0$
填方黄土	15.93	—	0.30	0.90	235.962	0.558	0.952	0.262	9.211	3.70	15
	20.39			0.95	302.924	0.522	0.830	0.220	4.307	12.90	16
	24.37			1.00	407.532	0.978	0.798	0.369	89.871	17.30	18
原始地基	—	16.15	0.30	—	—	—	—	—	—	20.50	12
路堤	—	50.00	0.25	—	—	—	—	—	—	—	—

下载: 导出CSV 

| 显示表格

计算结果如图9所示，由竖直方向（Y方向）位移云图可以看出，当填方坡度增加时，路堤右侧靠近填方区的部分沉降量明显增加。

图  9  不同填方坡度路基沉降图

Figure  9.  Settlement of the subgrade with different fill slope

下载: 全尺寸图片 幻灯片

将填方黄土的材料属性替换为ABAQUS内置的Mohr-Coulomb弹塑性本构，计算路堤瞬时沉降量。限于篇幅，本文仅在前后两种不同材料本构下，提取路堤顶部中点的沉降量进行对比，如图10所示。可以看出，路基的沉降量随填方黄土压实度的增加而减小，随填方坡度的增加而增加。不同的是ABAQUS内置的Mohr-Coulomb本构模型所得瞬时沉降偏小，计算结果较为保守，且对填方坡度的变化不如本文损伤本构模型更加敏感。

图  10  填方路基沉降量与填方坡度、压实度的关系曲线

Figure  10.  Relationship between the settlement of fill subgrade and the slope and compaction degree of fill

下载: 全尺寸图片 幻灯片

5.   结论

（1）不排水三轴剪切试验表明，饱和压实黄土表现出应变软化特性，峰值强度与初始变形模量随压实度的增加而增加。峰值强度随围压增大而增大，压实度低时峰值强度增幅达125.09%，压实度高时峰值强度增幅达27.7%。初始变形模量随围压增大而增大，压实度低时初始变形模量增幅为64.5%，压实度高时初始变形模量增幅达144.6%。

（2）建立了一种饱和压实黄土的Duncan-Chang损伤本构模型，开发了该模型的vumat子程序，通过ABAQUS软件验证了该模型正确性，该模型能够较好地描述饱和压实黄土受剪时应变软化特性，克服了Duncan-Chang本构模型难以反映应变软化的局限性，适用于黄土填方路基瞬时沉降计算分析。

（3）基于本文Duncan-Chang损伤本构模型与Mohr-Coulomb弹塑性本构模型的数值分析对比，揭示了填方坡度和压实度对填方路基瞬时沉降量的影响规律，发现路面沉降量随填方路基压实度增加而逐渐减少，而随填方坡度增加而逐渐增加；当填方坡度由60°增至73°时，路基沉降量最为显著。结果可为黄土填方路基沉降分析提供参考。

值得注意的是，本文虽然建立了一种黄土Duncan-Chang损伤本构模型，采用陕西洛川粉黄土进行了试验与计算验证，可为黄土地区填方工程沉降分析提供参考，但由于黏黄土、粉黄土和砂黄土的土性存在一定差异，模型的适用性有待后续深入研究。