高周次列车荷载作用下桩网结构路基桩土应力比计算方法及特性研究

李 晖，毕宗琦，黄建丹，宫全美,王炳龙

(1. 同济大学 上海市轨道交通结构耐久与系统安全重点实验室，上海 201804；2.道路与交通工程教育部重点实验室，上海 201804； 3.中国铁道科学研究院集团有限公司 高速铁路轨道技术国家重点实验室， 北京 100081；4.北京铁科特种工程技术有限公司， 北京 100081； 5.上海市城市建设设计研究总院(集团)有限公司，上海 200125)

近年我国高速铁路线路大规模建成通车，随之而来的技术重点是针对高速铁路基础设施长期服役性能的研究，其中路基的沉降预测是长期性能保障的关键问题之一。高速铁路桩网结构路基的长期沉降预测与桩土应力比直接相关。目前桩土应力比的确定主要是基于静力学的土拱理论研究成果，而对列车动荷载作用特别是动荷载长期作用的研究尚不完善，主要原因之一是利用常规的计算方法难以实现动荷载高周次计算。因此，提出合理有效的方法对高周次列车荷载作用下桩网结构路基的桩土应力比进行研究，在满足长期变形和应力变化预测精度的基础上仍能保持较高的计算效率，减少计算成本，以适应高速铁路及其轨下基础在国内的发展规模，满足桩网结构路基承载性能设计及其长期服役稳定性预测的需要。

我国对高速铁路轨下基础使用寿命的要求达到100年[1]，若以其寿命时间内列车通行所产生的循环荷载次数衡量，其数量级将达到百万甚至千万级，其间列车荷载的动力效应无法忽视[2]。而目前桩网结构路基设计时对动荷载的考虑基于拟静力法，无法考虑长期动荷载作用的累积效应。另一方面，依靠既有的动力加载室内试验或常规数值模拟方法需要大量的时间和成本，难以完整且高效地模拟这一累积过程，同时面临着循环次数增加导致的计算误差累积问题。因此，目前针对高周次动荷载作用下桩网结构路基中桩土应力比的研究较少。陈仁朋等[3]通过桩网加筋路堤足尺模型试验，研究速度324 km/h时列车120万次轮轴动荷载作用下，桩帽高度位置和3.2 m路基范围内动土压力随振次的变化规律，分析长期动荷载作用下土拱演化规律。韩高孝等[4]以4根桩为单元进行三维模型试验，研究在动荷载幅值为10～20 kPa、循环次数达18万次作用下的桩土应力变化。吕玺琳等[5]通过长期循环动态加载模型试验，研究最高循环次数达到 4 万次时路基沉降、桩土应力分担及桩身轴力分布随循环次数的变化规律。上述学者对高周次列车动荷载作用下的路基结构或桩土应力的研究取得了诸多有价值的结论，但主要限于试验规律的探索，缺乏与之对应的理论计算方法。

国内学者通过数值模拟，对桩网结构路基的研究主要集中在少量振次作用下桩土应力比的变化规律，如文献[6-8]等。现有有限元实现高周次加载的方法中，具有代表性的有张宏博[9]提出的基于安定性理论的砂土双硬化弹塑性本构模型、贾鹏飞等[10]的高周次荷载作用过程中土体塑性包络的弹塑累积模型、德国Ruhr-Bochum大学的HCA模型[11]以及康庄的显隐式转换计算方法[12]。总体上，前三种方法主要针对砂土，科学性及理论性强；而显隐式转换计算方法针对长期列车荷载和软土路基特性进行了拓展，具有计算效率高、工程应用性强的特点，但其中的经验公式存在一定的局限性。

