离散元-有限差分跨尺度耦合下桩与土体稳定性细观数值分析

高涛 陈云娟 闫甫 李艳龙 刘效智 王乐宁

摘要：为准确分析桩体宏观变形与周围土颗粒细观力学行为，采用离散元-有限差分的跨尺度耦合进行桩和土体接触过程中的稳定性分析，研究桩下沉过程中周围土体的细观变形、应力分布和桩体自身变化情况，通过FLAC3D建立桩和外部土体有限差分网格单元，对桩周围侧土体应用PFC3D离散元建立土颗粒微观结构模型。研究结果表明：离散元与有限差分耦合方案能够模拟桩周围土体细观力学行为;外部区域土体位移場呈包裹式C形对称分布，应力场以竖直应力为主导，离散元土颗粒接触力链网格以树根状向周围递减削弱分布;桩弹性单元体下沉时，应力场和位移场均有分层现象，应力场以压应力为主，但存在局部拉应力即反弹现象。

关键词：多尺度耦合;桩;位移场;应力场;接触力链;数值模拟

中图分类号：TP391.99;TU473.11

文献标志码：B

文章编号：1006-0871（2021）03-0024-07

DOI：10.13340/j.cae.2021.03.005

Abstract：To accurately analyze the macro deformation of pile and the micromechanical behavior of surrounding soil particles， discrete element-finite difference multi-scale coupling is used to study the stability analysis of the pile and soil in the contact process， and then the micro-deformation， stress distribution and changes of the pile itself during the pile sinking process is studied. The finite difference mesh of the pile and the external soil is established by FLAC3D. The soil particle microstructure model of the soil around the pile is established by PFC3D discrete element. The results show that the coupling scheme of discrete element and finite difference can be used to simulate the micromechanical behavior of the soil around the pile. The soil displacement field in the outer region is symmetrically distributed in a "wrapped C" shape， and the stress field is dominated by vertical stress. The soil particle contact force chain mesh is gradually weakened and distributed in the shape of tree roots. While the pile elastic element body sinks， the stress field and displacement field are both stratified， and the stress field is mainly compressive， but there is local tensile stress， that is rebound phenomenon.

Key words：multi-scale coupling;pile;displacement field;stress field;contact force chain;numerical simulation

0 前 言

我国很多地域的地层土质状况不稳定，由于雜填土成分复杂且持力层较深等因素，解决建筑施工地基层面问题常利用桩的径长和自身强度高、稳定性好等优点，但施工过程主要依靠大量实践衡量桩体的承载能力，桩下沉过程中底部土体的稳定性很难监测，桩土之间的载荷传递机理也不够完善。[1-3]

为进一步了解桩土受力机理，已对其相互作用关系展开许多室内试验，但大多数室内试验采取比例控制模型，成本花费高、地质情况单一，造成试验数据与现场监测差异较大。因此，国内外许多学者采用数值模拟技术进行此类研究，并取得一些进展。[4-6]周文龙等[7]应用有限差分数值软件FLAC3D模拟桩体在复合地基情况下的自身沉降量和地表土体沉降，得到桩土在不同工况条件下应力比，为实际工程确定最优垫层厚度。邱瑞成等[8]对单桩进行参数敏感分析，改变数值模型中黏聚力和内摩擦角参数，结果认为黏聚力对计算结果影响较显著，而内摩擦角参数对结果影响基本可忽略。季璇等[9]应用有限元MIDAS/GTS分析双排桩支护过程中的变形情况，分析其水平位移和桩身弯矩，选取单排桩支护最佳优化方案。李新伟等[10]从桩对土层影响的角度出发，结合桩周围土体应力场和位移场，分析施工过程中周围土体的变形。刘文白等[11]应用流离散元PFC2D模拟桩上拔过程中的土体承载力，研究桩周围土颗粒细观力学变化，分析桩上拔过程中土体结构单元的变形情况，并与物理模拟试验结果对比，认为颗粒流离散元软件对工程模拟具有可靠性。

