基于亚塑性本构模型的碎石路基渗透变形研究

孙子元，唐 扬，赵建永，樊 璇

(1.中国电建集团华东勘测设计研究院有限公司 浙江华东工程建设管理有限公司， 杭州 311122; 2.中南建筑设计院股份有限公司 工程数字技术中心，武汉 430071； 3. 长江科学院 水利部岩土力学与工程重点实验室，武汉 430010)

1 研究背景

路基作为线状工程中重要的组成部分，既承受着上部结构的荷载，同时还要承受交通工具在行驶时对路基表面上的动应力，故路基填料的力学特性对整个道路而言至关重要。碎石土具有抗剪强度高、抗变形能力好等特点，被广泛用于路基工程建设中[1-3]。但该种土体由于颗粒间的黏聚力小、孔隙较大，在水流的作用下容易出现渗透变形。一般来说，渗透作用对路基填料的不利影响主要表现为2个方面：一方面是导致路基边坡的失稳[4-6]；另一方面，渗透作用会改变填料内部结构组成，进而使填料本身发生渗透破坏[7-8]。目前，针对后者的研究还相对较少，因而，如何正确评价碎石路基在渗透作用之后的变形及强度变化已成为研究其水稳定性的关键问题之一，应该得到广泛重视。

部分学者[9-11]以三峡库区为研究区域，针对库区内库水位涨落导致碎石土强度参数弱化的问题进行研究，认为水对碎石土的弱化作用受到含水率、含石量和密度等因素的控制；同时，还提出了不同碎石含量下，抗剪强度参数随着含水率变化的弱化公式。此外，还有部分学者[12-14]采用有限元计算方法，结合Mohr-Coulomb强度准则来判断路堤边坡在渗流作用下的稳定性。但上述研究均没有考虑渗透作用对填料内部结构的改变。事实上，碎石土的强度参数，尤其是黏聚力，会因粒度成分的不同而出现很大差别[15]。此外，大多数本构模型由于表达式复杂、参数较多，难以进行实际应用。因此，本文针对碎石路基在渗透作用下的内部结构变化展开研究，采用室内试验的手段，探明渗透作用前后碎石土的粒度改变情况，以及变形和强度参数的变化情况。采用亚塑性本构模型嵌入ABAQUS软件中，描述碎石填料在颗粒流失前后的应力-应变关系，计算其颗粒流失前后的沉降变形。其结果可以为碎石路基的工程设计及渗透变形计算提供相关参考。

2 试样及试验方案

本次试验研究选用水布垭灰岩碎石料，其颗粒级配如表1、图1所示。试验开始前，通过重型击实试验确定了碎石料的最大干密度为2.271 g/cm3和最优含水量约为8%(图2)。考虑到实际工程的一般碾压情况，依照最大干密度指标，按照85%、88%和91%三个压实度开展后续的渗透及力学试验，对应干密度分别为1.93、1.99、2.05 g/cm3。

表1 灰岩碎石料试验颗粒级配Table 1 Test gradations of crushed rock aggregate of limestone

图1 灰岩碎石料试验颗粒级配曲线Fig.1 Test gradation curve of crushed rock aggregate of limestone

图2 灰岩碎石料击实试验曲线Fig.2 Test compaction curve of crushed rock aggregate of limestone

本次试验主要为了研究渗透作用前后碎石土的粒度改变情况，以及变形和强度参数的变化情况。首先对3种不同的碎石料进行渗透变形试验，获得试样在渗透作用下的颗粒流失情况，以及颗粒级配的变化情况；之后，再根据试验结果，对颗粒流失前后的碎石料分别进行三轴固结排水试验(CD)，获得其强度参数，以及应力-应变关系的变化情况。

3 渗透变形试验及大型三轴试验

3.1 渗透变形试验

渗透变形试验依照《土工试验方法标准》(GB/T 50123—2019)[16]进行。试验仪器如图3所示。按表1试验颗粒级配计算各粒径区间碎石料质量并准确称量，加入约2%水量并拌合均匀。将事先制备好的试样按相应的密度，均匀分2层装填在渗透仪中，每层高度为10 cm。各层试样均采用表面振动器进行振动击实。

图3 垂直渗透变形装置示意图Fig.3 Schematic diagram of vertical penetration deformation device