为突破上述实现高周次加载方法的局限性并研究桩网结构路基在相应列车荷载作用下的服役性能，本文利用自主研发的显隐式转换计算方法实现高周次加载，利用有限元Abaqus的二次开发平台UMAT实现桩间土、路基土的长期变形特性模拟；该算法考虑了土体在长期动荷载作用下的累积效应，避免了系统误差累积导致的计算结果失真问题；并通过京沪高速铁路徐沪段桩网结构路基断面模型与实测结果的对比验证了该算法的可靠性；最后，探讨了在高周次列车荷载作用下动荷载次数、动荷载幅值等因素对桩网结构路基桩土应力比的影响规律。

1 高周次有限元加载方法的实现

1.1 基本思路

(1)

本文中实现高周次加载的显隐式转换计算方法在德国HCA模型的基础上，将计算过程分为隐式计算过程和显式计算过程两部分，并在有限元中对其进行应用和深化。

如图1所示，隐式计算过程指考虑每一次循环荷载加载与卸载作用于土样单元上所产生的应力-应变作用。显式计算过程则是在循环荷载下土体永久变形的研究中，忽略土体的弹性变形，基于经验公式直接求取土体永久塑性变形曲线的有限元计算过程。

图1 显隐式转换计算方法的思路

在高周次加载的有限元计算中，单独采用隐式计算过程易造成系统冗余和成本过高，且产生的系统误差累积将导致最终计算结果失真。显式计算过程的应用性强，能反映高周次荷载下的土体变形特性，但在高周次下的计算精度相对较低，尤其在万次后；也存在计算所需的试验参数过多，实际工程无法直接采用，在通用有限元平台中应用性不足等问题。因此采取在显式计算过程中穿插隐式计算的策略，通过后者对前者的计算参数进行校核修正，在利用显式计算过程对高周次加载计算优势的同时兼顾了隐式计算的精度。

1.2 本构模型的选取与建立

本文隐式计算过程选用摩尔库伦弹塑性本构模型；显式计算过程的模型采用Li修正公式。算法在有限元平台的实现则借助Abaqus中蠕变计算模块进行二次开发。选取欧拉超前法求解蠕变应变，即

(2)

采用时间硬化模型，等效蠕变应变率为

(3)

积分可得等效蠕变应变计算公式

(4)

对比Li修正经验公式可以看出两者形式相似，式中的循环次数N与式(4)的作用时间t可以通过荷载作用频率建立关系，进而相互转化。

1.3 显隐式转换方法加载次数的判定准则

隐式计算过程需为后续的显式计算过程提供应力场，并在显式计算过程中适时打断，更新应力场。为避免起始点不规则循环及数据冗余带来的误差，在每一步的隐式计算过程中宜取3～4次循环进行加载过程的计算。

考虑显式计算过程本身的封闭计算特性，单步循环荷载次数Nn会影响该阶段的计算应变值。故根据试验数据及前后计算步应变关系，利用第一容许误差εtol,f和后续容许误差εtol,n确定每步显式计算过程的合理加载次数。

(1)第一容许误差εtol,f

设第1步显式计算过程中的循环荷载次数为N1，有限元显式计算过程所得累积塑性应变与实测累积塑性应变之间的误差为εtol

(5)

(2)后续容许误差εtol,n

在一定循环荷载加载次数后，为保证其后判断的可信性，需寻找误差参照系与试验数据之间的关联性。后续容许误差εtol,n界定了第n步显式计算过程所取的循环荷载次数Nn的限制条件。εtol,n可按以下递推关系参考取值

(6)

(7)

2 铁路路基计算断面的选取与计算结果验证

2.1 路基断面模型的选取

