带接头管线变形计算的传递矩阵法

程霖，杨成永，路清泉，马文辉，车敬珂

（1.北京交通大学 土木建筑工程学院，北京 100044；2.中国建筑第二工程局有限公司，北京 100160；3.北京市轨道交通建设管理有限公司，北京 100068；4.北京市市政工程设计研究总院有限公司，北京 100082）

地铁隧道施工使得邻近既有地下管线产生附加变形，造成管线破损，甚至引发工程事故[1].为保证隧道开挖不影响管线的正常使用，需要对管线进行安全评价，进而完善设计与施工方案.

地下管线的安全评价指标为地铁施工引起的附加变形[2].匀质管线变形的计算方法，一般将管线简化为刚度连续的匀质杆件[3]，承受隧道开挖引起的地层附加位移荷载.由于管线接头力学性能复杂，带接头管线的变形计算常采用经验方法，假设管节为完全刚性，接头能够自由转动，管线接头沉降与地层沉降一致，根据几何关系推算接头相对转角[2，4].经验方法忽略了管线与地层的相互作用，也没有考虑接头传递弯矩的能力.在理论计算方面，接头的存在使得管线刚度不连续，难以获得管线变形的解析解，因而多采用数值方法求解，解析计算方面的研究较少.Klar 等[5]采用边界积分法求解隧道施工引起的管线变形，将管线接头简化为具有转动刚度的弹簧铰单元.张陈蓉等[6]基于有限差分法，在管线接头两侧设置虚节点，给出了带接头管线变形微分方程的差分格式.程霖等[7]在管线变形微分方程中引入脉冲函数，得到了包含傅里叶级数的接头相对转角解析解，虽然考虑了接头处管线抗弯刚度的折减，但刚度折减系数不容易获得，因而不便于实际应用.

以上研究工作表明，目前带接头管线变形计算依赖于经验方法和数值方法，解析方法的研究尚不充分.为此，本文建立带接头管线变形的控制微分方程，从求解微分方程的角度，采用传递矩阵法[8]，推导了管节的场矩阵及管线接头的点矩阵，得到带接头管线状态变量的传递矩阵，进而求得管线变形和内力.通过与有限元解、经验方法、试验及实测数据的对比，验证了本文方法的正确性，进一步采用本文方法对管线变形的影响因素进行了参数分析.本文方法可用于计算带接头管线的变形和内力，为穿越管线施工的风险评估提供了理论参考.

1 管节和场矩阵

图1 为管线变形示意图，隧道开挖引起地层及管线产生沉降，接头处产生相对转角.

将管节视为弹性地基上的Euler-Bernoulli 梁，基底变形服从Winkler 假定，则管节变形的控制微分方程可写为[7]：

式中：Ep为管线弹性模量，Pa；Ip为管线截面惯性矩，m4；k 为地基系数，Pa/m；D 为管线外径，m；w 和S 分别为管线挠度和地层位移，向上为正，m.

地基系数k 通过式（2）进行计算[4]，即：

式中：Es为管线所在地层土体的弹性模量，Pa；νs为管线所在地层土体的泊松比.

地层位移用Peck 公式表示[9]，即：

式中：Smax为地层最大沉降，m；μ 为隧道中线的坐标，m；i 为沉降槽宽度系数，m.

如图2 所示，将一段管节划分为n 小段梁，每小段梁长度为ξ，当ξ 比较小时，每小段梁受到的荷载近似为均布荷载.

图2 管节计算简图Fig.2 Calculation diagram of pipelinesection

对图2 中任意一小段梁单独分析，将其受到的均布荷载简记为q，则该小段梁变形控制微分方程为：

将式（4）写成一阶微分方程组形式

式中：θ 为转角，以逆时针为正，rad；M 为弯矩，以使梁下侧受拉为正，N·m；Q 为剪力，以使截面右侧梁段顺时针旋转为正，N.

将式（5）改写为矩阵形式

其中：

对式（6）进行拉普拉斯变换，得到：

式中：s 为复变量；I 为4 阶单位矩阵.

对式（7）进行拉普拉斯逆变换，得到：

式中：L-1[·]为拉普拉斯逆变换运算符.

将坐标x=ξ 代入式（8），并将式（8）改写为增广矩阵形式，即：

其中：

式中：U 为场矩阵；V0、Vξ为状态向量，其元素为状态变量，下标表示状态向量和状态变量的计算位置；λ、β1、β2、β3、β4为系数，分别为：

由式（9）可见，状态向量V0通过矩阵U 传递至x=ξ 处，得到状态向量Vξ.通过上述方法获得各小段梁场矩阵，分别记为U1、U2、…、Un，则x=0 处状态变量可依次通过各小段梁场矩阵传递至x=xn处，即

