真实复杂地层大直径钢顶管三维数值模拟研究

张艳林， 曾天成， 甘 甜， 叶 阳， 王新顺

(1.中国水利水电第七工程局有限公司， 成都 610081； 2.中建交通建设集团有限公司， 北京 100142； 3.武汉大学土木建筑工程学院， 武汉 430072)

随着城市建设的发展，地面建(构)筑物越来越密集，地上交通越来越繁忙，在城市建成区进行市政管线建设时采用明挖法施工的可能性逐渐减少，取而代之的是采用对地面交通影响较小的非开挖法施工，其中顶管法是一种典型的非开挖市政管线施工方法，应用逐步增多。因此，研究顶管法施工对周围土体的扰动及地表沉降的影响已受到学术界和工程界的广泛关注。

目前中外学者通过现场测试[1-3]、理论分析[4-8]、数值模拟[9-14]等方法对顶管施工造成的地表沉降进行了研究。由于顶管施工中的管土相互作用是一个复杂的动态三维过程，因此仅开展单纯的理论研究难以解决具体的工程问题，而数值模拟并结合现场的实测结果进行分析，更有助于理解顶管施工对周围土体扰动的机理，也可以为工程施工中遇到的问题提供合理的解决方案。实际的顶管工程中往往穿越的是错综复杂的地层环境，而目前的顶管工程数值模拟往往在简化后的均匀地层模型中进行，与实际的工程地质环境不相符。参照盾构隧道模拟研究成果，采用符合实际的复杂地层模型开展数值模拟，更加符合工程实际[15-17]。

所谓复杂地层，不仅指顶管工程横断面上的地层分布的复杂性(成层分布、不均匀分布)，更重要的是指顶管工程纵断面上地层的不均匀分布。现有的数值模拟对于复杂地层的处理一般均是横断面层状均匀分布，而纵断面则分段层状均匀分布，导致各分段边界呈台阶状，与工程实际的差距仍然较大。为此，结合某市北线引水工程隐患整改工程，提出真实复杂地层模型建模方法，用于模拟顶管施工对周围地层扰动及地面沉降的影响，作为对比，同时建立常规简化均匀地层模型并进行数值模拟，通过数值模拟结果与现场实测结果的比较，说明采用提出的真实复杂地层模型的合理性和必要性。相关研究结果可供工程实际及数值模拟时参考。

1 数值建模及顶管施工模拟方法

1.1 工程概况

某市北线引水工程隐患整改工程E段采用顶管施工，采用NPD3 400泥水平衡顶管机顶进，顶管为DN3 400 mm×32 mm钢管。顶管外径为3.4 m，壁厚为0.032 m，顶管机头外径为3.46 m。选择的模拟区段位于工作井和接收井中间，且距工作井和接收井距离较远，可忽略工作井和接收井施工造成的土体扰动，因此地层变形及地表沉降量可认为均由顶管施工引起。

1.2 真实复杂地层三维数值建模方法

基于有限差分软件FLAC3D进行数值模拟。虽然FLAC3D6.0具备强大的运算能力和较好的后处理功能，但其前处理能力较弱，在建立复杂地层的三维数值模型时存在一定困难。因此提出采用有限元软件ABAQUS来建立复杂地层的三维数值模型。首先根据实际工程地质勘测确定的地层纵断面图，对数值模型进行空间划分，然后在ABAQUS中将模型剖分为四面体四节点网格单元。将地层信息和剖分网格数据导入FLAC3D6.0，然后进行边界条件设置、物理力学参数赋值及地应力平衡。最终得到真实物理力学参数及地层分布的三维数值模型。

为了对比分析真实复杂地层的准确性，在模型纵向采用各土层体积等量代换的方法对土层进行简化处理，生成纵向地层均匀成层分布的简化地层模型，并导入FLAC3D6.0进行顶管顶进施工模拟。

由于顶管工程的钢管每节长度为6 m，因此数值模型的纵向模拟范围取6 m的整数倍，取60 m，数值建模时将开挖区和注浆层单元均分为10段，每段长6 m。相应地，数值模型在横向模拟范围也取60 m，顶管两侧模拟范围约为8倍的顶管直径，可以认为模型两侧边界条件对模拟结果没有影响，故在两侧边界上均匀施加法向(水平)约束。垂直方向，模型上表面取为地面，顶管轴线埋深约12 m，取模型下部模拟范围为20 m，超过顶管直径的5倍，底部边界的影响可忽略不计，故模型上表面为自由边界，下表面施加法向(垂直)和切向(水平两向)约束。此外，在模型轴向，也在相对两表面施加法向(水平)约束。

