列车循环荷载下水泥土复合地基沉降的颗粒离散元计算

关振长，廖重辉，罗志彬，郭光照

(1.福州大学土木工程学院，福建福州350116;2.中国土木工程集团福州勘察设计研究院有限公司，福建福州350001)

列车循环荷载下水泥土复合地基沉降的颗粒离散元计算

关振长1，廖重辉1，罗志彬1，郭光照2

(1.福州大学土木工程学院，福建福州350116;2.中国土木工程集团福州勘察设计研究院有限公司，福建福州350001)

依托福州江阴港铁路支线水泥土搅拌桩软基处理试验工程，探讨利用颗粒离散元方法计算水泥土复合地基在列车循环荷载作用下的工后沉降问题。首先根据原场地淤泥质土的三轴压缩试验和水泥土的单轴压缩试验所得的应力—应变曲线，利用颗粒离散元建立三维数值模拟模型，并结合统计学习方法，对复合地基中的淤泥质土和水泥土颗粒的细观参数进行反演;然后按江阴港水泥土搅拌桩复合地基的实际尺寸建立缩尺的离散元数值模型，进行列车荷载下的沉降计算。结果表明:在500～1 000趟列车通过后水泥土搅拌桩复合地基的工后沉降趋于稳定，沉降量约为20 mm。这一结果与按现有规范计算所得沉降量较为接近。

水泥土搅拌桩 列车循环荷载 工后沉降 颗粒离散元

水泥土搅拌桩是一种利用水泥作为固化剂，就地将原位土和水泥浆液强制搅拌，使原位土硬结成具有整体性、水稳性和一定强度的水泥加固土的软基处理方法。该法因成本低廉、形式灵活和对环境影响小等优点，在软基处理工程中得到广泛应用。国内外许多学者通过室内试验、现场测试及数值模拟等方法，对水泥土搅拌桩的力学特性、加固机理及设计方法进行了大量的研究［1-3］。但受桩土相互作用的复杂性以及研究手段的制约，已有研究主要局限于考虑水泥土复合地基的承载力，而甚少探讨它的工后沉降，尤其是未能考虑在列车循环荷载作用下水泥土搅拌桩复合地基的工后沉降问题［4-6］。

对于复合地基的工后沉降分析，常面临如何模拟桩土接触这一难题。一般而言，对于刚性桩复合地基，在有限元等常规的数值方法中可采用接触面单元加以解决，其桩土接触面参数可参考室内试验或现场测试结果来进行选取［7-8］。但水泥土搅拌桩是原位搅拌固化成型的半刚性桩，受成桩过程中强力搅拌扰动的影响，搅拌桩与原位土之间的接触面常常并不明显，若简单套用刚性桩的桩土接触面参数，将给沉降计算带来误差甚至错误。反之，若不考虑桩土接触问题，即假设水泥土搅拌桩和原位软土始终共同变形，又与实际明显不符。另外，若要深入探讨水泥土搅拌桩复合地基在列车循环荷载作用下的工后沉降机理，还须设定水泥土与原位软土的动力本构及相应的动参数，而这方面可资利用的研究成果也甚为少见。因此，常规的数值方法尚存诸多不足和局限，有必要尝试新的思路和计算方法。

自上世纪七十年代Cundall提出适用于岩土力学的颗粒离散元方法以来，在学者们的推动下，该方法已成功应用于许多室内岩土试验的数值模拟与机理研究，并得以在岩土和地下工程领域广泛应用。与采用有限元等常规方法计算水泥土复合地基沉降的不同之处在于，在颗粒离散元方法中可设定水泥土颗粒、原位软土颗粒的接触本构及相应的接触参数，通过细观层面上不同颗粒之间的接触力与位移变化，来实现对水泥土复合地基宏观力学性能的模拟。此法可巧妙绕开半刚性桩复合地基中桩土接触模拟和动力参数选取的难题，同时在计算中采用显式时间差分动态求解技术，很适用于列车循环荷载之类的动力问题研究［9-12］。

福州江阴港铁路支线为Ⅱ级货运铁路，其线路穿越地区上覆有16～20 m厚的滨海相欠固结淤泥质软土层。线路下设计采用桩长12 m的水泥土搅拌桩进行地基加固，桩径0.5 m，梅花形布置，水泥掺量12%。本文以福州江阴港铁路支线水泥土搅拌桩软基处理工程为背景，尝试利用基于PFC3D数值平台的颗粒离散元方法对水泥土搅拌桩复合地基在列车循环荷载作用下的工后沉降问题进行深入探讨。