采用单一数值软件对桩和桩周围土体分析，桩和土体均为同一种计算原理，如有限元法、非连续性变形法、离散元法等，但实际上两者结构具有很大差异性，桩自身强度大、变形小，而土体性质复杂且呈颗粒离散状态，因此从不同角度研究桩和土体性质更为重要[12-13]。本文将FLAC3D和PFC3D进行耦合，将土体松散不稳定等特性定义为离散元体，桩体和周围宏观土层视为弹塑性材料，通过数值模拟分析桩体宏观变形与周围土颗粒细观力学行为。

1 FLAC3D-PFC3D耦合理论基础

1.1 耦合原理

在FLAC3D和PFC3D数值软件相互耦合分析中，FLAC3D从宏观上模拟连续区域内介质的力学行为，采用显示拉格朗日算法和混合离散分区原理，运用动态运动方程进行求解计算，拉格朗日算法将大变形问题转化为小变形本构关系，通过离散集成方法求解应力增量和不平衡力。PFC3D用于模拟细观离散元颗粒的微观力学行为，计算原理为颗粒间相互关系的离散元方法，颗粒间力的传递与颗粒变形通过二者接触刚度（法向刚度和切向刚度）进行计算，接触刚度由强度参数（摩擦系数、黏结和抗拉强度）控制[14]。两者相互耦合发生在域的接触边界上，不同域之间的相互关系主要借助Socket O/I接口进行数据传输与交换，FLAC3D-PFC3D耦合计算原理见图1。

1.2 耦合方案分析

为探讨FLAC3D和PFC3D的相互耦合作用，模拟连续区域内FLAC3D与离散介质PFC3D中球、簇的相互作用，在边界区域创建独立的耦合体域，耦合体域使球、簇等结构单元与实体单元相互作用，耦合体边界采用WALL单元作为相互交换媒介。WALL单元是由点构成三角形面区域的组成体，三角形区域位置可由时间函数确定，以此作为耦合边界上的媒介体。因此，在耦合分析中，WALL单元必须与实体单元表面协调一致。

初步确立耦合边界介质后，通过媒介体传递和接受单元体力、速度等信息，形成区域通道内的耦合逻辑（通过墙面接触力、力矩确定三角面顶点处的等效力系统）。确定耦合逻辑后定义相应的耦合工作，即取接触力和力矩与壁面，确定等效的力系统在面上的顶点，这些力与刚度发送信息至媒介体区域，相关的网格节点和单元体接受媒介体传递的信息。数值模拟耦合中力的传递示意见图2。

若没有附加约束条件，式（14）和（15）构成的欠正定方程组可能存在无数组解集。控制方程并没有提供明确的约束性质，因此欠正定方程组需得到一个特解，同时需要提供约束实体单元或结构单元的额外约束条件，通过该额外约束条件确立欠正定方程组的等效系统，等效系统需要的额外约束条件为

当点OP外推三角形区域的节点和结构单元触发更新时，其网格点与节点能够直接添加刚度属性，即分别对FLAC3D于PFC3D设置实体单元的相关参数，通过三角形区域这一介质实现传递力与速度信息。

2 连续-离散耦合的桩与土体数值模拟

2.1 耦合模型建立

以FLAC3D软件作为耦合计算方案的平台，通过设置软件Load PFC打开耦合开关系统，建立的数值耦合模型主要包括桩和土体2个模块，桩结构体采用FLAC3D实体单元，下部土体分为2个模块。外部土体采用FLAC3D实体单元作为“基床”，该实体单元利用有限差分数值原理计算周围大范围土体，处理大变形，较真实反映材料实际的动态行为，有效模拟外部土体随时间、空间变化产生的系统相关问题。与桩周围接触的土体采用PFC3D离散元结构体，离散元模型不仅不受变形限制，而且处理非连续性介质力学问题更精确，符合实际土颗粒物理力学性质。外部土体模型尺寸为15 m×10 m×8 m，土体单元采用FLAC3D中的Mohr-Coulomb计算模型，桩周围局部土体区域范围为5 m×5 m×5 m的 PFC3D颗粒（BALL），颗粒直径分布区间为0.2～4.0 mm。为使桩的数值模型与施工现场更好吻合，桩局部直径为0.6 m、高度为5 m，计算模型采用FLAC3D弹性体单元，整体系统数值计算模型见图3。从多个不同单元体模型体现不同的物理性质关系，可有效模拟土体结构分离、变形等非连续特性，优化分析土体与桩的细观力学行为。