显隐式转换计算方法对复杂的路基结构是否适用仍需证明，因此选取京沪高速铁路徐沪段某断面现场原位1∶1动载模拟试验[13]的路基横断面为计算断面，见图2。线路为板式轨道，地基处理方式为桩网结构，路基高度4.15 m，其中基床表层0.4 m，碎石垫层厚0.6 m，垫层内铺设一层土工格栅。桩为端承桩，桩径0.5 m，桩帽直径1 m，桩间距1.8 m，正方形布置，持力层为强风化角闪岩。地层分布从上而下依次为①黏土(6 m)和②2强风化角闪岩。现场试验分别在桩顶、桩间土中心及形心对应路基不同高度处布置传感器，包括静态土压力盒、动态土压力盒和沉降板各6层，静态土压力盒和沉降板用于测试路基填筑过程中静荷载和沉降的变化，动态土压力盒安装沉降板竖杆的位移计用于测试动载试验时路基的动应力和动变形，详细布置见文献[13]。

图2 现场试验路基横断面图(单位：m)

结合上述路基断面建立二维有限元模型，有限元模型的路基底部宽度外地基延伸宽度取15 m，深度延至地面下15 m，并在模型右边和底部设定无限元边界。土工格栅与垫层土之间采用嵌入接触方式，桩周、桩顶和桩底与土体的接触均采用面对面离散接触对的接触关系，桩土接触面法向接触模型采用硬接触，切向接触模型选用库伦摩擦接触，模型中土层间的接触均采用绑定接触。土工格栅釆用杆单元truss模拟，其余各部件均选用4节点缩减积分实体单元，混凝土桩均釆用线弹性材料模型。数值模型见图3。

图3 数值模型图(单位：m)

2.2 模型参数

2.2.1 材料参数

假定各种材料都是均匀、各向同性体，荷载作用过程中不考虑地下水的影响，不考虑加载过程中土体强度的变化且土工格栅强度随应变增加而减小。数值模型中仅考虑平面应变问题，将桩简化为桩墙，对其弹性模量采用式(8)进行折减。

(8)

式中：Esp为桩墙的弹性模量;Ep为桩的弹性模量；Es为土体弹性模量；D为桩径；l为桩间距。

有限元模型中各组成部分材料参数如表1所示。

表1 材料参数

2.2.2 本构模型选取

计算模型中，路基土、桩间土在隐式过程中釆用弹塑性材料，服从摩尔库伦屈服准则，在显式过程中采用Li修正经验公式

(9)

式中：εp为累积塑性应变，%;σs为土体静破坏偏应力;σd为循环荷载引起的动偏应力，采用式(10)进行计算;N为循环加载周期数;a为各土壤类型对应的材料参数。

(10)

式中：σ1和σ3为最大和最小主应力;J2为第二偏应力不变量。

静破坏偏应力计算式为

σs=2qf

(11)

式中：qf为不排水抗剪强度，可以利用固结不排水的总应力强度指标ccu和φcu来计算[14]，即

(12)

其中，K0为土的侧限系数;σz为土体的自重应力。文献[15]通过现场激振试验给出了Li修正公式中级配碎石、路堤填土、碎石垫层及底层土的参数a、m、b，如表2所示。本文路基土参数取值参考文献[15]。

表2 现场试验[15]推算得出的Li模型中材料系数a、m、b

2.2.3 土工格栅弹性模量设定

土工格栅的拉力和应变随时间变化过程见图4[16]。在Abaqus中使用USDFLD子程序模拟土工格栅弹性模量随应变增加而变小的特性。

图4 格栅拉力和应变随时间发展示意(抗拉刚度1 000 kN/m)

2.3 列车荷载模拟

有限元计算所施加的动荷载大小与现场动态模拟试验中相同，其动应力幅值处于10～20 kPa范围之间，频率为10 Hz，采用正弦荷载。荷载施加范围与现场试验的激振器下方承台宽度一致。

2.4 模型验证

2.4.1 施加动荷载前的应力分析

仅在路基填土荷载作用下，有限元计算得到的路基中竖向应力云图见图5。

图5 竖向应力云图

由图5可得，路基底部存在明显应力集中现象，与文献[17]中路基填土作用下土体的竖向应力云图进行对比，两者规律相近。在桩顶上方一定范围内，竖向应力出现应力集中现象，应力集中区域呈圆拱状，中心竖向应力约为桩侧未受影响土体中应力的6倍，路基中心处竖向应力集中较大，边坡处较小。