1 复合地基细观参数的反演

1.1 淤泥质土与水泥土的基本物性

在搅拌桩加固施工前，对场地淤泥质土进行了6个孔位的原位钻探取样，进行了常规土工试验和三轴压缩试验，试样直径D×高H=3.9 cm×8 cm，得到淤泥质土的基本物理性能指标如表1所示。

表1 加固处理前淤泥质土的物理性能指标

在水泥土搅拌桩成桩90 d后，对其中18根水泥土搅拌桩进行全长抽芯取样，进行了常规土工试验和单轴压缩试验，试样尺寸D×H=5 cm×10 cm，得到水泥土的基本物理性能指标如表2所示。

表2 加固处理后水泥土的物理性能指标

1.2 颗粒接触本构及参数

在颗粒离散元数值分析中，颗粒的接触本构及其参数的选择至为关键。对淤泥质软土和水泥土这两种具有黏聚力的材料选用粘结型颗粒接触本构是较为合适的。该本构模型假设:离散颗粒间的接触仅发生在接触点位置;颗粒间通过法向弹簧传递法向力，通过剪切弹簧传递剪切力;接触力的大小由颗粒的法向刚度kn、切向刚度ks、法向抗拉强度nb、切向剪切强度sb和颗粒间的摩擦系数μ这5个参数共同决定。当任意两个颗粒A与B相互接触时，总法向刚度Kn与总切向刚度Ks可由下式得到

在计算中，离散颗粒之间的法向接触力Fn与切向接触力Fs可根据上一计算步求得的颗粒间相对法向位移un与切向位移us由式(3)和式(4)求得。需要注意的是，若计算所得法向接触力Fn超过了阈值nb，可认为是法向黏结破坏，Fn被重置为0;若切向接触力Fs超过了阈值sb+μFn，则认为是切向黏结破坏，Fs被重置为μFn。

本文将根据室内试验所得应力和应变数据，通过颗粒离散元数值模拟与统计学习相结合的方法，分别反演得到淤泥质土颗粒和水泥土颗粒的kn，ks，nb，sb和μ这5个细观参数。

1.3 淤泥质土三轴压缩过程的数值模拟

本文以PFC3D作为数值分析平台，对淤泥质土UU三轴压缩试验进行全过程动态仿真模拟，具体思路和步骤如下:

1)根据淤泥质土三轴压缩试样的实际尺寸(即D ×H=3.9 cm×8.0 cm)，构建一个圆柱形的刚性模型箱，如图1所示。在箱内随机填充不同尺寸的圆形颗粒，其直径在0.8～1.2 mm之间均匀分布，以模拟淤泥质土试样的初始状态。根据常规土工试验测得的淤泥质土试样孔隙率，计算填充所需的颗粒数目，以避免填充过密或过疏。

图1 淤泥质土试样UU三轴压缩过程的数值模拟示意

2)设定淤泥质土颗粒之间的一组接触本构参数，即颗粒的kn，ks，nb，sb和μ 5个参数，在自重条件下计算至初始平衡。

3)将模型箱周围的圆柱形墙体设置为按应力伺服控制，即始终保持一定的围压。在每一计算步前，先求得各颗粒与模型箱墙体之间的平均接触应力，再控制本计算步内圆柱形墙体的外扩或内缩速度，以保证本计算步后颗粒与模型箱墙体平均接触应力即围压保持不变。

4)将墙体模型箱的上下底面墙体设置为按定常速度进行控制，即每一计算步压缩变形量相等，总计算步按试验过程中试样的总压缩量控制;同时记录模拟过程中土样中部的平均轴向应力，绘制主应力差与轴向应变曲线。此处平均轴向应力是指特定体积内所有颗粒间接触力之和在轴向上的集度。

1.4 基于统计学习方法的参数反演

本文采用数值模拟与基于神经网络的统计学习相结合的方法，对土颗粒的5个细观参数进行反演，即要找到一组最佳参数组合，使得由PFC3D数值模拟得到的轴向应力差与轴向应变曲线在最大程度上与试验结果相吻合［13-14］。参数反演的基本过程为:

1)利用误差反向传播神经网络构建一个统计学习机。学习机为一个5×5×5的神经网络，即输入层对应五维的输入向量(含5个神经元)，输出层和中间单隐藏层也设定为5个神经元;