2.2 模型参数标定

选取济南开源路地区粉质黏土力学参数，按照现场试验数据设置模型参数，在桩周围局部土体单元设置线性接触模型，对2种不同土体之间连接处设置WALL单元，作为FLAC3D有限差分网格与PFC3D离散元颗粒耦合的边界。对于离散元计算模型的细观参数，可通过一系列三轴数值试验绘制不同围压强度下的莫尔强度包络曲线，并进行多次标定计算，得到离散元模型细观参数，桩周围土体离散元模型材料参数见表1。FLAC3D区域内的材料参数根据现场测定得到，桩周围土体的材料参数见表2，桩体的变形模量为4.86 MPa，泊松比为0.15。

3 数值模拟结果分析

采用离散元（PFC3D）与有限差分（FLAC3D）耦合模拟桩与土体的接触变形关系，主要通过模拟桩接触土体颗粒，分析桩身和周围土颗粒的应力和变形情况。2种不同的单元体受力形态变化见图4，可验证PFC3D与FLAC3D相互耦合的可行性。PFC3D颗粒流单元的应力边界传递至FLAC3D实体网格中，同一工况加载状态下其应力云扩散性效果较为一致，PFC3D和FLAC3D的应力云梯度基本相同，均无网格单元体或球颗粒单元因应力传递问题导致应力点集中现象，二者在边界处具备耦合相互传递性能，说明边界连接体单元WALL对颗粒程序和有限差分网格的融接性好，同时也证明Socket O/I接口进行计算数据传输与交换的可行性。

3.1 桩周围外部有限差分模型土体变形特性

数值模拟加载后桩底周围土体的水平位移场见图5a。在桩加载下沉过程中，外部土体受到桩周围土体传递应力产生一定的位移，土体接触面附近水平方向产生4个较为显著的聚集性位移场，最大位移量为3.3 mm。产生该现象主要原因是桩在下沉过程冲击土体，造成桩周围土体被迅速密实挤压，产生冲击力侧压现象，使水平位移云图呈包裹式C形状分布;坑底侧土体在桩下沉挤压作用下呈两侧包裹圆位移，但从纵向切片图来看，底部较深处的一侧土体位移并非连续的，该情况与现场土体位移监测较为符合。桩体下沉受载非一点集中力，在桩的竖直方向集中力与土体接触过程中，产生中心线偏离引起桩身弯矩力效应，使下部土体位移不是等距变化。

数值模拟加载后桩底周围土体竖直方向的位移场见图5b。位移场基本呈现法向散状分布，土体沉降量自中心向四周减缓，桩底中心处土体压缩位移量约为四周的2倍左右，局部最大位移达19 mm左右，但从竖直方向位移切片图看，局部最大位移区域范围较小，说明越靠近樁底端部，土体压缩位移越大，竖直位移分布云图也可以更好地表明对土体建模划分的科学性，将土体压缩量较大区域用离散元颗粒表示更能贴合实际施工状态。

桩周围外部土体应力场变化云图见图6。

在外部上表层土体应力场中，水平应力场拉应力和压应力对称性分布，拉应力区域面积远大于压应力区域，整体以拉应力产生的土体变形为主导，但最大压应力数值比拉应力约大0.8 kPa;竖直应力场中表面土体均为拉应力且无压应力存在，主要原因是在桩竖向受载中，底部土体受压引起上部土体产生被动受拉效应，而在桩底部的土体中，水平应力场呈现“两压一拉”现象，且压应力区域远远超过拉应力区域范围。这说明同等时间段内桩底中心处土体较周围土体先压缩固结，周围土层固结量小，但中心竖向载荷不断增加，同一深度应力场发生改变。

3.2 桩周围离散元模型土体变形特性