试验开始前，需对试样进行饱和。调整水箱水位略高于试样底面位置，再缓慢提升水箱至一定高度，待水箱水位与试样中水位相等，并停留一段时间后，再提升水箱水位。随着水箱水位上升，让水由试样底部向上渗入，使试样缓慢饱和，以排除试样中的空气。

试验时提升供水桶，使供水桶的水面高出渗透仪的溢水口(上进水口)，保持常水头差，形成渗透坡降，测记测压管水位，并用量筒测量渗水量。若连续3次测得的渗水量基本稳定，又无异常现象，即可提升至下一级水头。试验中同时测读水温、气温，并观测记录试验过程中出现的各种现象，如水的浑浊程度，冒气泡，细颗粒的跳动，移动或被水流带出，土体悬浮、渗流，以及测压管水位的变化等。

3.2 大型三轴试验

三轴试验采用长江水利委员会长江科学院水利部岩土力学与工程重点实验室的应力-应变式高压大型静力三轴仪[16-18](图4)，试样尺寸为Ф300 mm×600 mm，最大围压3.0 MPa，最大轴向应力21 MPa，最大行程为300 mm。设备主体包括刚性反力架、轴向油压千斤顶、三轴压缩室，以及稳压系统(采用水气转换稳压系统稳定微调试验过程中的围压)。

图4 大型静力三轴试验仪Fig.4 Large-scale static triaxial test equipment

试验过程严格按照《土工试验方法标准》(GB/T 50123—2019)[19]进行，根据试验的要求进行级配配置，三轴试验中将试样分4层装入成模筒，然后振捣密实，再抽真空静置一段时间，如无漏气则安装三轴仪压力室，进行充水饱和，试样准备完成之后便可开展三轴固结排水剪切试验(CD)。

4 试验结果分析

4.1 渗透作用试验结果

渗透变形的结果是部分细颗粒流失、试样整体质量下降、颗粒级配发生变化，这也是后续开展三轴试验和压缩试验进行对比的重要依据。

图5是渗透试验前后不同密度的试样颗粒级配的变化情况。试验结束后，以<5 mm的颗粒作为细颗粒，计算并分析颗粒流失的情况：密度为1.93 g/cm3的试样，试验最大水力比降为2.75，细颗粒流失量占总质量的3.78%，其中，粒径在[0.25，0.1] mm颗粒流失量占原质量的36.6%，粒径<0.1 mm颗粒流失量占原质量的76.5%。密度为1.99 g/cm3的试样，试验最大水力比降为3.14，细颗粒流失量占总质量的5.61%，其中，在[0.25，0.1] mm颗粒流失量占原质量的50.4%，粒径<0.1 mm颗粒流失量占原质量的91.7%。密度为2.05 g/cm3的试样，最大水力比降为4.26，细颗粒流失量占总质量的2.06%，其中，粒径<0.1 mm颗粒流失量占原质量的53.6%。

图5 不同密度条件下渗透试验前后的颗粒级配变化Fig.5 Change of particle gradation before and after permeability test under different density conditions

本次室内渗透变形试验主要是为了获得碎石土在发生渗透变形下的颗粒流失情况，故试样在室内试验条件允许的情况下，尽可能地获得最大的细颗粒流失情况，所以3个试样试验结束时的水力比降有所差异，对最终的流失量有所影响。从试验结果来看，密度为1.99 g/cm3的试样，在渗透破坏的时候颗粒流失率最高。且排除试验击实及加压固结过程中的颗粒破碎，颗粒流失主要发生在<2 mm粒径区间。此外，由于干密度2.05 g/cm3在渗透变形试验下出现完全破坏，故不再考虑其颗粒流失之后的强度变化。

4.2 大型三轴试验结果

根据不同围压条件下的应力-应变曲线可以绘制相应的一组莫尔应力圆，进而获得相应的细颗粒流失前后路基填料的线性剪切强度指标。由于通过一组莫尔应力圆的切线来获得抗剪强度参数存在一定的困难性，故我们主要采用公式计算方法，结合最小二乘法拟合法，得到强度参数。

图6是不同密度的2种试样在颗粒流失前后抗剪强度参数的变化情况。试样密度为1.93 g/cm3时，细颗粒流失前后试样的内摩擦角均为39.5°左右，基本相同，但黏聚力从流失前的47.576 kPa降低到流失后的29.867 kPa；试样密度为1.99 g/cm3时，细颗粒流失前后试样的内摩擦角均为40.0°左右，基本相同，但黏聚力从流失前的44.056 kPa降低到流失后的28.790 kPa。

