路基压实粉质黏土动态回弹模量的四参数预估模型及有限元实现

董 城, 刘文劼, 李 亮

(1.湖南省交通科学研究院有限公司,长沙 410015;2.中南大学 土木工程学院,长沙 410075)

长期以来，我国在路基路面结构分析中采用是基于静力学的弹性层状体系理论，假定各层结构为匀质材料，每层材料只由回弹模量和泊松比描述[1-2]。该理论在低速、低轴载条件下可以被接受，但为了适应高速、重载公路的迅猛发展，开展基于动力学理论的路基路面结构设计方法成为了必然。在动力学框架下，结构内部的应力状态时刻发生变化，因此，采用受应力状态影响的动态回弹模量作为路基力学特性表征参数将能更好地与动力学分析相结合，更准确地反映路基实际变形情况[3]。

国内外许多学者和研究机构在深入分析应力状态对动态回弹模量影响规律的基础上，陆续建立了一些动态回弹模量预估模型以尝试定量描述路基土的变形行为，其中，考虑剪切影响的模型包括双线性模型、幂指数模型，考虑侧限影响的模型包括围压模型、K-θ模型，但这两类模型的主要缺点在于仅考虑了体应力或偏应力对路基土动态回弹模量的影响，而实际上大部分土与粒料的动态回弹模量是体应力、偏应力两者的函数，因此综合考虑体应力和偏应力的影响的复合类模型可以更真实、全面地反映材料的力学性状，Uzan等[4-5]都相继提出了复合类预估模型，美国于2004年在NCHRP 1-37A项目中提出公路路基设计的动态回弹模量模型(简称N37A模型)[6]，消除了模量不定值和量纲不一致问题，目前被工程界广泛接受。后续许多回弹模量预估模型均是在N37A回弹模型的基础上进行若干修正[7-8]。

针对具体路基土，建立合适的动态回弹模量预估模型并将其移植入有限元软件中，是利用动态回弹模量对路基路面进行结构有限元动力学分析的前提，但由于我国路面结构、所处的自然环境以及路基土材料属性的差异，一些国外成熟的路基土回弹模量预估模型不能直接运用于我国路基路面结构设计，仍需结合我国的地质条件和实验数据进行改进。本研究以湖南湘东地区路基粉质黏土为例，通过室内动三轴试验测定不同应力状态、压实度、含水率下粉质黏土的动态回弹模量，分析影响粉质黏土动态回弹模量的因素及相关规律，并在N37A模型的基础上提出含四参数的N37AP1模型，通过多元回归分析，得到了不同压实度下粉质黏土对应的N37AP1模型参数。在此基础上，严格推导出N37AP1模型的一致切线刚度矩阵，编写用户自定义材料子程序(UMAT)将其移植到有限元软件ABAQUS 中，通过单个单元的轴压、围压加载以及典型路基路面结构分析，验证所推导一致性切线刚度矩阵及UMAT编写的正确性，为粉质黏土路基—路面结构设计提供更真实有效的数值模拟方法，也可供其它种类路基填土的类似研究借鉴参考。

1 动态回弹模量试验及数据分析

1.1 试验概况

本文选取湖南湘东地区具有代表性的粉质黏土路基填料作为研究对象，其基本物理性质指标列于表1。

表1 土样基本物理性质指标

试验参照我国《公路路基设计规范》(JTG D30—2015)[9]附录A进行，将粉质黏土在不同含水率(wopt+3%、wopt、wopt-3%)和不同压实度K(93%、96%)条件下制备试样。路基土应力加载序列采用国内已有研究成果，动偏应力范围为30～105 kPa，围压范围为15～60 kPa，动态回弹模量测试中，某一加载级位下荷载循环次数为100 次，记录各级重复荷载作用下最后5 次循环的回弹变形平均值作为动态回弹模量的计算依据。

1.2 应力状态与回弹模量的关系

图1为w=wopt状态下粉质黏土动态回弹模量随应力状态的变化曲线。根据图1(a)和图1(b)，粉质黏土动态回弹模量随着动偏应力的升高而减小，且减小的速率逐渐增大；而根据图1(c)和图1(d)，动态回弹模量将随着体应力的升高而增大，基本呈线性变化或速率逐渐放缓。总体上看，压实度对动态回弹模量的影响也比较明显，如当体应力θ=210 kPa且σd= 30 kPa时，相对93%压实度，压实度为96%时动态回弹模量提高了16.65%。

