局部动荷载作用下土拱效应的离散元分析

2022-03-18李格烨沈盼盼张兴亚罗敏敏

长江科学院院报 2022年3期

关键词：模型试验填料幅值

李格烨，徐超，沈盼盼，张兴亚，罗敏敏

(1.同济大学土木工程学院，上海 200092； 2.中国长江三峡集团有限公司上海勘测设计研究院有限公司，上海 200434； 3.新城控股集团股份有限公司，上海 200092； 4.浙江大学建筑设计研究院有限公司，杭州 310027)

1 研究背景

土拱效应是土工结构中常见的现象。Terzaghi[1]通过Trapdoor试验证实了土拱效应的存在，并提出了“土拱效应”的概念。为了评价土拱效应的强弱，McNulty[2]定义了土拱率ρ，其表达式为

(1)

式中：σv为活动土体底部作用的平均竖向应力；γH为活动土体的自重应力，γ为活动土体的重度，H为活动土体的高度；q为土体表面作用的附加荷载。

自重荷载和表面静荷载作用下土拱效应的研究已有较多的成果，涉及土拱的变形特性、应力重分布规律及理论计算等。然而，在实际工程中，土拱不可避免地受到动荷载的影响。越来越多的学者对循环动荷载作用下的土拱效应开展研究。在模型试验方面，韩高孝等[3]以某桩-网结构断面为原型，开展了三维模型试验研究，结果表明在循环正弦波荷载作用下，先前由桩所承担的部分动应力会转移至桩间土，这意味着动荷载作用下土拱效应会发生退化。Xu等[4-5]利用自制的平面应变Trapdoor模型试验箱进行了一系列试验，结果表明，局部静动荷载均会弱化土拱效应，动荷载对土拱效应的退化程度大于静荷载。Bi等[6]采用砂土开展了Trapdoor试验，研究了均布循环荷载下土拱的结构形态、位移场、竖向应力的演变规律，再现了土拱从形成到完全塌陷的过程。Aqoub等[7]研究了循环荷载下低填方桩承式路堤的性能，表明土拱的破坏在荷载的前几个周期内发生，但是，路堤材料的振动密实以及路基软土的变形导致土体强度恢复，进而致使土拱效应部分或全部恢复。在数值计算方面，赖汉江等[8]通过离散元数值模拟，分析了动荷载作用下桩承式路堤中的土拱效应，结果表明，循环荷载作用下，低填方路堤中土拱效应先逐渐减弱，最终趋于稳定。梁自立等[9]采用 ABAQUS 有限元分析软件，分析了列车动荷载(循环正弦波荷载)对桩网结构路基土拱效应的影响，结果表明动荷载对土拱效应的影响与路基高度、桩间土性质密切相关；当路基内部形成完整土拱时，动荷载对土拱效应基本无影响，但会引起低矮路基中的非完整土拱发生失稳、退化。Pham等[10]利用有限元方法建立了桩承式加筋路堤三维模型，结果表明，循环荷载次数不同，土拱效应退化程度不同；交通荷载的速度越大，即荷载频率越大，土拱效应退化更明显。

上述研究成果集中于动荷载作用下土拱效应的退化甚至是破坏现象的认识，在作用机理、定量评价等方面有待进一步深入。本文基于Xu等(2019)[4-5]开展的室内Trapdoor模型试验，采用PFC2D建立了相应的离散元模型，研究了内摩擦角、荷载幅值和荷载频率对土拱效应的影响。

2 离散元模型

2.1 Trapdoor模型试验

Xu等(2019)[4-5]采用自制的Trapdoor模型试验箱开展了试验，其中填料采用横断面为圆形的铝棒相似土。铝棒相似土的单位重度为22.4 kN/m3，孔隙率约为17%。活动门的宽度B=128 mm。试验过程中先下移活动门至总位移达到8 mm，形成初始土拱，接着在填料表面施加动荷载。其他信息参见Xu等(2019)[4-5]。

2.2 离散元模型建立

2.2.1 试样制备方法

试样制备方法及过程决定了初始状态下试样中颗粒的孔隙及受力分布是否合理，因而对试验的结果有重要的影响。常用的试样制备方法有：半径扩大和重力沉积法(Bhandari等[11]、Han等[12])、改进的多层压缩法(Improved Multi-layer Compaction Method，IMCM)(Lai等[13])等。半径扩大和重力沉积法生成的试样内颗粒分布很不均匀，从而导致整个试样内都会出现较高的初始地应力。更关键的是，试样中的应力并不随高度呈现梯度变化，这种现象并不符合实际情况(Rui等[14])。IMCM试样制备的过程与模型试验过程相似，因此初始应力分布与实际大致相同。本文数值模拟基于PFC2D5.0软件，借鉴了IMCM的基本原理，提出了计算效率更高的多层收敛法，如图1所示。其主要过程如下：