2)根据计算经验分别确定淤泥质土5个细观参数的取值范围，并在各参数的上下限范围内随机取值生成一组参数组合，即输入向量;

3)将该输入向量代入PFC3D数值模型进行计算，记录轴向应变为1%，3%，5%，7%，9%时的主应力差，即输出向量;

4)把这样一组输入输出向量对作为一对学习样本，重复以上步骤共生成若干对学习样本;

5)将以上学习样本输入到统计学习机中，通过不断训练和调整自身权重，最终在统计意义上用学习机的映射过程取代PFC3D的数值计算;

6)利用训练完毕的学习机，结合遗传算法，不断优化输入向量，最终寻得一个最佳输入向量，使得其对应的输出向量尽可能接近试验结果;

7)将寻得的最佳输入向量再次代入PFC3D数值模型中，验证其对应的输出向量是否接近试验结果;

8)若相差较大则舍弃之，并按以上方法重新寻找一个最佳输入向量，直到满足精度要求为止。

1.5 淤泥质土细观参数的反演

利用上述方法，以深度5.5 m处的淤泥质土试样为例，在围压100 kPa下对试样的细观参数进行反演，最终得到一组最佳细观参数为:ks=49.1 MPa，kn/ks= 2.01，nb=4.9 kN，sb/nb=9.8，μ=0.05。将其代入PFC3D数值模拟中进行验证，并绘制三轴压缩试验和数值模拟所得的应力应变曲线，如图2所示。由计算结果可知，两条曲线在轴向应变为1%，3%，5%，7%，9%处的平均误差仅为3.4%，满足精度要求。

同样，对深度16 m处的淤泥质土试样，在300 kPa围压下利用数值模拟与统计学习相结合的方法，反演得到深部淤泥质土颗粒的一组最佳细观参数为:ks= 48.0 MPa，kn/ks=2.05，nb=5.1 kN，sb/nb=9.4，μ= 0.06。将其代入PFC3D数值模型中进行验证，绘制的应力应变曲线如图3所示。与试验曲线相比，数值模拟平均误差仅为2.9%，满足精度要求。

图2 淤泥质土试样UU三轴压缩过程的应力—应变曲线(100 kPa)

图3 淤泥质土试样UU三轴压缩过程的应力—应变曲线(300 kPa)

由反演结果可知，不同深度处淤泥质土颗粒的细观参数比较接近，这与原场地内淤泥质土均匀性较好的实际情况是一致的。

1.6 水泥土细观参数的反演

对于成桩90 d后的水泥土，也可进行类似的细观参数反演。以深度为2.0 m和8.0 m处的芯样为例，反演得到浅部水泥土颗粒的一组最佳细观参数为:ks=69.4 MPa，kn/ks=1.45，nb=13.6 kN，sb/nb=3.7，μ=0.75;深部水泥土颗粒的一组最佳细观参数为: ks=53.0 MPa，kn/ks=1.58，nb=8.1 kN，sb/nb=4.2，μ=0.63。

将此2组细观参数分别代入PFC3D数值模型进行验证，得到的应力应变曲线分别绘制于图4和图5中。与试验曲线相比，数值模拟平均误差分别为4.8%和5.3%，满足精度要求。

图4 浅部水泥土芯样单轴压缩过程的应力—应变曲线

图5 深部水泥土芯样单轴压缩过程的应力—应变曲线

由反演结果可知，不同深度处水泥土颗粒的细观参数差别较大，这与浅部软土水泥搅拌桩加固质量较好，而深部加固质量较差的实际情况相吻合。

2 列车循环荷载作用下的工后沉降

2.1 水泥土搅拌桩复合地基的离散元数值模型

在上文对淤泥质土和水泥土的细观参数进行反演时，颗粒粒径设定为0.8～1.2 mm之间均匀分布。若按这一粒径来模拟实际尺寸的水泥土搅拌桩复合地基(桩径0.5 m，桩长12 m，桩间距1.2 m)，对其进行沉降分析，则需要上亿个离散元颗粒，这大大超出了微型工作站的运算能力。因此本文根据相似理论进行缩尺的数值模拟［15］。由量纲分析法获得的各物理量的相似比如表3所示。根据前人研究结果［7-10］，当模型最小几何尺寸与颗粒最大几何尺寸的比值≥10时，可忽略颗粒粒径对数值模拟结果造成的尺寸效应，因此缩尺数值模拟中n的取值为100。即水泥土搅拌桩复合地基的缩尺数值模型，其内部水泥土搅拌桩的直径为5 mm，深度为120 mm;其周边淤泥质土的直径为12 mm，深度为180 mm。淤泥质土和水泥土颗粒的5个细观参数按上节选取，但重力加速度根据相似关系取为980 m/s2。