(a)

(b)

(c)

(d)

2 动态回弹模量模型改进及评价

2.1 已有模型的改进及回归分析

在建立动态回弹模量的预估模型时，首先应尽可能准确地按本构定律建立起应力-应变的关系模型，而其他因素的影响则可通过模型参数予以反映。已有研究和本文试验均表明，路基土的动态回弹模同时受体应力和偏应力的双重影响，考虑该现象的复合式模型中具有代表性的N37A模型形式为

(1)

对于不同路基土，剪应力对路基土回弹模量的影响程度也是不同的，据此考虑在N37A模型形式基础上，引入k4项用于调节偏应力对回弹模量的贡献比例，得到的改进模型(简称N37AP1)具体为

(2)

当k4取很小值时，该模型逼近于K-θ模型[10]；当k4取1则成为N37A模型。

分别采用N37AP1模型、N37A模型和UZANP1[11]模型对对本次试验结果进行了回归分析，得到了粉质黏土在不同条件下的预估模型参数和对应的决定系数R[12-13]，回归结果见表2。

表2 预估模型参数回归结果(k4=0.1)

由表2可见，在不同的含水率和压实度下，UZANP1模型的决定系数最低，N37A模型次之，而N37AP1模型的决定系数在同条件下都大于其它两者，对于试验结果的回归效果明显更好。

2.2 不同模型预估结果对比

为了验证N37AP1模型的预估效果，测试了更多应力状态下粉质黏土的动态回弹模量，表3给出了w=wopt，K=96%条件下粉质黏土动态回弹模量实测数据和三种模型预估值，总体上改进模型(k4=0.1)的每项误差值都最小，而且在剪应力较小时误差都在5%以下，当剪应力增大时误差略有增加，最大为8.3%。

综上分析可知，N37AP1模型对最佳含水率下的粉质黏土的动态回弹模量的预估能力明显更加突出，误差大体上是N37A模型的50%～60%，是UZANP1模型的20%～45%。N37AP1模型通过调整k4取值的大小，改变偏应力对粉质黏土动态回弹模量的贡献，进而改变了动态回弹模量的预估结果，使得预估模型能更好的拟合试验数据。

3 N37AP1一致性切线刚度矩阵的推导

建立动态回弹模量预估模型的主要目的之一为了将其有效应用于有限元计算，本文利用ABAQUS的UMAT接口将N37AP1模型移植入该软件中。ABAQUS /standard分析模块的每一步非线性求解都是通过采用切线刚度矩阵进行的，因此要实现N37AP1模型的二次开发就必须求得其全应力状态下的一致切线刚度矩阵，求解过程如下：

线弹性本构关系可表示为

(3)

式中：E为应变张量；S为应力张量；E为线弹性模量；ε=tr(E)为体应变；v为泊松比；α=ν/(1-2ν)。

类似于线弹性本构关系，非线性弹性本构关系可以利用刚度矩阵Mr写成

表3 动态回弹模量预测与试验结果对比(k4=0.1)

(4)

令

(5)

式中：k=k1/(1+ν)；k1,k2,k3,k4为N37AP1模型参数；θ为体应力；Pa为标准大气压；τoct为八面体剪应力。

结合式(4)和式(5)，动态回弹模量本构关系可以简写为

S=C(θ,τoct)(αεI+E)

(6)

由式(6)可得体应力表达式

(7)

体应力和偏应力可以简写为

(8)

(9)

将C(θ,τoct)的表达式代入式(9)，则式(9)可以写成

(10)

显然θ,τoct与ζ,γ之间存在函数关系

(11)

因此，C=C(θ,τoct)可以用应变表示为：C=C(ζ,γ)，进而，式(6)可以表示成

S=C(ζ,γ)(αεI+E)

(12)

对式(12)求导可以得到材料的切线刚度矩阵D=∂S/∂E

(1+αI⊗I)+(αεI+E)⊗▽EC

(13)

由式(12)和式(13)可得

(14)

由定义υ=|ε|，可得

(15)

(16)

(17)

由式(5)可得

(18)

由式(10)可得