图1 试样制备过程Fig.1 Flow chart for generating sample

(1)在指定区域内按照填料的孔隙率和质量比生成Ball颗粒层，如图1(a)所示，其颗粒接触如图1(b)所示，并对其进行求解直至收敛，结果如图1(c)所示。需要指出的是，在试样达到指定高度前，颗粒重力加速度g和表面摩擦系数μ均设置为0。

(2)在已经形成的颗粒层上方创建新的墙体形成新的空间，按照步骤(1)生成新的Ball颗粒层，如图1(d)所示，这时颗粒接触如图1(e)所示。删除2层颗粒间的墙体并再次计算至收敛，此时颗粒间接触消失，如图1(f)所示。

(3)重复第(2)步，直至达到试样的目标高度。

(4)对所有颗粒赋予重力加速度g及合理的表面摩擦系数μ，并计算至收敛。

2.2.2 Trapdoor 离散元模型

基于Trapdoor模型试验[4]，建立了如图2所示的Trapdoor离散元模型，其中，模型箱和活动门(B1—B7)均采用Wall单元模拟，填料利用2.2.1节中多层收敛法生成的Ball颗粒模拟。每次生成厚度为20 mm的填料(与模型试验相同)，直至试样高度达到2B。在活动门矩形区域上方有均匀分布的Ball颗粒来模拟物理模型中的垫层，垫层厚度h=22.4 mm。采用Clump单元模拟加载板，在Clump上施加动荷载模拟表面荷载。

图2 Trapdoor离散元模型Fig.2 Trapdoor discrete model

2.3 细观参数标定

填料的细观特性通过数值双轴试验标定，其细观力学参数列于表1。与物理试验相同，数值双轴试验的试样是将3、4、5 mm三种直径的颗粒以质量比1∶1∶1混合而成，采用50、100、150 kPa三级围压，试验所得偏应力及体积应变随轴向应变的变化关系如图3所示。数值双轴试验结果验证了文中颗粒填料细观参数取值的合理性，求解得出颗粒填料的宏观内摩擦角为24.8°，与实际双轴试验得出的结果25°十分接近。

表1 填料的细观力学参数Table 1 Micromechanics properties of fillings

图3 数值及物理双轴试验中偏应力和体积应变随轴向应变的曲线Fig.3 Curves of deviatoric stress and volumetric strain versus axial strain obtained from numerical simulation and laboratory biaxial tests

3 离散元模拟结果与分析

选取Xu等(2019)[4]中的正弦波CT4开展数值模拟，即填料高度H= 2B，加载板的长度L= 2B。数值试验过程与模型试验过程相同。

3.1 自重荷载下Trapdoor试验

在数值试验中，通过Fish语言编程记录活动门B4每下移0.1 mm时作用在活动门上的应力，绘制土拱率ρ随Trapdoor归一化位移(Trapdoor位移与Trapdoor宽度的比值)的变化关系曲线，即地基反应曲线(Ground Reaction Curve,GRC)。图4给出了模型试验和数值模拟对应的GRC。整体而言，模型试验与数值模拟中GRC的重合度较高，均分为初始拱、最大拱、应力恢复及极限状态4个阶段。最大拱对应归一化位移为1.25%，模型试验中的最大拱对应的土拱率，即最小土拱率是0.31，数值模拟中该值为0.32；极限状态开始阶段对应的归一化位移均>10%。因此，自重条件下数值模型较好地还原了物理模型。

图4 模型试验和数值模拟的GRCFig.4 GRCs obtained from model tests and numerical simulation

3.2 局部动荷载作用下Trapdoor试验

以幅值为8 kPa、频率为0.5 Hz，加载时间为200 s的正弦波(CT4)作为动荷载，施加在数值模型的填料表面。表2中，8 mm位移是指活动门B4下降8 mm，形成初始土拱，并未施加表面动荷载的阶段。8 kPa峰值是指在表面荷载8 kPa作用下，最后一个循环周期峰值所对应的土拱。8 kPa谷值是指最后一个循环周期谷值所对应的土拱。从表2可以看出，在局部动荷载加载阶段(荷载峰值)作用下，土拱效应有一定程度的退化。在卸载阶段(荷载谷值)，土拱效应有一定程度的恢复。对比模型试验和数值模拟的结果，可以发现各阶段土拱率均很接近。从应力角度来看，数值模型对物理模型的还原度较好。

表2 模型试验及数值模拟对应的土拱率Table 2 Soil arching ratios obtained from numerical simulation and model tests