土颗粒的水平位移场和竖直位移场见图7。由此可以看出，位移量之间的传递具有较好的衔接性，土颗粒位移基本与外侧土体位移趋势相同。离散元土颗粒水平位移量集中在土体表面两侧，最大位移量局部相差2 mm左右，竖直方向位移呈现中间大两边小的形态，颗粒位移递减量约35 mm。

为进一步研究土体位移过程中的细观变化，提取离散元模型的土体速度矢量云图，见图8。

从2个方向的速度矢量场可知，水平方向的速度矢量运动轨迹基本与竖直位移场相同，竖直方向的颗粒运动轨迹与其水平位移场相同。速度矢量场这一现场可解释在上部桩载作用下，桩载压缩土体导致水平方向的土体颗粒向中间靠拢，而中心处土体在桩直接压缩下先固结，因此周围土颗粒向下位移产生拉应力。这种运动轨迹符合竖直方向位移场中间大两边小的规律。

3.3 离散元模型土体细观变形特性

通过离散元颗粒模型构建桩周围土体，能够更好地从细观角度分析颗粒间受力形态的变化。应用FLAC3D与PFC3D耦合技术可得到接触模型颗粒细观接触力，见图9。颗粒与颗粒间通过接触力链传递应力和位移，接触力链网格约35 000个，其接触力链相互连接形成交叉网格呈树根状分布，在桩体与颗粒接触轴向处，这些接触力链网格分布密集并向周围散射。

桩底部竖向载荷的传递，使桩底侧颗粒出现强度高但数量少的接触力链网格。随着深度的增加，高强接触力链网格转化为数目更多的低强度力链网络。在高强度接触力链网格周围，水平方向的接触力基本以拉应力为主，但压应力出现最高峰值82.41 N;竖直方向的接触拉应力和压应力数量相同，但竖向接触力链应力值远大于水平方向，接触应力最大值为1.97 kPa。这说明桩体下沉中竖向应力起主导作用。

接触力链网格分布与应力传递形态可揭示桩体下沉中土体的细观特征，桩顶载荷通过桩身弹性体传递至桩底侧土颗粒，传递的应力以竖向应力为主导作用，形成树根状接触力链网格，其中心处土体固结速度大于周围土体压缩量，使周围产生水平方向的接触力链网格。这些水平方向接触网格相对强度较低，从而对竖向接触力链网格起到辅助支撑作用。若这些土颗粒无法对竖向接触力链网格起约束作用，则会导致桩体下沉过程中出现力链坍塌错落现象，甚至引起桩体宏观坍塌破坏;若接触力链对竖向力链约束程度过强，则会导致水平应力集中，使土体表面出现隆起，导致原有土层破坏。

3.4 桩有限差分网格单元变形分析

应用FLAC3D切片功能，将桩有限差分网格单元沿纵向作切片处理，得到桩单元不同方向的位移场和应力场，见图10。桩身上半部分水平位移场与应力场保持相对均匀，主要原因是桩单元模型采用弹性体模型，该模型自身强度高、变形较为一致;下半部分桩体位移场与应力场不同，主要原因是土体产生的水平方向的束缚应力和桩自身摩擦阻力，使桩下半侧水平位移场呈相对挤压的位移，竖直位移场具有明显的层次阶段变化，且越靠近桩底端部，位移场分布区域越聚集于角落。

因此，针对桩身各个部位网格单元进行位移变化监测，结果见图11。由此发现，越靠近桩底端部，水平位移越小，桩的竖直方向位移变化越大，说明弹性体桩单元符合胡克变形准则，弹性体变形基本以轴向为主导作用，水平位移量变化主要是由于土颗粒间的接触力链约束引起的。

图10b说明桩身应力场变化主要集中于桩底侧单元区域，水平应力场呈轴对称分布形态。水平方向土颗粒间受接触力链约束作用，竖直方向的应力场上部为拉应力作用，但拉应力数值较下部压应力小（压应力最大值为238 kPa，是拉应力近20倍），桩身上半部的拉应力主要是自身在压缩受力过程中发生的应力回弹现象，因此上部出现较小的均匀拉应力分布。

4 结 论