(19)

(20)

结合式(9)将切线刚度矩阵写成偏应力的形式

(21)

为方便编程，将切线刚度矩阵写成如下形式

(22)

按照以上推导结果，在UMAT中以当前应力张量计算出一致切线刚度矩阵(雅可比矩阵)[14]，然后结合当前应变增量得到应力增量，最后根据当前应力和应力增量更新当前的应力，将更新后的应力传入下一步计算，即可实现预估模型在ABAQUS软件中的二次开发(见图2)。

图2 UMAT实现流程图

4 预估模型有限元实现验证

4.1 预估模型的单元验证

为了验证所推导一致性切线刚度矩阵及UMAT编写的正确性，在ABAQUS中先采用单个单元进行验证，验证过程中k4取不同值。采用材料参数见表4。

表4 N37AP1模型的材料参数

有限元模型采用单个单元C3D8，约束相邻三个面的法向位移，在其余三个自由面上施加法向荷载。轴力和侧压力均为渐变形式加载，初始荷载：σ1=0 kPa，σ2=σ3=0 kPa；最终荷载：轴力σ1=100 kPa，σ2=σ3=10 kPa,如图3所示。

图3 单元加载示意图

应变的解析解可表达为

(23)

提取单个单元回弹模量随轴力变化曲线以及轴向应力应变和侧向应力应变曲线，如图4～图6所示。

图4 Mr与σ1的关系曲线

图5 ε1与σ1关系曲线

图6 ε3与σ3关系曲线

由图4～图6可知，在渐变加载工况下，三组参数所得的应力应变曲线和理论值十分接近，因而推导的一致切向刚度矩阵和编写的UMAT是正确的。同时，也可以看出k4在改进模型中所起的作用，能调整八面体剪切应力在预估模型中的比重，在相同应力状态下，k4值越大回弹模量预估值越小，预估应变越大，路基土剪切软化效果越明显；反之则回弹模量预估值越大，预估应变越小。

4.2 预估模型的结构验证

为进一步验证预估模型，将其移植入了路基—路面有限元模型中，该模型尺寸参照高速公路粉质黏土路基—沥青路面结构，取 1/4模型进行研究，如图7所示，模型参数参考湖南省某高速公路的设计参数，如表5所示，其中路基土的材料参数采用w=wopt时粉质黏土的动态回弹模量改进模型参数。车辆荷载形式为单轴双轮组，轴载大小为100 kN，为了便于设置轮胎接地区域，将该区域等效为与网格形状一致的矩形[15]，考虑行车道和超车道同时存在车辆的情况。

表5 路基路面结构参数

图8为在自重作用下路基路面结构路基土回弹模量的分布云图，由于路基土采用了应力相关的动态回弹模量，路基内不同位置的回弹模量值各不相同。

图7 整体网格

图8 自重作用下路基土回弹模量云图

图9为路基各点回弹模量在竖向和横向上的分布曲线。从图9(a)可以看出，在自重作用下路基土回弹模量随深度增加而有所增加，最终趋于稳定，路床的回弹模量值在70～100 MPa，路基本体的回弹模量值在60～85 MPa，由于K=0.96和K=0.93路基土的模型参数不同，回弹模量在路床底有了突变，体现了上下层压实度差异造成的影响。从图9(b)可以看出，路基的回弹模量由路基中央向两侧衰减，但深度越大，衰减效果越不明显。

图10为施加车辆荷载后，由于应力状态变化形成的回弹模量相对变化率云图。对于行车道和超车道下方的路基，由于偏应力的明显增加，回弹模量出现了衰减，最大衰减幅为-1.3%；路床层的车轮间隙区域出现回弹模量的增加，最大增幅约2.0%。由此可以看出，采用动态回弹模量预估模型有效体现了路基各点回弹模量因车辆荷载作用产生的时间变化。

(a)竖向分布

(b)横向分布

图10 施加车辆荷载后路基土回弹模量的相对变化率云图

5 结 论

(1)结合粉质黏土动态回弹模量试验数据，对已有的动态回弹模量预估模型进行了改进，提出综合N37A模型、UZANP1模型和K-θ模型的四参数改进模型N37AP1，确定了N37AP1模型对不同含水率和压实度的粉质黏土的预估参数，通过与实测值的比较，验证了改进模型对粉质黏土动态回弹模量的预估效果更为理想。

