估算竖向曲线顶管顶进力的整体位移控制有限单元模型建立的难点在于每个管节经过顶管线路转折处会形成1个平动和1个转动，在建立模型时，除在起始管节施加1个管节厚度的位移外，还需要对经过每个转折处的管节施加1个管节厚度的平动位移和1个转折角度的转动位移。对模型要求更加苛刻，网格划分更加精细，但是计算量降低，计算过程简化。
1 工程概况
郑州市轨道交通4号线的会展中心站，为4号线与1号线的换乘车站，其中该站的1号出入口下穿商务内环路，采用竖向曲线矩形顶管法施工，长度91.6 m。从始发井到接受井依次为：坡度-3%、长度27.0 m；坡度-2%、长度9.0 m；坡度-1%、长度9.0 m；坡度0、长度46.6 m。轴向坡度分布如图1所示。顶管隧道围岩主要为素填土、砂质黏土以及砂土。顶管穿越土层主要为砂质黏土，局部区段为砂性土。地下水位位于地表以下11 m，对顶管的影响可以忽略。顶管施工采用土压力平衡法，采用高黏度膨润土浆液软化减阻。施工选用带顶进力、端头土压力和触变泥浆注浆压力监测功能的矩形顶管机，记录施工过程中顶进力、端头土压力和触变泥浆注浆压力的变化。该过街通道平面布置图如图2所示。
2 有限元模型设计
2.1 几何尺寸与单元划分
有限元模型建立的目的是通过在顶管起始环形端面施加位移荷载计算不同顶进距离下顶管顶进位置的顶进力，实现对顶管顶进力的估算。首先利用Abaqus有限元软件建立曲线顶管的整体有限元模型（图3），模型的实际长度为91.5 m，按照实际顶管的轴向坡度分布进行顶管管节排布，顶管模型如图4所示。考虑边界效应，依据圣维南原理，设模型土体的宽为39.1 m，高为26.0 m，长为91.5 m。模型单元为8结点六面体单元。顶管转折处、机头顶进位置处的网格适当加密。网格数量既要满足计算精度，又要满足计算效率，确保每个模型的计算时间为30 min左右。顶管的几何参数如表1所示。整体模型的单元、结点数量如表2所示。
2.2 材料参数的选择以及边界条件的设置
顶管实际穿越的土层主要为砂质黏土，数值模拟土体材料简化为同一种材料，模型土体重度为20 kN/m3、弹性模量为20 MN/m2、摩擦角为0.2°、流动应力比为0.935、膨胀角为10.4°。土体采用Drucker-Prager模型，模型参数参考砂质黏土室内土工试验结果。机头的材料为钢材，本模型机头部件为实心，而实际的机头内部为非实心，为了使模型与实际施工过程保持一致，令模型中机头的密度为实际钢材的1/3，即2.62×103 kg/m3。顶管的材料为钢筋混凝土，顶管及机头的材料参数如表3所示。整体模型地面为自由边界，模型中与x轴垂直边界施加x方向的约束，与y轴垂直的边界施加y方向的约束，底部边界施加x、y、z方向的约束。
表3 顶管及机头材料参数
Tab. 3 Material parameters of pipe jacking and nose
[部件 弹性模量 /
（GN/m2） 泊松比 重度 / （kN/m3） 顶管 33.4 0.2 26.0 机头 200.0 0.3 26.2 ]
2.3 顶管顶进以及接触问题
为了获得不同顶进位置竖向曲线顶管对应的顶进力，建立不同顶进位置的顶进力估算分析模型。计算步骤如下：
（1）初始地应力平衡，消去由于土体自重应力引起的位移。
（2）消除机头端面厚1.5 m、长9.1 m、高5.5 m的土体单元，同时对起始管节端面施加1.5 m的位移荷载；如果顶进位置超过顶管轴向转折位置，在对起始管节端面施加1.5 m位移荷载的同时，还需要对经过转折处的管节施加角度为顶管坡角大小的刚体转动荷载。
（3）顶进完成后，提取起始管节端面的轴向单元结点力或应力，对轴向单元结点力求和或轴向应力对顶管截面面积求积分获得对应顶管顶进位置的顶进力。
（4）重复步骤（1）~步骤（3），计算得到不同顶进位置的顶进力，得到顶进力随顶进距离的拟合曲线。
顶管顶进过程中，顶管侧壁与围岩之间形成相对滑动位移，产生滑动摩擦力，该滑动摩擦力形成顶管的侧摩阻力，为顶管顶进力的主要组成部分。影响侧摩阻力的主要因素为接触面的摩擦系数和法向注浆压力。因此，在模型中机头、管节侧壁与围岩之间的接触面设置滑动接触，机头、管节侧壁为主接触面，围岩土体表面为次接触面。通过设置不同的滑动摩擦系数和管节侧壁法向注浆压力得到不同条件下的顶管顶进力。对于大尺寸矩形截面顶管，很难形成完整的泥浆套，仅在顶面和侧面形成浆液减阻面，通常认为顶管底面与土体直接接触。为了增加模型的收敛性，采用等效滑动摩擦系数［16］将底面的滑动摩擦力等效到矩形顶管的顶面、侧面和底面，提高模型计算的收敛性。本研究计算了等效滑动摩擦系数为0.20、0.30、0.38条件下的顶进力。
同时，对于矩形顶管，注浆压力通常只存在于顶管的顶面和侧面，因此，模型中的注浆压力仅施加于顶管的顶面、侧面以及与之接触的围岩土体内表面。根据实际工况，起始区段（0~27.0 m）的注浆压力为0.15 MPa，中间区段（27.0~45.0 m）的注浆压力为0.10 MPa，接收区段（45.0~91.5 m）的注浆压力为0.05 MPa。对应区段顶进位置采用对应区段的注浆压力。
2.4 计算结果分析
每次顶管顶进计算完成时，顶管在顶进过程中会与周围的土体产生新的应力场。提取每个模型起始管节截面的轴向应力或轴向单元结点力分量，并进行相应的积分或求和，得到对应位置的顶管顶进力。若选择轴向应力分量，按式（1）对所有轴向应力分量取平均值，并对截面面积求积分计算得到顶进力。若选择轴向结点力分量，按式（2）对所有结点力分量求和得到顶进力。
[Fs=(σin1+σ′in2)×A] （1）
[Fs=Fn] （2）
式中：Fs为顶管顶进力（kN），[σi]为负结点正应力分量（kPa），A为顶管管节截面面积（m2），n1为负正应力分量的总结点数，[σ′i]为正结点正应力分量（kPa），n2为正正应力分量总结点数，Fn为对应结点编号的结点力分量（kN）。
在等效滑动摩擦系数为0.20、0.30、0.38的条件下，提取整体位移控制有限单元模型结点力，计算得到不同顶进距离的顶进力，如表4所示。
将整体位移控制有限单元法与文献［15］提出的分段位移控制有限单元法在侧摩擦系数为0.20、0.30、0.38条件下进行顶进力对比分析，结果如图5所示。
图5显示，相同侧摩擦系数条件下，采用整体位移控制有限单元法计算得到的顶进力通常大于采用分段位移控制有限单元法计算得到的顶进力。当侧摩擦系数为0.20、0.38时，整个拟合区段，前者顶进力均大于后者；当侧摩擦系数为0.30时，中前区段（顶进距离≤30 m）前者顶进力大于后者，中后区段（顶进距离>30 m）前者顶进力小于后者。分段位移控制有限单元法通过轴向坡度将竖向曲线顶管划分为若干区段，分别建立直线位移控制有限单元分析模型，区段与区段之间通过设置特定工况条件下的位移边界条件建立区段与区段间的力学联系。分段位移控制有限单元法仅考虑顶管顶进到轴向转折处引起的转折处所在截面上各单元结点的力学响应量的改变，未考虑管节经过转折处刚体转动角度对截面上各单元结点力学响应量的改变。而整体位移控制有限单元法，无须设置特有的边界条件，区段与区段之间的力学联系是自然存在的，同时，自然考虑了通过轴向转折处管节转动引起的顶进力的变化，计算精度更高。计算结果表明，在计算顶进力时不能忽略管节经过转折处的刚体转动位移，该位移对于顶进力的影响较大，一定程度上增加了大截面竖向曲线顶管的顶进力。
3 结 论
（1）与分段位移控制有限单元法相比，整体位移控制有限单元法估算大尺寸矩形截面竖向曲线顶管顶进力具有计算过程简单、力学关系清晰及计算精度高的特点。该模型将竖向曲线顶管视为一个整体，通过对经过轴向转折处施加坡角大小的转动位移和起始管节截面的一个管节宽度的平动位移实现竖向曲线顶管的整体顶进，顶进完成时，对起始管节截面单元结点力的轴向分量求和叠加，该叠加值作为顶管顶进到某个位置的顶进力。
（2）与分段位移控制有限单元法相比，相同侧摩擦系数条件下，整体位移控制有限单元法计算得到的顶管顶进力通常稍大于分段位移控制有限单元法，两者的整体变化趋势一致。这是因为竖向曲线顶管管节经过轴向转折处的刚体转动角度对顶进力的影响所致。分段位移控制有限单元法仅考虑了曲线顶管的区段坡度对顶进力的影响，未考虑经过轴向转折处的管节转动过程对顶进力的影响，这是两者差异的主要原因。当侧摩擦系数为0.38时，2种计算方法的计算结果均能反映实际顶进力的变化过程，但是整体位移控制有限单元法的计算模型简单，工程应用更有优势。