循环荷载下加筋土黏弹塑力学模型与数值分析

匡希龙 ，周志刚，李雨舟

(1. 长沙理工大学 公路工程教育部重点实验室，湖南 长沙，410014；2. 长沙学院 土木系，湖南 长沙，410003)

现代加筋土路堤技术是20世纪60年代发展起来的加固路基土体稳定性的一项新技术[1-3]。Sawicki[4]把加筋土体看成均匀的各向异性材料，假设其宏观应力由土的微观应力和筋材的微观应力组成，筋材与土间不产生相对滑动，建立了加筋土复合刚塑性模型。所谓宏观应力是指在物质宏观体积间相互作用的力，而微观应力则指作用在物体微观组织间的应力。Uzan等[5]在 3个参数黏弹性模型基础上，继续拓展一个新的黏塑性元件，并认为其黏滞系数是应变的函数，得到了加筋土黏弹性模型。本文作者在前人研究的基础上将Bingham黏塑性体改进为黏弹塑性体，试图构建一种更能反映加筋土复合体力学响应的黏弹性塑性计算模型，并通过数值分析来揭示其合理性。

1 基本假定

本文认为加筋土体同样满足Mohr-Coulomb准则，通过几个基本假定，分析宏观加筋土单元中土和筋材的微观应力，建立并分析基于黏弹塑特性的加筋土本构模型。

假定加筋土体是宏观均匀的各向异性复合材料，将加筋体中的土当作弹塑性材料，将筋材视为黏弹性材料，加筋土体的宏观应力保持不变，如图1所示。

图1 复合材料宏观与微观应力关系示意图Fig. 1 Relationship of macro-stress and micro-stress of composite materials

根据 Sawicki的研究，加筋土复合体的宏观应力及应变同土和筋材中的微观应力及应变之间存在下列关系

式中：n为加筋土中的体积含筋率，即筋材在加筋土中所占的体积比；上标s和r分别表示土和筋材。

(2) 假定土与筋材之间处于完全黏结状态而没有相对滑动，两者变形相协调。

当筋材与土处于完全黏结，同时不考虑筋材在厚度方向的压缩和弯曲变形时，加筋土复合体、筋材和周围土的宏观应变及微观应变为：

式中：e和Δh分别为筋材的厚度和筋材层间距。

2 加筋土黏弹塑性演化方程

2.1 加筋土初期的黏弹性效应

加筋土体在施工完成并通车的开始阶段处于黏弹性阶段。基于加筋土的受力特性按三参数黏弹性模型来表述，总应力及总应变为：

通过Laplace变换，应力应变关系式为：

式中：E1和E2为弹性模量；η为黏性系数。逆变换后，得

再令 0＜t＜T，其中 T为加筋土体塑性到达的时间，对式(9)进行Laplace变换可得：

式(10)即为加筋土体开始阶段的黏弹性计算公式。

2.2 加筋土黏弹性塑性力学特征转变响应

随着外部荷载作用时间的进一步增大至T之后，加筋土体开始向黏弹性塑性阶段转变。文献[5]只单独考虑筋材的黏弹塑特性，在3个参数黏弹性模型基础上拓展一个新的黏塑性元件，并认为其黏滞系数是应变的函数。假设筋材瞬时弹性应变为εe、黏弹性应变为εev以及黏塑性为应变εvp，其总应变表示为：

黏塑性应变速率可表达成：

然而，关于筋材与土共同构成的加筋土复合体的黏弹塑性本构模型，国内外的相关研究尚未达成共识[6-7]。针对加筋土的力学特点，本文效仿国内大量有关加筋沥青黏弹塑特性的研究成果[8]，并综合文献[5]再继续创建一个弹性元件和一个塑性元件，如图2所示，用6个参数来全面研究加筋土复合体的黏弹塑特性。采取的研究方法是基于 St.Venant流变模型(即由Hooke弹性体、Kelvin黏弹性体及Bingham黏塑性体串联而成)，并在此基础上将Bingham黏塑性体改进为黏弹塑性体，便得到加筋土黏弹塑性流变模型，对应的总应变直观地表示为：