图6 颗粒流失前后莫尔圆及强度参数Fig.6 Mohr’s circles and strength parameters before and after particle loss

碎石料在经历细颗粒流失之后，黏聚力有明显的减小，但内摩擦角几乎无变化，表明细颗粒流失对碎石路基粗粒组骨架的受力影响不大，主要影响体现在粗细颗粒的黏结力上。

从应力-应变曲线(图7)来看，不同条件下的应力-应变关系基本呈现应力强化型，即当轴向应变进一步增大时，偏应力无明显峰值，而是始终保持在接近峰值的状态，未出现明显的应变软化现象，表明该种碎石料的工程特性较好。但从曲线的初始阶段来看，细颗粒流失会对变形特性产生影响，细颗粒流失之前，曲线的初始阶段较陡，而细颗粒流失以后，曲线的初始阶段较为平缓，即初始的切线模量在颗粒流失之后会变小。

图7 颗粒流失前后不同围压条件下试样的 应力-应变关系Fig.7 Stress-strain relations under different confining pressures before and after particle loss

5 亚塑性本构模型及其数值实现

5.1 亚塑性本构模型简介

亚塑性理论是20世纪中晚期诞生的一类全新本构理论，该理论以连续介质力学为基础，采用张量工具，摒弃了传统弹塑性力学中人为的假定概念，如屈服面、硬化规律、塑性势及流动法则等，直接建立应力率和应变率之间的张量函数关系。

Wu等[20]根据亚塑性理论的基本原理，针对颗粒材料提出了一种亚塑性本构方程，其基本方程为

(1)

(2)

式中δij是二阶单位张量。

式(1)等号右边的前3项用来描述材料可逆的表现，后面1项则用来描述材料的不可逆表现以及卸载过程中刚度的转变。该方程最突出的优点就是形式简单，参数易获取，可以很好地反映土体的非线性力学行为。

如果需要考虑黏聚力对于材料强度的影响，我们可以在Cauchy应力上增加一个结构张量，即将式(1)中的Cauchy应力张量替换为σc，且按如下表达式计算，即

σc=σ-pcδij。

(3)

式中pc是与黏聚力c和内摩擦角φ相关的参数，表达式为

(4)

因此，Wu亚塑性本构方程可以重新表达为

(5)

5.2 亚塑性本构模型在ABAQUS软件中的实现

ABAQUS是一款功能强大的通用有限元软件，由于岩土工程问题涉及较多的非线性问题，而该软件在这方面有较为突出的贡献，故该软件非常适合用于分析岩土工程问题。尽管ABAQUS中的材料模型非常丰富，但仍有一些岩土工程常用的本构模型未包含在内，限制了ABAQUS在岩土工程的应用。但该软件提供给了方便的二次开发接口(UMAT)，可以让用户灵活地创建自定义材料模型。

在每个增量步的每一次迭代中，在每个积分点处UMAT均会被调用，以获取材料信息，形成刚度矩阵，求解应力和相关的一些变量。由于这样频繁地调用用户子程序，故我们需要通过一些积分方法来提高计算的效率和准确性。本研究使用FORTRAN语言，以改进欧拉法作为积分算法，编写亚塑性本构模型的数值实现代码[21]，由于篇幅有限，不再给出程序代码。

6 算例验证分析

6.1 三轴压缩试验数值计算

为了验证模型二次开发的正确性以及在路基碎石料中的适用性，采用亚塑性本构模型进行三轴试验的数值计算，并将试验数据与数值计算结果进行对比分析。亚塑性本构模型计算参数如表2所示，内摩擦角和黏聚力根据上述三轴试验结果取值。

表2 数值模拟计算参数Table 2 Numerical simulation parameters

图8是基于亚塑性本构模型的数值试验计算曲线与实际试验数据之间的对比结果。从图8来看，亚塑性本构模型反映了土体的非线性力学行为，且数值模拟计算的三轴压缩曲线和试验数据基本吻合，证明该模型能很好地描述碎石料在不同围压下的应力-应变关系。

图8 不同围压下颗粒流失前后数值计算应力-应变 曲线与试验数据对比Fig.8 Comparison of stress-strain relations between test and numerical simulation before and after particle loss under different confining pressures

6.2 碎石路基数值模拟