按图5的虚线处提取桩顶、桩间土上方的竖向应力，与现场数据对比，得到图6。有限元和现场试验的竖向应力沿深度分布规律相近，数值基本相同。

图6 竖向应力沿深度分布

2.4.2 施加550万次动荷载的应力分析

为与现场原位试验进行的550万次激振试验结果进行对比分析，利用显隐式转换方法对上述路基断面不同位置处动应力、桩土应力比随循环荷载次数的变化进行计算，结果如下。

(1)动应力随振次的变化

现场试验和有限元得到的桩顶、桩间土上方的动应力随动载作用次数的关系见图7。

图7 现场试验与有限元的桩顶、桩间动应力和动载次数关系

从图7看出，现场试验的桩顶上方动应力集中在3.5 kPa，桩间土上方动应力从3 kPa减小到2.5 kPa。有限元计算的规律与现场试验类似。

(2)动应力在深度方向的变化

有限元与现场试验的桩顶、桩间动应力沿路基深度分布见图8，图中动应力是动载施加50万次结果。

图8 动应力沿路基深度的分布结果

从图8看出，动应力随深度衰减，桩间土上方动应力小于桩顶，有限元计算结果的动应力随深度变化的规律与现场试验类似，所测得的数值差异较小。

(3)桩土应力比在不同循环次数下的变化

图9是现场试验和有限元计算的桩土应力比与作用次数的关系曲线。有限元计算结果比现场试验略大，但整体上有限元计算与实测值规律较为统一，且量值差异较小，证明计算方法的可靠性，可用于后续的影响因素研究。

图9 桩土应力比与作用次数关系曲线

3 高周次动荷载对桩土应力比的影响

在前述路基高度、地基处理参数以及动荷载频率不变的情况下，分别针对动荷载作用次数、动荷载幅值等影响进行研究。为方便后续对比，此处采用无量纲参数动静荷载比η表示

(13)

式中：σd为在路基面施加的动应力；h为路基高度；γ为填土重度。

3.1 动荷载对大主应力的影响

施加不同动荷载幅值后，路基填土中的大主应力云图见图10。

通过图10施加动荷载前后对比看出，随着动荷载施加，桩顶上方应力增大区域减小，土拱高度下降，表明动荷载作用下，桩土上方的应力重新调整，一部分桩顶上方的应力逐渐向桩间土上方转移，土拱效应发挥的作用被削弱，与韩高孝等[4]的试验结果及魏平等[7]的数值模拟结果相符。桩顶上方大主应力受动荷载影响区域，横向为路基中心到混凝土底座边缘外3 m范围，即从混凝土底座边缘正下方以1∶1的坡度扩散至路基面。随着η的增加，受影响范围内桩顶上方应力增大的区域不断变小，表明土拱效应不断减弱；当η=1.67时，桩顶上方的大主应力云图形状变化最明显，桩顶上方应力增大的区域最小，说明动静荷载比η越大，土拱效应减弱越明显。

3.2 动荷载对竖向应力的影响

提取桩顶、桩间土上方的应力(如图5虚线处)，得到不同动荷载幅值作用后，桩顶及桩间土上方竖向应力沿深度的分布，见图11。

图11 桩顶及桩间竖向应力沿深度的分布(10万次动荷载作用后)

从图11可以看出：(1)动荷载作用后原有的虚拟土拱效应仍存在，但土拱效应的范围减小。(2)在相同振次下，随荷载比η的增加，桩顶及桩间土的竖向应力分叉点逐步下移。(3)η为0.33、0.67时，分叉点位置变化较小；η≥1时，分叉点下移的程度增大。(4)η=1.67时，桩顶及桩间土竖向应力沿深度分布趋于相同，可认为土拱效应基本消失；但是由于土工格栅的存在，使得仍有部分路基填土的荷载和动荷载传递到桩上，所以分叉点依旧存在，此时分叉点的高度与土工格栅所在位置接近。