根据实际工程地质勘测资料，数值模拟段地层由上至下分别为素填土、粉质黏土、粉砂岩残积土、全风化风砂岩，各地层沿顶管纵向分布不均，尤其是素填土层和粉质黏土层层厚差距大。在模型纵向模拟范围内，素填土层厚2.2～11.4 m，粉质黏土层厚3.2～16.8 m，粉砂岩残积土层厚4.1～ 10.1 m，全风化风砂岩层厚7.4～10.0 m；采用各土层体积等量代换的方法对土层进行简化处理后，在模型纵向模拟范围内各地层等效厚度分别为素填土层厚7.0 m，粉质黏土层厚8.0 m，粉砂岩残积土层厚9.1 m，全风化风砂岩层厚7.9 m。

建立的真实复杂地层数值模型如图1(a)所示，相应简化均匀地层数值模型如图1(b)所示。图1中，①～表示不同的监测断面。两模型的地层土体均剖分为四面体单元，顶管管体剖分为壳单元，注浆层厚度为0.03 m，也剖分为四面体单元。单元剖分后，真实复杂地层数值模型的单元总数为 696 633，结点总数为117 133；简化模型的单元数总为648 311，结点总数为109 305。

图1 计算模型网格Fig.1 Computational model grid

本构模型采用各向同性弹塑性模型，屈服准则为Mohr-Coulomb准则。具体物理力学参数取值根据工程地质勘察现场试验确定，如表1所示。钢顶管采用线弹性本构模型，材料重度为77 kN/m3，弹性模量为200 GPa，泊松比为0.3。注浆层的作用，主要是防止管周地层因为管土之间空隙的存在而发生过大变形，数值模拟时主要考虑其防止地层变形的作用，以空单元和作用在地层上的反力来模拟注浆层的作用[18]，该反力大小与地层压力、开挖后的地层应力释放系数等有关。

表1 各地层参数

数值模拟过程中，沿着顶管顶进方向每6 m设置一个监测断面，共设置11个监测断面，相应地地面上有11条监测线，如图1所示。在每条监测线上设置23个监测点，中间点位于顶管轴线正上方，两侧监测点对称分布于顶管轴线两侧，相邻监测点的间距为2 m。

采用的顶管机是泥水平衡顶管机，开挖面处的泥水压力与地下水压力平衡，并且钢顶管与周围土体之间充满膨润土泥浆，可以有效控制施工对地下水位的扰动。触变泥浆的使用可以显著降低管土摩擦力，从而消除摩擦力对土体造成的扰动。因此，在充分考虑泥水平衡顶管机的工作原理的情况下，对数值模型建模和计算作如下假设：①在模拟顶管顶进过程时不考虑土体变形的时间效应；②作用于掌子面上土体的压力为圆形均布荷载，其值为实测的顶管机头泥水舱压力；③不考虑顶管施工过程中的管-土摩擦力；④不考虑地下水的影响。

1.3 施工过程模拟

顶管工程施工过程的模拟步骤如下。

(1)初始地应力平衡，得到初始地应力场，并将初始位移归零。

(2)将第一段的开挖区单元和注浆层单元挖去，进行1次循环，然后得到注浆层最外层的节点上的不平衡力，将这些节点的不平衡力乘以一个小于0的系数，并作用于对应的节点上，以模拟注浆压力对管周土体的作用[18]。该系数最理想的取值为-1，取-1时表示注浆压力正好平衡了该处的不平衡力，对应无土体损失的情况。然而，实际施工中无法实现这种理想情况，通过反复试算，系数取 -0.85 时，模拟的结果与实际工程的监测值拟合较好。添加 shell弹性结构单元以模拟顶管，同时，根据项目实测数据，在掌子面处加上0.12 MPa均布荷载模拟泥舱压力。

(3)以6 m为一个开挖步，依次开挖至全部施工完毕。施工过程模拟示意图如图2所示。简化地层模型的施工过程与真实复杂地层模型完全一致。

图2 施工过程模拟示意图Fig.2 Schematic diagram of construction process simulation

2模拟结果与分析

2.1 地表沉降结果与分析

图3为顶管顶进54 m时各典型监测断面监测的地表沉降曲线。从图3可以看出，两种模型中各典型断面的地表沉降最大值均出现在顶管轴线正上方，并沿管道横向两侧逐渐减小，沉降曲线类似于倒正态分布，与Peck得出的沉降变化规律一致[19]。结合图1(a)可知，真实复杂地层模型中顶管上覆土层中素填土厚度排序依次为：第4监测断面<第6监测断面<第8监测断面。对于真实复杂地层模型结果，各典型监测断面沉降值的排序为：第4监测断面<第6监测断面<第8监测断面，且差异非常明显，与素填土厚度排序一致；在简化模型中，3组监测值之间差异很小，但由于边界效应的影响，第4监测断面的地表沉降比第6监测断面和第8监测断面的地表沉降略小。将地表沉降小于0.3 mm的区域视为无影响区域，则两种模型的地表横向受影响范围为-18～18 m即(-5.3～5.3)D(D为顶管直径)，表明本文模型的横向模拟范围取值是合理的。