图5(a)显示的是8 mm阶段所对应的模型试验和数值模拟的填料竖向位移云图。其填料的变形模式与塔形变形模式(Rui等[14])相同。图5(b)和图5(c)分别是8 kPa峰值和8 kPa谷值所对应的填料位移云图。如图5所示，模型试验和数值模拟的变形形态及变形大小二者均非常接近。从变形角度来看，数值模型很好地还原了物理模型。

图5 模型试验和数值模拟填料的竖向位移云图Fig.5 Vertical displacement contours of the backfill of model tests and numerical simulation

为了定量衡量动荷载作用下土拱效应的退化程度，定义了土拱退化参数α，如式(2)所示。

(2)

式中:ρload代表动荷载作用结束后的土拱率；ρinitial代表动荷载施加前的土拱率。当α=0时表明土拱效应在动荷载作用下没有退化；当α=100% 时表明土拱效应消失，即ρload=1。在8 kPa动荷载作用下，土拱退化参数α=64%，表明土拱效应退化了64%，即动荷载作用下土拱的发挥程度是初始土拱的36%。

填料颗粒间接触力是土拱效应的细观表现。图6展示的是加载前后接触力的变化规律，接触力线条颜色越深代表接触力越大。图6(a)是8 mm位移阶段接触力的分布情况，可以看出，稳定木块上方(B3和B5)的接触力明显大于活动门(B4)，且力链沿高度发生偏转，形成明显的拱形结构。图6(b)展示的是8 kPa谷值所对应的接触力分布，稳定木块上方(B3和B5)的接触力与活动门(B4)差别不明显，拱形结构基本消失，表明土拱效应退化明显。

图6 接触力分布Fig.6 Distribution of contact forces

4 参数分析

4.1 内摩擦角

通过改变表1中颗粒表面摩擦研究试样内摩擦角对土拱发挥程度的影响。颗粒表面摩擦系数μ=0.5、0.8、1.2对应填料的内摩擦角分别为24.8°、29.7°、32.6°。数值模型中填料高度H= 2B，荷载幅值为8 kPa，荷载频率为0.5 Hz,加载范围L= 2B。

图7展示了无荷载、动荷载峰值对应的土拱率与颗粒表面摩擦系数间的关系。填料的土拱率均随试样内摩擦角的增大而减小，表明土拱效应增强。但是，动荷载作用下土拱率的值均大于无荷载情况，说明在动荷载作用下土拱有一定程度的退化。土拱退化参数α与颗粒表面摩擦系数的关系如图7所示，随着颗粒表面摩擦系数的增大，土拱退化参数α逐渐减小，说明土拱效应的退化程度降低，抵御动荷载的能力增强。

图7 土拱率、土拱退化参数与颗粒表面摩擦系数间的关系Fig.7 Variations of soil arching ratio and degradation parameter of soil arching against surface frictional coefficient of particle

4.2 荷载幅值

数值模型中填料的高度H=2B和H=3B，动荷载加载频率为0.5 Hz，加载范围L=2B，颗粒表面摩擦系数μ=0.8。动荷载作用下H=2B和H=3B的填料土拱率随荷载幅值的变化(见图8)。荷载幅值为0对应活动门B4下降8 mm，形成初始土拱，并未施加表面动荷载的阶段。从图8可以看出，动荷载作用下填料的土拱率均随荷载幅值的增加而增大。在相同大小的动荷载作用下，高度较高(即H=3B)填料的土拱率明显小于高度较低(即H=2B)的土拱率，说明填料高度越高，土拱效应更明显。如图8所示，随着动荷载幅值的增大，土拱效应的退化越来越明显。相同荷载幅值作用下，高度较高(H=3B)的土拱退化参数α均小于高度较低(H=2B)的情况，如荷载幅值为12 kPa时，填料(H=2B)的土拱退化参数α为77%，H=3B的土拱退化参数α约为47%，前者是后者的1.6倍，说明增加填料高度可以明显减少动荷载作用导致的土拱效应退化。

图8 土拱率、土拱退化参数随动荷载幅值的变化Fig.8 Variations of soil arching ratio and degradation parameter of soil arching with cyclic loading amplitudes

4.3 荷载频率

选用4种频率的动荷载，分别为0.1、0.5、2.0、10.0 Hz，并与静荷载(频率为0 Hz)的结果进行了对比。数值模型填料的高度H=2B，颗粒的表面摩擦系数μ=0.8，荷载幅值为8 kPa，加载范围L=2B。

如图9所示，填料的土拱率随动荷载频率的增大而增大，表明局部荷载作用下土拱效应明显退化，并且动荷载作用下土拱效应的退化程度大于静荷载；当荷载频率为0时，土拱退化参数α值最小也说明了这一点。此外，图9表明随着荷载频率的增大，土拱退化程度越来越明显，在荷载频率为10 Hz时，土拱效退化参数α为72%，约是静荷载的2倍。