表3 缩尺数值模拟的相似比

2.2 列车循环荷载的施加

为便于研究，本文根据文献［16］将列车竖向荷载表示为

式中:P0为列车恒载，根据《铁路路基设计规范》(TB 10001—2005)，江阴港Ⅱ级货运铁路支线P0取为60.1 kPa;f为荷载频率，根据文献［16］，当列车速度低于250 km/h时，其振动频率主要集中在4 Hz左右。根据相似关系，数值模拟中所采用的荷载频率为400 Hz。

由此作用于复合地基顶面的荷载

图6 水泥土搅拌桩复合地基工后沉降的经时变化

2.3 工后沉降计算

按500 000个计算步进行模拟，每个计算步代表的模拟时间为0.000 21 s(此为程序根据颗粒数量、应力水平等因素自行设定的计算参数)，则相当于进行了约105 s的列车循环加载的动态模拟。由此所得水泥土搅拌桩复合地基的工后沉降与模拟时间的关系曲线如图6所示。由该图可以看出，经历模拟时间100 s后(按相似关系即实际时间为10 000 s后)，水泥土搅拌桩复合地基的工后沉降趋于稳定，其中搅拌桩桩顶沉降为0.202 mm，按相似比其实际沉降为20.2 mm;搅拌桩中部6 m深度处沉降为0.107 mm，实际沉降为10.7 mm，搅拌桩底部12 m深度处的沉降仅为0.01 mm，实际沉降仅为1 mm。

需要说明的是，实际时间是指列车荷载实际作用在搅拌桩复合地基上的时间，假定一趟列车通过某处复合地基的所需时间为10～20 s，那么实际时间10 000 s则相当于有500～1 000趟列车通过。

2.4 应力面积法的沉降计算结果

为便于比较，根据《铁路路基设计规范》(TB 10001—2005)和《建筑地基基础设计规范》(GB 50007—2011)，计算了场地水泥土搅拌桩复合地基在列车荷载下的沉降。将列车视为宽度3.5 m，集度60.1 kPa的条形荷载，则水泥土搅拌桩复合地基的压缩模量Esp

式中:m为置换率，根据设计取为0.136;Es为淤泥质土的压缩模量，根据室内土工试验结果取为3.4 MPa; Ep为水泥土的压缩模量，根据芯样的单轴压缩试验，其平均值取为50 MPa。

经计算求得复合地基压缩模量Esp为9.69 MPa。最后采用应力面积法，得到福州江阴港铁路支线水泥土搅拌桩复合地基在列车荷载作用下的总沉降量为14.1 mm，与上节的计算结果较为接近，这表明采用颗粒离散元方法进行沉降计算是可行的。

3 结论

针对有限元等常规数值方法在计算水泥土复合地基沉降时遇到的桩土接触参数和动力本构参数选取等诸多不易解决的问题，本文提出了利用颗粒离散元方法研究水泥土搅拌桩复合地基在列车循环荷载作用下的工后沉降计算问题的整体思路和方法:

1)首先根据室内地基土的试验数据，对淤泥质土颗粒和水泥土颗粒的细观参数进行反演计算;

2)再按相似原理，对水泥土搅拌桩复合地基建立缩尺的数值模型;

3)最后将反演所得的细观参数应用于缩尺模型的数值模拟，求得水泥土搅拌桩复合地基工后沉降的经时变化。

针对福州江阴港Ⅱ级货运铁路支线软基处理的实际，通过离散元数值模拟得到在500～1 000趟列车通过后，场地水泥土搅拌桩复合地基的工后沉降约20 mm，并趋于稳定。这一结果与采用传统的应力面积法计算所得的工后沉降值较为接近，证明了采用颗粒离散元方法研究水泥土搅拌桩复合地基在列车循环荷载作用下的工后沉降问题的可行性。

［1］HORPIBULSUK S，CHINKULKIJNIWAT A，CHOLPHATSORN A，et al.Consolidation Behavior of Soil-cement Column Improved ground［J］.Computers and Geotechniques，2012，43(1):37-50.