图3 顶管顶进54 m时的地表沉降曲线Fig.3 Surface settlement curve when pipe jacking is pushed into 54 m

图4为工程顶管施工现场地表沉降实测结果与数值模拟对应地表沉降结果的比较。从图4可以看出，在第3监测断面处，数值模拟结果均大于现场实测结果，其中真实复杂地层模型的模拟结果比较接近实测结果，而简单地层模型的模拟结果则与实测结果相差较大，其最大沉降比实测最大沉降大约60%。在第7监测断面处，真实复杂地层模型的模拟结果比实测值稍大，其最大沉降与实测值的相对误差不超过10%，而简化地层模型的模拟结果明显比实测结果小，其最大沉降与实测最大沉降值的相对误差约为30%，相比较而言，真实复杂地层模型模拟获得的地表沉降曲线更接近于实测地表沉降曲线。此外，从图4还可看出，无论第3监测断面处还是第7监测断面处，现场实测地表沉降值都比真实复杂地层模型模拟结果稍小，原因是该监测区段的地表为水泥路面，而水泥路面的强度比一般土体高，抵抗变形的能力更强，导致实测地表沉降比数值模拟结果小。由此可见，采用所提出的真实复杂地层模型模拟的结果与地表沉降实测结果吻合较好，而采用简化地层模型模拟的结果则与地表沉降实测结果相差较大。

图5为真实复杂地层模型中顶管机刀盘(开挖掌子面)与监测断面不同距离时的地表沉降曲线，以顶进方向为前方，L=6 m表示顶管机刀盘(开挖掌子面)在监测断面前6 m(已通过监测断面)；L=-6 m表示顶管机刀盘(开挖掌子面)在监测断面后6 m(尚未通过监测断面)。由图5可见，当顶管机机头由L=-6～6 m段顶进时，对于监测断面地表沉降的影响均较大，因此选取L为-6～6 m顶进段进行分析。当顶管机刀盘距第4监测断面-6 m时，第4监测断面顶管轴线处地表竖向位移约为-0.5 mm；当刀盘通过第4监测断面且距离为6 m时，第4监测断面顶管轴线处地表竖向位移约为-3.0 mm。顶管机机头从L=-6～6 m的过程中产生的顶管轴线处地表沉降增量为2.5 mm，约占施工全过程中第4监测断面轴线处最终沉降的66.7%。从图5(b)可见，当顶管机刀盘距第8监测断面-6 m时，第8监测断面顶管轴线处地表竖向位移为-1.0 mm；当刀盘通过第8监测断面且距离为6 m时，第8监测断面顶管轴线处地表竖向位移为-7.0 mm。顶管机机头从L=-6～6 m的过程中产生的顶管轴线处地表沉降增量为6 mm，约占施工全过程中第8监测断面轴线处最终沉降的75%。由图5可知，顶管机机头穿越第8监测断面产生的沉降量增量及其占总沉降的比例均大于机头穿越第4监测断面时的相应值，说明顶管穿越上覆软弱土体越厚的区段，顶管施工过程中的地表沉降反应越剧烈。由图5可知，在数值模拟过程中地表都没有产生明显的隆起，原因是工程的泥舱压力控制在较小值，若泥舱压力较大，则顶管机前方土体地表应当会产生一定程度的隆起。从现场监测数据来看，整个施工过程中，地表土体很少出现隆起，且当隆起出现时，数值也较小，监测结果与真实复杂地层模型的模拟结果基本一致。

图4 数值模拟结果与现场实测值对比Fig.4 Comparison between numerical simulation results and field measured values

图5 顶管机刀盘与监测断面不同距离时的地表沉降曲线Fig.5 Surface subsidence curve at different distances between cutter head and monitoring section

2.2 管周土体扰动结果与分析

图6为真实复杂地层模型和简化地层模型中顶管顶进48 m时，过顶管轴线的垂直截面z方向(垂直方向)位移云图。由图6可知，由于土体开挖卸载、地层损失等因素造成顶管上方土体产生明显的下沉，顶管下部土体产生一定程度的抬升，且上部土体下沉变形影响范围远大于下部土体抬升变形的影响范围；过轴线垂直截面上越靠近顶管的土体变形越大，而远离顶管的土体变形逐渐减小；过轴线截面的下沉变形最大，偏离轴线截面的地层变形逐渐减小，符合沉降槽的基本特征；真实复杂地层模型的地层z向变形沿顶管纵向呈不均匀分布，局部变形明显较大，而简化地层模型的地层z向变形沿顶管纵向的分布则相对比较均匀，说明纵向地层的等效简化有可能导致分析结果变小，使工程偏于不安全。