3.3 动荷载对桩土应力比的影响

不同荷载比作用下桩土应力比的变化曲线见图12。从图12看出：(1)随着循环次数的增加，桩土应力比先降低后逐渐稳定。当η为1、循环次数小于10万次时，桩土应力比降低较快，超过10万次后，其降低幅度变小，逐渐趋于一定值。说明动荷载作用次数有一定敏感范围，这一敏感范围随动静荷载比η而变化。(2)动静荷载比较小时，动荷载对桩土应力比的影响可以忽略。如η为0.17时，桩土应力比变化可忽略；随着动静荷载比η的增加，桩土应力比随振动次数减小的梯度增加，如η分别为1、1.33时，桩土应力比的变化率分别达到32%、40%。(3)动静荷载比增加到一定值，如η≥1.67时，5万次动荷载作用后，桩土应力比降低为1.1，桩顶及桩间竖向应力沿深度分布趋于相同，表明土拱效应消失。

图12 不同荷载比η时桩土应力比的变化曲线

同时引入土拱效应稳定系数β来评价高周次荷载作用下桩土应力比变化对土拱效应稳定性的影响。定义土拱效应稳定系数β为

(14)

式中：n后为动荷载作用后的桩土应力比；n前为动荷载作用前的桩土应力比。得到高周次动荷载作用下土拱效应稳定系数β与振次、η的关系，见图13。

图13 高周次动荷载作用下土拱稳定系数β与振次、η的关系图

从图13看出，当桩土刚度比Kp=14时，η=1.33、η=1.67把图划分成土拱效应稳定区、过渡区、消失区三个区域。其中，土拱效应稳定区是指动荷载作用后，桩土应力比n后有一定减小，但在一定作用次数后保持不变，土拱效应的稳定性有所降低，但能长期保持稳定；土拱效应消失区是指桩土应力比n后迅速减小至接近1，此时土拱效应消失。过渡区是指在稳定区与消失区之间存在临界动静荷载比η。

通过上述研究，在本文计算条件(路基高度3 m，桩间距1.8 m，桩土刚度比Kp为14)下，当η≥1.67时，5万次荷载作用后，桩土应力比迅速减小，可能会威胁土拱的稳定性。

4 结论

本文基于HCA模型，借助Abaqus有限元平台的二次开发实现了显隐式转换算法，并通过现场测试数据进行验证，为高周次加载研究提供了一种高效简便的有限元计算思路。利用该算法，研究了高周次列车荷载作用下动荷载作用次数、幅值等因素对高速铁路桩网结构路基中桩土应力比的影响。根据本文模型分析结果，得到以下结论：

(1)高周次动载影响下，桩土应力比随振次的增加呈现先减小后趋于稳定的规律，伴随着桩顶上方应力集中区域的减小、竖向应力分叉点的下移，其影响范围与桩顶上方大主应力受动荷载影响范围一致，即从混凝土底座边缘正下方以1∶1的坡度扩散至路基面。

(2)动荷载幅值大小影响着桩土应力比随循环次数的变化规律，动静荷载比大于1.67时，桩土应力比迅速减小至土拱效应的消失，动静荷载比小于0.33时桩土应力比保持稳定。

(3)高周次动载影响下，桩顶及桩间土的竖向应力分叉点随动静荷载比的增加呈现逐步下移的规律。动静荷载比η小于1时，土拱高度随η增加而平缓减小；当动静荷载比η大于1时，土拱高度随η增加而急剧增大。

(4)根据动静荷载比η，可将高周次动荷载作用下土拱效应稳定系数β与振次、η的关系图划分成土拱效应稳定区、过渡区、消失区三个区域，且在稳定区与消失区之间存在临界动静荷载比η。