2 接头和点矩阵

将接头简化为铰点，如图3 所示，图中，ks为接头转动刚度，N·m/rad；Xj为接头坐标，j 为接头编号.将管线接头分为两类，一类不能传递弯矩，两侧管节能够自由转动，将这种接头简化为自由铰，采用此类接头的管线简称“自由铰管线”；另一类能够传递部分弯矩，将这种接头简化为具有一定转动刚度的弹簧铰，采用此类接头的管线简称“弹簧铰管线”.

图3 管线接头的简化Fig.3 Simplification of pipeline joints

铰点左侧状态变量通过铰点传递至右侧，即：

铰点两侧管线的挠度、弯矩和剪力均相同，管线转角将产生突变.

对于弹簧铰，转角突变Δθ 等于弯矩与接头转动刚度的比值，即弹簧铰两侧转角的关系为：

因此，任意一个弹簧铰的点矩阵为：

对于自由铰，设第j 个接头右侧转角为：

式中：ηθ为各状态变量以及外力对管线转角的影响系数，则第j 个接头的点矩阵可表示为：

根据式（11），第j 个接头两侧状态向量关系为：

将式（14）代入式（10），得到第j+1 个接头左侧的状态向量

解方程组（16）可得影响系数ηθ的值，进而得到第j 个接头的点矩阵，即：

3 边界条件和状态变量求解

对于计算长度为L、管节数为N+1、接头数目为N 的管线，如图4 所示，管线左端（x=0）的状态变量可依次通过管节的场矩阵和接头的点矩阵传递至管线右端（x=L），即：

图4 管线计算简图Fig.4 Calculation diagram of pipeline

式中：Ue为总体传递矩阵.

计算范围内管线两端距离隧道穿越中线较远，受隧道开挖影响可以忽略不计.

对于弹簧铰管线，根据其接头特点，令管线两端转角为0，剪力为0，则左端的状态向量为：

其中，w0、M0为管线左端未知状态变量.通过式（17）将V0传递至管线右端，利用管线右端转角为0，剪力为0 的边界条件，得到关于未知状态变量的方程组

对于自由铰管线，根据其接头特点，令管线两端沉降为0，弯矩为0.则左端的状态向量为：

其中，θ0、Q0为管线左端未知的状态变量.

求解方程组（18）或方程组（19）得到x=0 处未知状态变量，再通过场矩阵和点矩阵计算管线各点的变形和内力.隧道开挖引起的管线相对转角等于接头两侧管节的转角差值，即

从求解过程看，与有限元法相比，传递矩阵法的未知量只有管线两端的未知边值，无需求解大型方程组，计算量小，且可通过增减分段数达到任意计算精度.从所需计算参数看，与文献[7]的傅里叶级数法相比，传递矩阵法采用的接头转动刚度可通过现有资料获得[10]，便于实际应用.

4 案例计算与试验验证

4.1 案例计算

采用MATLAB 编写本文方法的计算程序.为验证本文方法的正确性，采用ANSYS 建立有限元模型，其中管节采用梁单元（BEAM4）模拟，地基弹簧和管线接头采用弹簧单元（COMBIN14）模拟.地基弹簧单元的一端与梁单元节点相连，另一端施加由隧道开挖引起的地层位移荷载.

案例1 为文献[4]算例，隧道垂直下穿管线，计算参数为：地层最大沉降13.6 mm，沉降槽宽度系数2.6 m，地基系数2.38 × 107Pa/m；管线弹性模量70 GPa，外径0.5 m，壁厚0.018 m；管节长度5.49 m，管线接头简化为自由铰.

案例2 为北京地铁14 号线某区间工程垂直下穿铸铁污水管线，测得地层最大沉降12.4 mm，沉降槽宽度系数2.5 m，地基系数2.66×107Pa/m；管线弹性模量100 GPa，外径1.462 m，壁厚17.1 mm；管节长度6 m，接头转动刚度1.79×107N·m/rad.

图5 为管线沉降和弯矩的计算结果，同时也给出了与有限元结果及实测数据的对比.可见，本文方法计算结果与有限元方法结果一致，与实测数据趋势基本符合，证明了本文方法的正确性.

图5 计算结果Fig.5 Calculation results

文献[4]给出的经验方法未考虑管土相互作用，简单假定管线接头处的沉降与同水平处地层自由沉降一致，得到最大相对转角计算结果为4.42×10-3rad（0.25°）.由于管节本身具有一定的抗弯刚度，因而管线接头沉降与地层沉降并不相同，本文方法考虑了管土相互作用，所得隧道正上方接头沉降大于同水平处地层自由沉降，最大相对转角的计算结果为4.96×10-3rad（0.28°），由此可见，本文方法与经验方法相比更为安全、合理.