［2］张伟丽，蔡健，林奕禧，等.水泥土搅拌桩复合地基荷载传递机理的试验研究［J］.土木工程学报，2010，43(6):116-121.

［3］向玮，刘松玉，经绯，等.深长变径搅拌桩荷载传递规律的试验研究［J］.岩土力学，2010，31(9):2765-2771.

［4］赵春风，李永刚，钱涛.水泥土搅拌桩复合地基桩土应力比的解析算法［J］.中南大学学报(自然科学版)，2012，43 (6):2390-2395.

［5］MAEDA Y，ICHIKAWA A，TSURUKUBO S，et al.Study on Bearing Capacity Characteristics and Practical Application of Composite Ground Pile Foundation［J］.Japanese Geotechnical Journal，2008，3(1):37-54.

［6］安关峰，张洪彬，刘添俊.旋喷群桩复合地基承载特性的数值分析［J］.岩土力学，2012，33(3):906-912.

［7］张洪彬，安关峰，刘添俊，等.接触面刚度对双排桩支护结构的影响分析［J］.地下空间与工程学报，2012，8(2): 301-306.

［8］VOOTTIPRUEX P，SUKSAWAT T，BERGADO D，et al. Numerical Simulations and Parametric Study of SDCM and DCM Piles under Full Scale Axial and Lateral Loads［J］. Computers and Geotechnics，2011，38(3):318-329.

［9］罗勇.土工问题的颗粒流数值模拟及应用研究［D］.杭州:浙江大学，2007.

［10］周健，陈小亮，王冠英，等.开口管桩沉桩过程试验研究与颗粒流模拟［J］.同济大学学报(自然科学版)，2012，40 (2):173-178.

［11］ELMEKATIA，SHAMYU.APracticalCo-simulation Approach for Multiscale Analysis of Geotechnical Systems［J］.Computers and Geotechnics，2010，37(4):494-503.

［12］Itasca Consulting Group Inc.PFC3D User Manual Version 3.0［R］.Minnesota:Itasca Consulting Group Inc.，2009.

［13］PICHLER B，LACKNER R，MANG H.Back Analysis of Model Parameters in Geotechnical Engineering by Means of Soft Computing［J］.International Journal of Numerical Methods in Engineering，2003，57(11):1943-1978.

［14］GUAN Z，JIANG Y，TANABASHI Y.Rheological Parameter Estimation for the Prediction of Long-term Deformations in Conventional Tunneling［J］.Tunnelling and Underground Space Technology，2009，24(3):250-259.

［15］杜延龄，韩连兵.土工离心模型试验技术［M］.北京:中国水利水电出版社，2010.

［16］王常晶.列车移动荷载作用下地基的动应力及饱和软粘土特性研究［D］.杭州:浙江大学，2006.

Particle discrete element calculation for subsidence of cement mixing pile composite foundation under train cyclic load

GUAN Zhenchang1，LIAO Chonghui1，LUO Zhibin1，GUO Guangzhao2
(1.College of Civil Engineering，Fuzhou University，Fuzhou Fujian 350116，China;2.China Civil Engineering Construction Corporation，Fuzhou Survey and Design Institute Co.，Ltd.，Fuzhou Fujian 350001，China)

Based on the soft ground treatment with cement mixing piles project in Jiangyin railway sub-line，the particle discrete element method(DEM)is used to explore the post-construction subsidence of cement mixing pile composite foundation with train cyclic loads.According to the stress-strain curve from triaxial compression and uniaxial compression lab tests，the micro-parameters of silty and cement soil particles are back-analyzed through an integration methodology，which consists of statistical analysis method and 3D numerical simulation model based on particle discrete element.T hen built the reduced scale DEM model according to the real scale of cement mixing pile composite foundation in Jiangyin Port to calculate the subsidence with train load.From the result，the postconstruction subsidence of cement-mixing pile composite foundation tends to stable around 20 mm after 500～1 000 times of train cyclic load，and this subsidence value is close to the subsidence calculated by existing standard.

Cement mixing pile;T rain cyclic load;Post-construction subsidence;Particle discrete element method

TU472.6

A

10.3969/j.issn.1003-1995.2015.12.28

1003-1995(2015)12-0104-06

(责任审编周彦彦)

2015-06-21;

2015-08-11

教育部留学回国人员科研启动基金(LXKQ0904);国家自然科学基金(51008082)

关振长(1980—)，男，副教授，博士。