通过离散元-有限差分尺度耦合对桩与土体的稳定性进行细观数值分析，得到以下结论。

（1）应用离散元和有限差分原理进行耦合计算，可实现PFC3D与FLAC3D数值模型耦合计算，其数值模型的位移场、应力场等变形传递在不同属性单元中具有良好的融合性，为跨尺度耦合数值模拟提供参考。

（2）桩周围外部土体位移场呈包裹式C形对称性分布，竖直应力场起主导作用，且约为水平方向应力7倍;内部离散元土颗粒位移场与外部土体基本一致，土颗粒的速度矢量场可解释土颗粒在固结作用下的运动变化，其细观接触力链网格分布呈树根状并向周围递减，接触应力密集处为桩底部土体附近。

（3）对桩弹性体单元作切片处理可知，桩下沉应力场和位移场有分层现象，应力表现以压应力为主但有较小的拉应力反弹现象，水平位移随桩身增加而增加，而竖直位移则相反。

参考文献：

[1] 刘汉龙， 赵明华. 地基处理研究进展[J]. 土木工程学报， 2016， 49（1）：96-115.

[2] STAHL M， KONIETZKY H. Discrete element simulation of ballast and gravel under special consideration of grain-shape， grain-size and relative density[J]. Granular Matter， 2011， 13（4）：417-428. DOI：10.1007/s10035-010-0239-y.

[3] 朱怡， 魏祥， 赵登坚. 桩体贯入大变形计算分析[J]. 建筑结构， 2018， 48（S2）：838-841. DOI：10.19701/j.jzjg.2018.S2.171.

[4] 寇海磊， 张明义， 张吉坤. 层状粘性土及砂土地基中静力压桩连续贯入的数值模拟[J]. 工程力学， 2012， 29（12）：175-181.

[5] HENKE S. Influence of pile installation on adjacent structures[J]. International Journal for Numerical and Analytical Methods in Geomechanics， 2010， 34（11）：1191-1210. DOI：10.1002/nag.859.

[6] GAO G， MEGUID M A. Effect of particle shape on the response of geogrid-reinforced systems：insights from 3D discrete element analysis[J]. Geotextiles and Geomembranes， 2018， 46（6）：685-698. DOI：10.1016/j.geotexmem.2018.07.001.

[7] 周文龙， 刘玉柱， 王小鹏. 应用FLAC3D分析确定碎石桩复合地基的褥垫层厚度[J]. 施工技术， 2014， 43（S1）：76-78.

[8] 邱瑞成， 艾健森. 基于FLAC3D的单桩静载模拟桩土接触面参数敏感性研究[J]. 路基工程， 2019（2）：164-169. DOI：10.13379/j.issn.1003-8825.2019.02.34.

[9] 季璇， 高全臣， 杨卓， 等. H型双排桩支护结构三维数值模拟分析与工程应用[J]. 施工技术， 2016， 45（7）：70-73. DOI：10.7672/sgjs2016070070.

[10] 李新偉， 曾昌. 桩基础施工对地层影响的三维数值模拟[J]. 四川建筑， 2015， 35（1）：95-97. DOI：10.3969/j.issn.1007-8983.2015.01.036.

[11] 刘文白， 周健. 上拔荷载作用下桩的颗粒流数值模拟[J]. 岩土工程学报， 2004， 26（4）：516-521. DOI：10.3321/j.issn：1000-4548.2004.04.018.

[12] 谭鑫， 冯龙健， 赵明华， 等. 刚性基础下筋箍碎石桩复合地基桩土应力比计算模型[J]. 中国公路学报， 2019， 32（9）：42-50. DOI：10.19721/j.cnki.1001-7372.2019.09.004.

[13] 傅金阳， 谢佳伟， 房雅楠， 等. EPB盾构开挖面稳定性的PFC-FLAC耦合分析[J]. 华中科技大学学报（自然科学版）， 2019， 47（5）：116-121. DOI：10.13245/j.hust.190522.

[14] 石崇， 张强， 王盛年.颗粒流（PFC5.0）数值模拟技术及应用[M]. 北京：中国建筑工业出版社， 2018.

（编辑 武晓英）