本研究主要关注碎石路基运行期间由于渗透作用造成的颗粒级配变化对其变形的影响，因此，计算时不考虑孔隙水压力。数值计算模型及网格剖分如图9所示，边坡的坡度为1∶1.67，模型左右边界和下边界固定，模型顶部施加一个30 kPa均匀分布的荷载，用来模拟路基上的行人和车辆。计算参数通过室内三轴试验获取，路基填料的物理力学参数如表3所示，其对应的亚塑性本构模型参数如表2所示。

图9 碎石路基数值模拟计算模型网格Fig.9 Meshes of numerical simulation model of gravel subgrade

表3 路基填料的物理力学参数Table 3 Physical and mechanical parameters of subgrade soil material

图10、图11是不同密度碎石路基在颗粒流失前后的竖向变形云图。

图10 密度1.93 g/cm3的碎石路基颗粒流失前后 运行期竖向位移云图Fig.10 Vertical displacement distribution of gravel subgrade before and after particle loss with density of 1.93 g/cm3

图11 密度1.99 g/cm3的碎石路基颗粒流失前后 运行期竖向位移云图Fig.11 Vertical displacement distribution of gravel subgrade before and after particle loss with density of 1.99 g/cm3

从图10、图11可以看出，变形较大的部位主要集中在远离坡面的位置。对比颗粒流失前后的路基变形情况，发现颗粒流失之后，碎石路基的变形沉降会出现增大的趋势，对于密度为1.93 g/cm3的碎石路基，颗粒流失之前的最大竖向变形为1.24 cm，而颗粒流失之后，最大竖向变形增大到2.07 cm。对于密度为1.99 g/cm3的碎石路基，颗粒流失之前的最大竖向变形为1.1 cm，而颗粒流失之后，最大竖向变形增大到1.79 cm。

图12是颗粒流失前后碎石路基运行期地表不同位置的沉降曲线。从图12可以看出，路基表面在受荷载作用时，路基中心的沉降较大，越远离中心点处，其沉降越小，这符合路基沉降的一般规律。2种路基填料在颗粒流失之后，路基的地表位移均较流失之前增大，且路基填料的密度越大，沉降越小。

图12 颗粒流失前后碎石路基运行期地表沉降Fig.12 Ground settlement of gravel subgrade surface in operation period before and after particle loss

7 结 论

(1)不同密度试样在经历渗流作用后，均会出现不同程度的细颗粒流失现象。排除试样击实及加压固结过程中的颗粒破碎，颗粒流失主要发生在<2 mm的粒径区间。其中，密度为1.99 g/cm3的试样，在渗透破坏的时候颗粒流失率最高。

(2)选取85%(1.93 g/cm3)和88%(1.99 g/cm3)2个压实度开展渗透破坏前后的大型三轴压缩试验，其分别代表压实度较低和细颗粒流失率较高2种相对不利情况。据试验结果分析，碎石填料应力-应变关系呈应力强化型，填料工程特性较好。

2个试验压实度条件下，细颗粒流失前后的内摩擦角基本相同，分别约为39.5°和40.0°，但流失后试样的黏聚力有较明显的下降。密度1.93 g/cm3时，黏聚力从流失前的47.576 kPa降低到流失后的29.867 kPa；密度1.99 g/cm3时，黏聚力从流失前的44.056 kPa降低到流失后的28.789 8 kPa，表明细颗粒流失对碎石路基粗组构骨架的整体受力工程特性影响不大，不利影响主要体现细颗粒流失改变了试样内部粗细颗粒的粘结力。

(3)使用FORTRAN编程语言，结合改进欧拉方法，将亚塑性本构模型通过UMAT用户子程序嵌入ABAQUS有限元分析软件中。采用该模型进行三轴数值计算以及碎石路基颗粒流失前后的变形沉降计算。三轴数值计算结果表明，亚塑性本构模型可以很好地描述碎石填料在颗粒流失前后的非线性应力-应变关系。碎石路基变形沉降计算结果表明，颗粒流失之后，碎石路基的变形沉降会出现增大的趋势。且路基中心的地表沉降较大，越远离中心点处，其地表沉降越小。

(4)渗透作用导致碎石路基的颗粒级配变化，虽然对路基的强度无明显影响，但会对路基的变形产生一定影响，由此会产生路基的不均匀沉降，故保持路面的原颗粒级配是工程后期养护的关键。建议定期为已发生沉降的路基区段填补细粒土，或修筑防护措施保护细粒土不被流失，在设计和施工时，必须考虑将地下水和地面水排出和拦截在路基用地范围之外，并防止地面水漫流和下渗。