图2 黏弹塑性流变模型Fig. 2 Model of viscoelastoplastic

首先应求出土体到达塑性状态所需要的时间 T，根据Mohr-Coulomb准则，按下式确定：

当考虑σz和sxσ为主应力时，

当t≥T时即加筋土体进入塑性阶段后，其应力表现为恒力，其变形体现出蠕变特性。

则加筋土复合体黏弹塑本构模型为：

为简化计算，假设以进入塑性阶段开始计时，即当t=0时，应力σ=σ0＞σs，同时保持σ不变，则有：

同时可得到：

将式(20)，(21)和(22)求和得到处于黏弹塑状态下的加筋土简化后的力学模型为：

上述关于加筋土复合体黏弹塑性本构关系的模型理论具有表述方便、直观等优点，全面反映了及通车后的加筋土路基在车辆的反复碾压下所具有的黏弹塑性材料全部特征，因此，可利用上述黏弹塑性力学理论来研究循环荷载作用过程中加筋土的力学响应。

3 高填方路堤加筋土黏弹性塑性有限元理论分析

随着计算机技术的飞速发展，用来进行数值分析的有限元软件也越来越多,应用领域也越来越广泛[9-10]。本文欲借助有限元通用软件ANSYS10.0进行数值模拟黏弹塑性高填方路堤加筋土的长期工作性能，并与土工离心模型试验进行对比分析，更好地验证黏弹塑性模型的合理性。

3.1 加筋土黏弹塑有限元计算模型

本文基于黏弹塑性有限元理论，选用 Drucker-Prager弹塑性模型来模拟填土及Kelvin黏弹性模型来模拟筋材，结合上述St.Venant流变模型，并在此基础上将Bingham黏塑性体改进为黏弹塑性体，将增量法与初应变迭代法相结合，模拟逐层填筑过程和分级加载非线性过程，以达到ANSYS 2次开发的目的。

在计算过程中将弹性-黏弹性-黏塑性一体化单元加入ANSYS单元库，利用ANSYS建模、网格划分，对结构刚度、质量矩阵进行非线性的迭代计算。在矩阵数据运算时计入随时间的黏性，再将这些计算结果返回ANSYS处理，实现结构的图形显示和数据提取。

在复杂应力状态下，总应变速率分为弹性应变速率和黏塑性应变速率：

式中：{ε}为总应变速率；{εe}为弹性应变速率；{εvp}为黏塑性应变速率。

应力速率{σ}可以表示为：

式中：[D]为弹性矩阵。

开始产生黏塑性变形的屈服条件可以表示为：

对于黏弹塑性体的应力计算，在时间段Δtn=tn+1-tn内产生的应力增量为：

式中：[B]为单元的应变矩阵。

数值模拟分析流程步骤如图3所示。

图3 有限元分析计算流程Fig. 3 Calculation process of numerical analysis

根据填筑过程与筋材布置方式建立高填方路堤加筋土有限元计算模型，土体采用四节点等参单元模拟，土工网采用杆单元模拟，在填土与筋材之间设置界面单元。节点总数为3 513，单元总数为1 906，其中土体单元数为1 128，筋材杆单元数为225，界面单元数为 553。加筋土底部为铰支座，路堤中央土体的边界采用水平约束条件，路堤外侧和土侧的边界为自由边界条件。

在非线性有限元计算中，将整个填筑与加载过程分成2个阶段，即填筑完成阶段和分级施加均布荷载阶段。对于上述各个指定阶段的变形与应力状态进行分析，按照与试验中相同的处理方式，在顶部沉降与边坡侧向变形中扣除了填筑完成时的计算变形值。

3.2 试验与理论对比分析