4.2 离心模型试验

试验采用交通运输部天津水运工程科学研究院的TK-C500 型土工离心机，设计加速度为80g（由π定理可知模型的几何相似数为80[11]），模拟盾构隧道垂直下穿既有管线，如图6 所示.

如图6 所示，采用LVDT 位移传感器测量管线沉降及管轴线同一水平处土层沉降.设置11 个弯矩测量断面，每个测量断面用4 个应变片连接成全桥电路.

图6 试验装置及传感器布置（单位：mm）Fig.6 Arrangement of test device and sensors（unit：mm）

试验管线为承插式铝合金管.试验隧道开挖采用外套钢套筒的液压油缸模拟，套筒可沿油缸纵向滑动，推出套筒将引起地层损失从而使地层及管线产生变形.将管线与隧道的模型及原型试验参数列于表1.

表1 管线与隧道试验参数Tab.1 Test parameters of pipeline and tunnel

管线模型包含8 个管节，总长922.5 mm，管线模型中心接头与隧道开挖中线重合.在管节承口内部粘贴橡胶圈，将管节插口插入橡胶圈粘牢，如图7 所示，室内试验测得接头旋转刚度为4.76×108N·m/rad.试验用土为丰浦砂，室内试验测得地基系数为7.37×106Pa/m.

图7 试验管线的接头构造（单位：mm）Fig.7 Detail of joint of the test pipelines（unit：mm）

试验时，逐级增大离心机加速度至80g，待管线与土层沉降稳定，将钢套筒按40 mm/min 的速度匀速推出，模拟隧道开挖，管土沉降稳定后逐渐降低离心机加速度至停止.

将隧道开挖引起的管轴线同一水平处土体自由沉降监测数据按式（3）拟合，得到Smax=73.31 mm，沉降槽宽度i=6.60 m.采用地层沉降拟合数据（已换算为原型沉降数据）计算管线沉降和弯矩，与试验结果进行对比，如图8 所示，可见，管线沉降和弯矩的理论计算结果与试验数据趋势一致，相互能够进行较好的印证.

图8 试验结果与理论计算结果的对比Fig.8 Comparison of test results with calculation results

5 影响因素分析

5.1 接头与隧道中线相对位置的影响

设接头与隧道中线距离为d，管节长度为Lp，则案例计算相当于考虑了d=0 和d=0.5Lp两种位置关系.为更全面了解接头与隧道中线相对位置对管线变形和内力的影响，令d=0、0.75 m、1.5 m、2.25 m、3.0 m，结合案例2 计算参数，计算管线变形和弯矩.

将位置坐标x、管线沉降w 和弯矩M 归一化，得到无量纲量（x-μ）/i、w/Smax和Mi2/（EpIpSmax），计算结果如图9 所示，可见，d=0 时，隧道中线处管线沉降大于地层沉降，两侧相邻管节变形呈现刚体转动，管线弯矩较小.随隧道中线向管节中部移动，管线沉降减小，隧道中线附近管节弯曲变形增大，管线弯矩增大.图9（a）表明，距离隧道中线最近的管线接头产生的相对转角最大，该处接头相对转角是由管线沉降造成的，本文称为“沉降角”；距离隧道中线2.5i～3.5i范围内，管线产生一定的隆起变形，该处管线接头也产生明显的相对转角，本文称为“隆起角”.

图9 d 不同取值下的计算结果Fig.9 Calculation results under different values of d

图10 为管线“沉降角”和“隆起角”的最大值随d/Lp的变化规律.由图10 可见，“沉降角”和“隆起角”随d/Lp的增大都呈减小的趋势，因而，对接头变形来说，隧道中线与某一接头位置重合为最不利工况.d/Lp=0.5 时，“隆起角”大于“沉降角”，说明新建隧道中线位于管节中心时，应关注沉降槽边缘（沉降槽半宽约为3i[4]）管线接头产生的“隆起角”.

图10 接头相对转角随d/Lp 的变化Fig.10 Variation of relative angle of joints with d/Lp

5.2 地层变形及管线参数的影响

管线的力学响应受地层变形及管线自身参数的影响，包括接头转动刚度ks、管线抗弯刚度EpIp、地基系数k、管节长度Lp、沉降槽宽度系数i 和地层最大沉降Smax.

令管线接头处在隧道中线正上方，即最不利位置.将接头最大相对转角归一化，即Δθmaxi/Smax，结合案例2 计算参数，讨论以上参数对其影响.