(2)基于广义虎克定律推导了N37AP1模型的精确一致切线刚度矩阵，通过用户子程序UMAT成功将该模型嵌入到有限元软件ABAQUS。

(3)单个单元测试结果表明，k4在改进模型中能调整八面体剪切应力在预估模型中的比重，k4=0.1、k4=1、k4=10时所得的应力应变曲线都和理论值十分接近，在相同应力状态下，k4值越大回弹模量预估值越小，预估应变越大，路基土剪切软化效果越明显。

(4)对典型粉质黏土路基—沥青混凝土路面结构的分析结果表明，对路基土采用N37AP1模型可以有效反映路基各点在自重作用下回弹模量的空间差异，且可以体现压实度的影响；同时能反映车辆荷载作用下回弹模量的演变，实现了模量与应力状态的动态耦合。四参数预估模型在有限元软件中的嵌入及运用为路基路面结构设计提供了更为真实有效的数值模拟方法。

参 考 文 献

[1] 罗志刚. 路基与粒料层动态模量参数研究[D].上海:同济大学, 2007.

[2] 邓学钧.路基路面工程[M].北京:人民交通出版社,2010.

[3] 曹海莹,朱毅,刘云飞,等. 车辆荷载作用下双层路基层间动应力响应试验研究[J]. 振动与冲击, 2017, 36(5): 6-10.

CAO Haiying,ZHU Yi,LIU Yunfei,et al. Testing research on dynamic stress response at interlayer of two-layer roadbed under vehicle load[J]. Journal of Vibration and Shock, 2017, 36(5): 6-10.

[4] LEKARP F, ISACSSON U, DAWSON A. State of the art I: resilient response of unbound aggregates[J]. Journal of Transportation Engineering,2000,126(1):66-75.

[5] LI D, SELIG E T. Resilient modulus for fine-grained subgrade soil[J]. Journal of Geotechnical Engieering, ASCE, 1994,120(6):939-957.

[6] 董城,冷伍明,李志勇.粉土动态回弹模量试验研究[J].中南大学学报(自然科学版),2012,43(12):4834-4839.

DONG Cheng, LENG Wuming, LI Zhiyong. Dynamic resilient modulus of silt[J]. Journal of Central South University(Science and Technology), 2012,43(12):4834-4839.

[7] 陈声凯,凌建明,罗志刚. 路基土回弹模量应力依赖性分析及预估模型[J].土木工程学报, 2007, 40(6): 95-99.

CHEN Shengkai, LING Jianming, LUO Zhigang.Stress-dependent characteristics and prediction model of the resilient modulus of subgrade soils[J]. China Civil Engineering Journal, 2007, 40(6): 95-99.

[8] 凌建明,苏华才,谢华昌,等. 路基土动态回弹模量的试验研究[J]. 地下空间与工程学报,2010,6(5): 919-925.

LING Jianming, SU Huacai, XIE Huachang, et al. Library research on dynamic resilient modulus of subgrade soil[J]. Chinese Journal of Underground Space and Engineering, 2010,6(5): 919-925.

[9] 中华人民共和国交通运输部公路路基设计规范:JTG D 30—2015[S].北京:人民交通出版社,2015.

[10] WITCZAK M W,UZAN J. The universal airport pavement design system[R]. Report I of V: Granular Material Characterization, Department of Civil Engineering, University of Maryland, College Park, 1988.

[11] 董城.路基土动态回弹模量与临界动应力及路基路面性能研究[D].长沙:中南大学,2014.

[12] TOMMY C, HOPKINS, TONY L, et al. Resilient modulus of Kentucky solis[R]. Lexington: University of Kentucky, 2004.

[13] 刘明. 线性回归模型的统计检验关系辨析[J]. 统计与信息论坛, 2011, 26(4): 21-24.

LIU Ming. Analysis of the relationship among the statistical tests on linear regression model[J]. Statistics & Information Forum,2011, 26(4): 21-24.

[14] 费康,张建伟.ABAQUS在岩土工程中的应用[M].北京:中国水利水电出版社,2010.

[15] 王仕传.路床加筋力学行为与设计方法的研究[D].上海:同济大学,2007.