近年来，利用土工加筋技术处理软土路堤的方法被广泛应用到实际工程中，但其作用机理和加筋路堤的整体性能还有待进一步探究，其中移动荷载下土工加筋路堤的动力响应是人们非常关心的问题，很多学者对此进行了研究．在解析研究方面，Maheshwari等[1]提出将土工加筋路堤看作双梁模型，并推导了条形荷载作用下的解析解；在此双梁模型的基础上，Zhao等[2]将土工加筋层看作Timoshenko梁，考虑了梁所受的剪切作用．Zhang等[3-4]建立了考虑地基土阻尼系数的双欧拉梁模型，并推导了在移动荷载作用下梁的解析解．在数值模拟研究方面，郑俊杰等[5-6]建立了三维离散元筋土模型，系统研究了拉拔阻力、格栅应变、局部孔隙率等力学响应并揭示了其发展规律；张玲等[7]建立了土工加筋路堤的三维数值模型，分析了其在移动荷载下的应力和变形规律．上述成果对于研究交通荷载作用下土工加筋路堤的动力响应特征具有十分积极的意义．但是，这些研究均将地基土体视为均质，不能反映出岩土材料自身所具有的随机性与不确定性的特征．Lumb等[8]首次提出了岩土参数空间变异性的概念，指出由于各种不确定的外力等因素导致土性参数随着空间分布呈现出变异特征．Vanmarcke等[9]提出土性参数随机场的概念及模型，用于描述土体参数的空间变异性．在随机场模型的基础上，国内外学者针对边坡稳定性[10-11]、地基承载力[12]、隧道变形特征[13-14]、水泥稳定土变异性[15-16]等问题开展了广泛的研究．其中大量研究是针对土体弹性模量变异性展开的，结果显示空间变异性可能会对工程结构的安全程度及力学响应造成影响，因此进行结构响应及性能研究时考虑土性参数的随机性是十分必要的，尤其是弹性模量及强度参数等的变异性．在道路工程中的随机场研究方面，LUA等[17]建立了机场路面的二维有限元模型，首次将随机有限元分析用于计算路面响应；Caro等[18]研究了二维柔性路面中内部孔隙分布的随机性对路面响应的影响，通过建立弹性模量与孔隙率之间的关系，将孔隙分布的随机性转换到了弹性模量的随机性上；Castillo等[19]研究了柔性路面沥青层的弹性模量空间变异性对力学响应的影响．上述研究集中在对路面面层材料变异性的实现，缺少对地基土体材料弹性模量变异性的研究，但天然土体形成时具有沉积特性，其变异性不可忽视．此外，未有研究在加筋路堤动力响应中考虑空间变异性，变异性是否对加筋路堤也有显著影响须进一步研究．受计算效率等限制，很少有基于三维路面模型的大量随机动力响应研究，此前研究中数值计算组数只有几十组[20]，无法反映随机性的一般规律，须提升计算数量级并进行大量模拟计算．基于以上的分析，本研究建立了土工加筋路堤的三维有限元模型，结合随机场理论和蒙特卡罗模拟，分析了地基土弹性模量(E)的变异性对移动荷载作用下加筋路堤变形特征的影响，分别探讨了随机场的变异系数和自相关距离对路堤横纵截面挠度的影响规律，并利用概率分析方法总结了挠度的分布特征．1 随机场理论岩土工程中的土体性质具有显著的不均匀性和随机性．即使同一类土体，其性质随着空间位置的不同往往也不同，即具有空间变异性．岩土参数的空间变异性不仅体现为单一的随机性，还具有很强的相关性，距离越近相关性越大，距离越远相关性越小．随机场理论中常用自相关函数来描述这一相关性，选用三维高斯型自相关函数，计算公式为      ρ(a,b)=expxa-xblx2+ya-ybly2+za-zblz2, (1)式中：a和b分别为空间中的两点；(xa，ya，za)和(xb，yb，zb)分别为a和b点的坐标；lx，ly和lz分别为x，y和z方向的自相关距离，自相关距离越大，参数的空间相关性越强[10]．选用连续性好且展开效率高的K-L级数展开法实现随机场的离散．岩土参数通常可视为符合对数正态分布的均质随机场，对数正态随机场的离散结果可表示为ω(x,θ)=expμln+σln∑i=1Mλifi(θ)ξi(θ)，(2)式中：μln和σln分别为对数正态随机场的均值和标准差；M为展开项数；ξiθ为相互独立的标准正态随机变量；λi和fiθ分别为自相关函数的第i个特征值和特征向量，λi和fiθ通过求解如下的Fredholm二类积分方程得到∫Ωρ(x,x')fi(x')dx'=λifi(x)，(3)其中，Ω为积分区域，x和x'为随机场中的任意两个坐标点[19]．对于离散的随机场，理论上式(2)中的展开项数M应为无穷大，但在保证计算精度的前提下，通常只截取前n项以保证计算效率，采用场期望能比率因子ε来判别K-L展开式的n[10]，表达式为ε=∑i=1nλi/∑i=1∞λi，(4)n应尽可能地保证比率因子ε接近于1．首先在ABAQUS中建立确定性分析模型，根据该模型生成相应的土体参数随机场；然后通过Python编程二次开发实现每一个土体单元参数的重新编写，并生成新模型；最后批量提交计算．2 数值模型2.1　几何尺寸及参数在文献[7]的数值模型基础上进行空间变异性分析，从上至下分别为路堤顶层(高度为h1)、路堤填土(高度为h2)、加筋垫层(高度为h3，由碎石垫层和土工加筋合成材料组成)和地基土(高度为h4)，路堤顶面宽度为L1，路堤底面宽度(即加筋垫层宽度)为L2，地基土体宽度为L3，路堤边坡坡率为m，路堤纵向长度为L．基于ABAQUS有限元软件，建立加筋土路堤的3D数值模型，模型的几何参数如下：h1=0.3 m，h2=3.0 m，h3=0.4 m，h4=6.0 m，L1=10 m，L2=19 m，L3=35 m，L=10 m，m=1.5．路堤各结构层参数取自文献[7]，土工合成材料采用线弹性模型，其他材料采用Mohr-Coulomb本构模型，各材料参数见表1．10.13245/j.hust.238925.T001表1材料参数取值材料密度/(kg∙m-3)弹性模量/MPa泊松比内摩擦角/(°)黏聚力/kPa路堤顶层2 5003 0000.20路堤填土1 850150.302818碎石垫层2 100300.25300地基土体1 800100.302520土工格栅1 0000.002.2　初始条件与边界条件以路堤填筑完成后的地应力平衡状态为应力初始状态，在ABAQUS中通过设置地应力分析步及自带的重力加载进行地应力平衡．地应力平衡分两次完成，路堤填筑前进行一次平衡，填筑完成后再进行一次平衡．加筋路堤的分层填筑仅考虑重力的叠加作用，暂不考虑排水固结．移动荷载下路堤的边界条件为：在地基土体的底部设置为固定约束，约束横向(U1)、竖向(U2)和纵向(U3)的位移；约束沿路堤纵向的两个横截面水平方向的位移，即固定U3；约束路堤坡脚远端的地基土体边界水平方向的位移，即固定U1．2.3　单元类型与网格如图1所示，采用八节点线性六面体单元(C3D8R)对路面顶层、路堤填土、碎石垫层及地基土体进行网格划分[7]，车辆作用在路面顶层，这一区域的范围内应力应变较大，为确保计算精确度并保证计算效率，对路面顶层，尤其是在车道位置处的网格进行加密．采用二维线单元模拟土工加筋材料，使用T3D2桁架单元对其进行网格划分，并且组合形成网状加筋结构．10.13245/j.hust.238925.F001图1土工加筋路堤网格2.4　移动荷载移动荷载形式采用单轴双轮组标准轴载(BZZ-100)，轴载大小为100 kN．为便于模拟，将轮胎作用面简化为0.22 m×0.16 m的矩形[21]，荷载移动速度设置为72 km/h．采用Fortran子程序编写荷载作用范围、大小及移动速度．移动荷载如图2所示．10.13245/j.hust.238925.F002图2移动荷载示意图3 随机性分析3.1　随机有限元模型用本模型计算文献[3]中的上梁挠度和下梁挠度，即本模型中结构层上表面位移和垫层位移，并将二者进行对比，数值模型计算结果与Zhang等[3]的结果具有一致的趋势，且在数值上具有较高的可信度，因此可以保证本模型计算的准确性．由于所实现的是弹性模量随机场，弹性模量变异性对结构变形的影响较大[13-14]，因此只以路面的变形特征作为研究切入点．变异系数(ν)和土体的自相关距离(θ)是表征土体随机场的两个最主要参数，决定了生成的随机场的特征，因此分析部分围绕这两个参数展开．取土体弹性模量均值10 MPa，θ=1，5，10，50 m；ν=0.1，0.3和0.5进行计算．所选模型为三维模型，单次模型计算量较大，综合考虑计算效率与计算准确度，选取每组随机计算次数为200．采用Intel Xeon W-2255 CPU@ 3.7 GHz，64 GiB内存服务器进行计算，单次随机计算时间约1.5 h．对随机场进行离散，并对每一个土体单元的属性进行重新赋值后，得到了路堤的随机有限元模型，如图3所示．地基土被分成1.24×104个单元，不同的颜色表示单元具有不同的弹性模量值(不表示弹性模量值之间的差异大小)．一次随机实现的单元弹性模量频数分布直方图符合对数正态分布．10.13245/j.hust.238925.F003图3土工加筋路堤随机有限元模型3.2　变异系数的影响图4为θ=5 m，ν=0.1，0.3和0.5，当荷载作用在路面中心时，200次随机计算得到的路堤表面挠度(γ)曲线，l为距路面中心的距离；黑色曲线为确定性分析结果，灰色曲线为随机计算结果，蓝色和粉色虚线为其中的两次计算结果．可以看到：随机性曲线在确定性曲线附近波动，且曲线形状类似．随着ν不断增大，地基土的不均匀性增大，随机性曲线趋于离散，表现为随机性曲线组成的包络线范围在不断增大，经计算得到3种变异系数对变形的最大影响分别为3.36 %，8.44 %，17.32 %．10.13245/j.hust.238925.F004图4路堤表面挠度曲线(θ=5 m)当ν=0.1时，随机计算得到的曲线形态与确定性分析几乎一致；随着ν的增大，随机性曲线形态发生了明显变化：确定性分析曲线沿着路面中心为对称分布，而随机计算曲线中出现了不对称的现象．图5给出了确定性分析和一次随机计算的挠度云图，确定性结果沿路堤中心对称，而随机计算结果中出现了明显的偏移．此外，所有的随机计算结果并非一直处于确定性曲线的上方或下方，而是会与之出现交叉、重合的现象．10.13245/j.hust.238925.F005图5确定性分析和随机计算的挠度云图对比图4种不同ν值的曲线位置可以发现：随着变异系数的增大，随机计算曲线表现为多数出现在确定性曲线的下方，即计算所得挠度趋于偏大，这是由于较大的变异系数导致土体内产生了低弹性模量区域，称作低刚度区域，低刚度占优效应[13]导致路堤整体的变形增大．低刚度占优效应即当结构或系统中出现刚度较低的区域时(表现为低弹性模量区域)，变形或裂缝会趋于向低刚度区域发展，从整体上看会导致结构的整体承载能力下降，变形增大．图6给出了θ=5 m时，不同ν下的路堤最大挠度(γk)的累计概率(P)分布曲线．对于每一个给定的控制挠度，可以由曲线得到超过该控制挠度的概率．随着ν的增大，路堤顶面最大挠度的概率分布趋于分散，同时超过同一控制挠度的概率也逐渐增大，即路堤结构的安全程度随着变异性的增强而相应降低．图7为θ=5 m时，200次随机计算所得的路面中心点位平均挠度(γ¯)时程曲线．当ν=0.1时，随机计算曲线平均值与确定性分析曲线几乎重合，随着ν的增大，平均挠度逐渐增大，土体的空间变异性对挠度的平均影响作用逐渐增大．10.13245/j.hust.238925.F006图6不同ν下路堤最大挠度的累计概率分布曲线(θ=5 m)10.13245/j.hust.238925.F007图7路面中心点平均挠度程曲线(θ=5 m)3.3　自相关距离的影响图8给出了ν=0.3，θ=1，10，50 m时，200次随机计算所得的加筋土路堤表面挠度曲线，黑色曲线为确定性分析结果，灰色曲线为随机计算结果，蓝色和粉色虚线为其中的两次计算结果．结合图4中θ=5 m的结果可见：随着自相关距离的增大，曲线分布的离散程度也相应增加，通过计算随机结果的挠度值与确定性结果的差异，得到3种自相关距离对挠度的最大影响分别为3.85%，20.83%，34.41%，这是因为自相关距离越大，在一定范围内的土体参数的相关性越大．在本模型中，当θ=50 m时，整个路堤的土体参数都在某一值附近波动，造成土体整体刚度同时增大或同时减小，因此位移曲线波动范围就大；而当θ=1 m时，土体之间的相关性变小，地基土中会同时出现高刚度区和低刚度区，多个区域之间会相互抵消各自对整体位移曲线的影响．10.13245/j.hust.238925.F008图8加筋土路堤表面挠度曲线(ν=0.3)可以看到：θ=1 m与θ=50 m时曲线形态基本上都与确定性分析相类似，挠度沿横纵截面的分布曲线基本沿路面中心对称；当θ=5，10 m时曲线形态发生了变化，挠度分布曲线不再沿路面中心对称．这两种现象产生于3种不同的原因：a．当自相关距离较小时，土体内部出现多个小范围的低刚度区和高刚度区，这些刚度区相互抵消了各自的影响，因此对整体曲线形态的影响较小；b．当自相关距离很大时，整个地基土体属于同种刚度，土体中不会出现较大的刚度差，分析情况类似于确定性分析，只不过是所确定的模量值有所不同，因此虽然量值会产生差异，但曲线形态不变；c．当自相关距离为1~20 m时，地基土区域内会形成较大区域的贯穿刚度区，使地基土中产生明显的刚度差，从而对挠度曲线的形态造成影响．图9给出了当ν=0.3时，不同θ下顶面最大挠度的累计概率分布曲线．随着θ的增大，随机分析得到的最大挠度分布曲线趋于分散，4条曲线相交于一点，当γk≥1.58 mm时，θ=50 m处超过控制挠度的概率最大；而当γk1.58 mm时，θ=1 m处超过控制挠度的概率最大．10.13245/j.hust.238925.F009图9不同θ下顶面最大挠度的累计概率分布曲线(ν=0.3)图10绘制了ν=0.3时，200次随机计算所得不同θ下路面中心点平均挠度时程曲线，随着θ的增大，平均挠度不断增大，土体空间变异性对挠度的平均影响程度逐渐增大，与图7对比可得，自相关距离的对挠度的平均影响程度小于变异系数．10.13245/j.hust.238925.F010图10不同θ下路面中心点平均挠度时程曲线(ν=0.3)4 结论以土工加筋路堤地基土的弹性模量变异性为切入点，借助随机有限元和统计分析方法，系统研究了地基土体弹性模量的变异系数和自相关距离对移动荷载下土工加筋路堤的变形特性的影响，主要结论如下．a．土体弹性模量的空间变异性对路堤的动力变形特征有着重要的影响，在研究路堤的动力问题中不应忽略．这一变异性对路堤变形的影响表现在变形量值和变形形态两方面．b．当自相关距离固定时，随着随机场变异系数的增大，随机计算得到的挠度曲线趋于离散，同时曲线形态逐渐发生明显变化．c．当变异系数固定时，随着随机场自相关距离的增大，随机计算得到的挠度曲线趋于离散．当自相关距离较小或较大时，随机计算对挠度曲线的形态影响较小；当自相关距离位于两者之间且在路堤实际尺寸以内时，随机计算对曲线形态有明显影响．d．随着变异系数和自相关距离的增大，最大挠度累计概率分布曲线均趋于离散，且弹性模量的变异性对挠度的平均影响程度均逐渐增大．