为了研究加筋土之间的相互作用并合理地解释加筋土复合机理，国内外学者进行了大量的加筋土实验研究。离心模型试验因其独特的优势已成为验证计算方法和解决土工问题的一种强有力的手段[11-12]。本文在长沙理工大学土工离心模型试验中心进行了加筋土离心模型试验，图4所示为土工离心模型试验。

所用的离心机型号为 TLJ-150A，最大容量 150 g/t，旋转半径3.5 m，最大加速度150 g；选用的大模型箱尺寸：长0.9 m、宽0.7 m、高0.7 m，模型比尺n=30，结构原型如图5所示。

图4 土工离心模型试验Fig. 4 Beotechnical centrifuge modeling test

图5 结构原型Fig. 5 Structure prototype

实验选用的研究对象为红砂岩填料，其物理参数为：wL=31%，wP=17%，Ip=14，γ=18.91 kN/m3，w=4.7%，干容重为 γd=21.7 kN/m3，泊松比 μs＝0.35，c=17.1 kPa，φ=22°。

本文旨在借助上述有限元理论，以离心模型试验获得的实测值为衡量标准，对加筋土黏弹性模型和黏弹塑模型进行定性的对比，从而揭示黏弹塑模型理论值的变化趋势更接近实测值。为此，通过基于Levenberg-Marqurdt和全局优化问题的最优化求解理论来求解其模型参数：E1=165kPa，E2=397kPa，E3=566 kPa，η2=4 954 kPa·d，η3=1 325 kPa·d。第 2，3 和 4 层加筋土体的屈服强度分别为σs2=77.32 kPa，σs3=51.54 kPa，σs4=25.80 kPa。

经数值分析得到加筋土理论计算与实测结果的对比曲线，如图6～7所示，并将理论值与实测值求解误差归纳成表1～2。

表1 2种模型针对加筋土体侧向变形求解误差Table 1 Calculation error of lateral deformation on two models

表2 2种模型针对加筋土体竖向变形求解误差Table 2 Calculation error of vertical deformation on two models

从表1～2及图6～7可知：同实测值相比，黏弹塑性模型与黏弹性模型所得到的理论值均不同程度的存在误差，但总体来看，黏弹塑性模型的理论值对比实测值，其拟合误差可减少 8%以上，在一致性方面也明显比黏弹性模型要好，由此说明同以往的黏弹性模型相比，本文构建的黏弹塑性模型更能反映加筋土复合体的力学特性和变形特点。

图6 各观测点侧向变形的理论值与实测值对比Fig. 6 Contrast of theoretical value and actual value of lateral deformation

图7 各观测点竖向变形的理论值与实测值对比Fig. 7 Contrast of theoretical value and actual value of vertical deformation

4 结论

视加筋土为土与筋材组成的复合材料，从而建立了黏弹塑特性的加筋土体本构模型。通过ANSYS10.0数值分析和离心模型试验相互验证，充分表明拓展成的黏弹塑性力学响应模型，更能反映循环荷载作用下加筋土复合体的变形特征。

(1) 加筋土体在施工完成并通车的开始阶段处于黏弹性阶段。基于加筋土的受力特性按3个参数(E1，E2，η)黏弹性模型来表述初期的应力-应变关系式

(2) 当作用循环时间 t到达土体开始进入塑性状态所需要的时间T后，按6个参数法(E1，E2，E3，η2，η3，σs)构建加筋土黏弹塑本构模型，即：

(3) 将增量法与初应变迭代法相结合，模拟逐层填筑过程和分级加载非线性过程，运用ANSYS对黏弹塑性有限元理论进行二次开发和数值模拟分析,充分证实了黏弹塑模型的合理性。

[1] 周志刚, 郑健龙, 李强. 土工加筋材料处治填挖交界路基非均匀沉降设计方法研究[J]. 中国公路学报, 2003, 16(1): 27-31.ZHOU Zhigang, ZHENG Jianlong, LI Qiang. Study of the design method of using geonets to treat the unequal settlement of subgrade at the cut to fill location[J]. China Journal of Highway and Transport, 2003, 16(1): 27-31.