图7为两种模型中顶管顶进54 m时的过顶管轴线垂直截面z方向位移云图。从图7可以看出，与图6相比，顶管轴线截面z方向变形特征更加明显，便于二者的比较。图7中简化地层模型的土体竖向位移沿顶进方向分布较为均匀，而真实复杂地层模型的土体竖向位移沿顶进方向分布极不均匀。由表1可知，土层的强度大小和变形模量大小排序均为：素填土<粉质黏土<粉砂岩残积<全风化粉砂岩。对比图7(a)、图8(a)可知，真实复杂地层模型中土体沉降和抬升变形大小与地层的纵向分布密切相关，沿顶管顶进方向，地层沉降最大处位于顶管上方素填土最厚处，而顶管下部粉质黏土最薄处土体抬升变形最小。

图6 顶管顶进48 m时的土体z方向位移云图Fig.6 Cloud map of soil displacement in z direction when jacking pipe is jacked in 48 m

图7 顶进54 m时过顶管轴线纵断面的土体z方向 位移云图Fig.7 Cloud map of z-direction displacement of soil passing through longitudinal section of jacking pipe axis when jacking 54 m

图8 土层分布纵断面Fig.8 Vertical profile of soil layer distribution

图9、图10为两种模型中顶进54 m时各典型监测断面的x方向(水平方向)位移云图。从图9、图10可以看出，由于开挖卸荷，导致顶管两侧的土体均向顶管中心移动，且顶管上覆土体沿x方向的位移量明显大于顶管下部土体，整个位移云图呈“蝴蝶状”。简化模型中，3个典型监测断面的x方向位移云图基本相同，而真实复杂地层模型中，顶管下部土体的x方向位移排序依次为第4监测断面>第6监测断面>第8监测断面，且差异明显。结合图11可知，顶管下部软弱土层厚度越大，顶管施工过程中造成的下部土体x方向位移越大。真实复杂地层模型的第4监测断面和第6监测断面，上覆土体部分的x方向位移云图无明显差异，而第8监测断面的①区和②区产生了远离顶管中心的x方向位移。结合图11可知，与第4监测断面和第6监测断面不同，第8监测断面的上覆土体及部分开挖区均为素填土，由表1的数据可知，素填土的抗剪强度低于粉质黏土，因此，在③区和④区土体的剪切滑移和顶管注浆压力的共同作用下，①区和②区的素填土产生了远离顶管中心的x方向位移。对比第4监测断面和第6监测断面位移云图可知，该现象仅在顶管穿越软硬土层交界处时发生。

3 结论

依托某市北线引水工程隐患整改顶管工程项目，提出了真实复杂地层三维数值建模方法，开展了钢顶管施工数值模拟，分析了顶管施工过程中管周土体扰动变形特征和地表沉降规律，并与简化地层模型模拟结果及现场地表沉降实测结果进行了比较，得出如下主要结论。

图9 真实复杂地层模型中顶进54 m时各典型监测断面的x方向位移云图Fig.9 X-direction displacement cloud diagram of each typical monitoring section when jacking 54 m in the real and complex strata

图10 简化地层模型中顶进54 m时各典型监测断面的x方向位移云图Fig.10 X-direction displacement cloud diagram of each typical monitoring section when jacking 54 m in the simplified and uniform strata

图11 真实复杂地层模型中各典型监测断面的土层分布Fig.11 Distribution of soil layers in each typical monitoring section in the real and complex strata

(1)提出的真实复杂地层三维数值建模方法能够准确反映真实地层纵向非均匀空间分布规律。与常规简化均匀地层模拟结果相比，真实复杂地层模拟得到的地表沉降纵向分布规律与现场监测结果吻合较好，能够揭示管周土体扰动沿顶进方向的真实变化规律。

(2)顶管施工会引起其上方地层的沉降变形和下方地层的抬升变形，但其上方地层的沉降变形的影响范围明显大于下方地层抬升变形的影响范围。地层的变形与地层土体性质和分布有关，土体强度和变形模量越低，其变形越大；顶管上部或下部软弱土层越厚，其变形越大。

(3)顶管上覆土体横向位移基本朝向顶管轴线，其变形受土体空间分布的影响较小，但顶管下部土体的横向变形受土体空间变化影响较大，随着顶管上覆软弱土层厚度的增加而增大。在顶管穿越软硬两种土体(上软下硬)的交界处时，顶管上部两侧的部分土体会产生一定程度的背向顶管轴线的横向位移。

(4)数值模拟结果显示，地表沉降的横向受影响范围为±5.3D(D为顶管直径)，且受土层分布影响较小；但地表沉降受地层分布影响较大，顶管上覆地层中软弱土层越厚，地表沉降越大，顶管施工过程中的地表沉降反应越剧烈。因此，顶管施工时应当着重关注软弱土层较厚的区段的施工。