图11 给出了归一化最大相对转角随接头转动刚度ks的变化规律.由图11 可见，接头转动刚度为0 时，接头不承受弯矩，此时接头相对转角最大；随接头转动刚度的增大，接头承受弯矩的能力增强，接头相对转角减小.

图11 接头转动刚度的影响Fig.11 Influence of the joint rotation stiffness

图12～图14 分别给出了两类管线归一化最大相对转角随地基系数k、管节抗弯刚度EpIp和管节长度Lp的变化规律.由图12 可见，接头相对转角随地基系数的增大呈先增大后减小的趋势.这是因为地基系数较小时，管节弯曲变形由接头转动释放，管节变形呈刚性，在接头两侧产生转角差异，此时地基系数增大使管线沉降增大，进而接头两侧管节的转角差异增大；地基系数较大时，随地基系数增大，管节的弯曲变形增大，接头两侧管节的转角差异减小.对于自由铰管线，相对转角在k=1×106～2×108Pa/m 范围内（大部分土质地基的地基系数位于此范围[12]）变化不大，可以推断，在土质地层中，不论土体压缩性如何，自由铰管线均易于产生较大的接头相对转角.

图12 地基系数的影响Fig.12 Influence of the foundation coefficient

由图13 可见，两类管线受管节抗弯刚度的影响相同.管节刚度较小时，管线整体变形与地层变形接近，变形曲线平缓，接头相对转角较小；管线刚度较大时，管节呈刚体转动，管线变形与地层变形差异较大，接头产生较大的相对转角.

图13 管节抗弯刚度的影响Fig.13 Influence of the bending stiffness of pipeline section

由图14 可见，两类管线受管节长度的影响相同.管线被接头划分为若干小段，管节长度较小时（Lp＜1.6i），整体变形呈现柔性，变形与地层变形接近，在这个阶段，随管节长度增长，管线整体性增强，与地层变形差异增大，接头相对转角也增大；由于管节具有一定的抵抗弯曲变形的能力，管节长度较大时（Lp＞1.6i），随管节长度增长，管线沉降减小，接头相对转角也减小；管节长度为1.6i 时，接头相对转角达到峰值.在实际工程中，应对管节长度满足这一特征的管线予以重视，加强防护和监测.

图14 管节长度的影响Fig.14 Influence of the length of pipeline section

由图12～图14 可以看出，归一化相对转角存在极限值，不会随某一计算参数的增长而无限增长.为得到归一化相对转角的极限值，计算自由铰管线归一化最大相对转角随地层变形参数Smax和i 的变化规律，如图15、图16 所示.

由图15 可见，归一化相对转角不随地层最大沉降的变化而变化，说明接头相对转角与地层最大沉降呈正比.由图16 可见，随沉降槽宽度系数的增大，归一化最大相对转角呈先增大后减小的趋势，这是因为沉降槽宽度较小时，管线由于自身刚度的原因，变形小于地层沉降；随着沉降槽宽度增大，管线变形超过地层沉降，接头相对转角增大，但沉降槽宽度的增大使得地层沉降逐渐平缓，管线变形逐渐与地层沉降趋于一致，接头相对转角减小.通过图16 可以看出，管线归一化相对转角的极限为1.1，当实际工程中缺乏设计资料时，可取Δθmax=1.1 Smax/i 作为保守的接头相对转角估计值.

图15 地层最大沉降的影响Fig.15 Influence of the maximum soil settlement

图16 沉降槽宽度系数的影响Fig.16 Influence of the width of settlement trough

6 结论

将带接头管线按接头传递弯矩的特性分为“自由铰管线”和“弹簧铰管线”两类.采用传递矩阵法求解了管线变形的控制微分方程，推导了管节的场矩阵和两类管线接头的点矩阵，给出了管线未知边值的求解过程.通过与有限元计算结果、离心模型试验结果和实测数据对比，本文方法是可靠的.

采用本文方法针对管线接头与隧道开挖中线距离、管线参数和地层变形参数进行了影响因素分析.结果表明，对接头变形来说，接头与隧道中线位置重合为最不利工况.若隧道中线位于管节中心，则距开挖中线2.5i～3.5i 的管线接头将因此处的管线隆起而产生较大的相对转角.相同条件下，自由铰管线的接头相对转角大于弹簧铰管线.在土质地层中，自由铰管线可产生比较大的接头相对转角，且相对转角基本不受地基系数影响.通过分析管节长度的影响规律，可知管节长度为1.6i 时，接头相对转角出现峰值，实际工程中应对符合这一长度特征的管线采取相应的保护措施，并加强监测.参数分析表明，管线接头相对转角的极限值为1.1Smax/i，当缺乏设计资料时，可用该值保守估计管线相对转角.