[2] 李丽华, 王钊, 陈轮. 考虑筋材蠕变特性的加筋土流变模型[J]. 岩土力学, 2007(8): 34-37.LI Lihua, WANG Zhao, CHEN Lun. Rheological model considering creep of geosynthetics for reinforced soil[J]. Rock and Soil Mechanics, 2007(8): 34-37.

[3] 介玉新, 李广信. 加筋土数值计算的等效附加应力法[J]. 岩土工程学报, 1999, 21(5): 614-616.JIE Yuxin, LI Guangxin. Equivalent additional stress method for numerical analysis of reinforced soil[J]. Chinese Journal of Geotechnical Engineering, 1999, 21(5): 614-616.

[4] Sawicki A. Rheological model of geosynthetic- reinforced soil[J].Geotextiles and Geomembranes, 1999, 17(1): 33-49.

[5] Uzan J A, Sides Perl. Viscoelastoplastic model for predicting performance of asphalt mixture[J]. Transportation Research Record, 1985: 78-79.

[6] Radhey Sharma, CHEN Qiming. Analytical modeling of geogrid reinforced soil foundation[J]. Geotextiles and Geomembranes,2009: 63-72.

[7] 栾茂田, 肖成志, 杨庆, 等. 土工格栅蠕变特性的试验研究及黏弹性本构模型[J]. 岩土力学, 2005, 26(2): 187-192.LUAN Maotian, XIAO Chengzhi, YANG Qing, et al.Experimental study on creep properties and viscoelasticity constitutive relationship for geogrids[J]. Rock and Soil Mechanics, 2005, 26(2): 187-192.

[8] 陈骁, 应荣华, 郑健龙, 等. 基于 MTS压缩试验的热态沥青混合料黏弹塑性模型[J]. 中国公路学报, 2007, 20(6): 25-30.CHEN Xiao, YING Ronghua, ZHENG Jianlong, et al.Viscoelastic-plasticity model of hot asphalt mixtures based on MTS compaction test[J]. China Journal of Highway and Transport, 2007, 20(6): 25-30.

[9] 刘汉龙, 吴维军, 高玉峰. 土工织物加固堤防非线性有限元分析[J]. 岩土力学, 2003, 24(1): 79-82.LU Hanlong, WU Weijun, GAO Yufeng. Nonlinear finite element analysis of reinforced dike with geotextile[J]. Rock and Soil Mechanics, 2003, 24(1): 79-87.

[10] 钱劲松, 凌建明, 黄琴龙. 土工格栅加筋路堤的三维有限元分析[J]. 同济大学学报: 自然科学版, 2003, 31(12):1422-1425.QIAN Jinsong, LING Jianming, HUANG Qinlong. 3D finite element analysis of embankment reinforced by geogrid[J].Journal of Tongji University: Natural Science, 2003, 31(12):1422-1425.

[11] 胡红蕊, 陈胜立, 沈珠江. 防波堤土工织物加筋地基离心模型试验及数值模拟[J]. 岩土力学, 2003, 24(3): 389-394.HU Hongrui, CHEN Shengli, SHEN Zhujiang. Centrifugal &numerical modeling of geosynthetic reinforced embankments on soft clay[J]. Rock and Soil Mechanics, 2003, 24(3): 389-394.

[12] 胡黎明, 劳敏慈, 濮家骝, 等. LNPALs在非饱和土中迁移的离心试验模拟[J]. 岩土工程学报, 2002, 24(6): 690-694.HU Liming, LAO Minci, PU Jialiu, et al. Centrifuge modeling of LNAPLs migration in unsaturated soil[J]. Chinese Journal of Geotechnical Engineering, 2002, 24(6): 690-694.