地铁隧道开挖使周围土层产生位移，引起土层中的管线产生附加变形，管线纵向变形是地下管线安全性的重要评价指标[1]．地下管线纵向变形的解析方法主要基于弹性地基梁理论，将地下管线视为土层中的梁，承受隧道开挖引起的土层附加荷载．文献[2]对比了Winkler地基模型和连续弹性地基模型求解得到的管线变形计算结果，给出了适用于分析隧道穿越施工引起管线变形的地基模量．文献[3]基于连续弹性地基模型给出了隧道施工引起管线变形的计算方法，提出了管土相对刚度，分析了管线弯矩与管土相对刚度的关系．文献[4]认识到地基模量与管土相对位移方向有关，根据Winkler地基模型推导了管线变形的控制微分方程，求得微分方程的解析解．文献[5-6]假设管线为细长杆件，建立了管线大变形的控制微分方程，求解了地层沉降引起的管线轴力．文献[7]采用基于傅里叶变换的传递矩阵法求解了管线变形的控制微分方程，分析了管线柔度系数、刚度、埋深等因素对管线变形的影响．上述研究均假定管线与土体保持接触，未考虑管土脱空．实际上，当土体位移较大时，管土脱空是可能发生的，脱空区内管线周围土体不稳定，加之地下水渗流及管线渗漏的影响，可能引起地表沉陷等工程事故，严重影响施工区周边环境安全[8]，故在管线变形计算中应对管土脱空予以考虑．目前，考虑管土脱空的管线变形计算依赖数值模拟方法[9]，其计算结果受材料本构模型、单元参数、网格质量及尺寸等因素影响较大，涉及计算参数相对较少且简便可靠的理论方法方面的研究相对较为缺乏．本研究建立了考虑管土脱空的双层弹簧弹性地基梁理论模型，推导了管线变形的控制微分方程，采用传递矩阵法进行求解[10]．将理论解与有限元解及离心模型试验结果进行对比，验证了本文方法的正确性．1 计算模型及控制微分方程传统管线变形计算方法基于单层弹簧弹性地基梁模型，不能考虑管土脱空．建立图1所示的双层弹簧弹性地基梁模型，图中：ku和kd分别为管线上部和下部土层的地基系数；L为管线计算长度；x轴和y轴分别为水平及竖向坐标轴．10.13245/j.hust.220405.F001图1双层弹簧弹性地基梁计算模型采用如下计算假设：a. 管线为匀质连续长梁，变形服从平截面假定，基底变形服从Winkler假定；b. 管土之间不存在拉力；c. 管线只承受土层荷载的作用，忽略隧道施工期间其他荷载的影响；d. 管线周围土体质量较好，管土共同变形过程中管线基底未达地基极限承载力．管线变形前受到的初始土压力通过令弹簧产生预压缩量来体现，设土层竖向压力为p，则上层和下层弹簧的预压缩量分别为δu=p/ku和δd=p/kd．在管线变形后，若管土沉降差异超过弹簧的预压缩量，则管线与土层产生脱空．根据文献[11]的试验结果，总结出管线与地基的三种接触状态．S和w分别为土层和管线沉降(以向上为正)，当w-Sδd时，管线与下部土层脱空；当S-wδu时，管线与上部土层脱空；其余情况下管线与上下部土层均保持接触．管线变形的控制微分方程[5]可写为EpIpd4wdx4=-keDw+qe，(1)式中：Ep为管线弹性模量；Ip为管线截面惯性矩；D为管线外径；ke为等效地基系数，与管土接触状态有关；qe为土层荷载．ke和qe在管土不同接触状态下的取值分为以下几种情况．a. 管土接触区段(S-δuwS+δd)为：ke=ku+kd;qe=keDS．(2)b. 管线与上部土层脱空区段(w≤S-δu)为：ke=kd;qe=keD(S+δd)．(3)c. 管线与下部土层脱空区段(w≥S+δd)为：ke=ku;qe=keD(S-δu)．(4)2 传递矩阵法将计算长度为L的管线划分为n段，每段梁长度为ξ，如图2所示．由图2可见：当ξ足够小时，每段梁所受荷载可视为均布荷载．根据弹性地基梁的初参数解，可以得到梁段的传递矩阵，即10.13245/j.hust.220405.F002图2传递矩阵法计算简图Vj=UjVj-1，(5)式中：Vj和Vj-1分别为第j段梁右端和左端的状态向量；Uj为第j段梁的传递矩阵．有Uj=β112λjβ22λj2ηe,jβ3λjηe,jβ4A-λjβ4β12λj3ηe,jβ22λj2ηe,jβ3qe,jλjηe,jβ4-ηe,j2λj2β3-ηe,j4λj3β4β112λjβ2qe,j2λj2β3-ηe,j2λjβ2-ηe,j2λj2β3-λjβ4β1qe,j2λjβ200001,式中：A=-qe,j(β1-1)/ηe,j；qe，j为第j段梁所受荷载；ηe，j为第j段梁地基系数ke,j与D的乘积，即ηe，j=ke，jD；λj和β1～β4为系数，λj=[ηe,j/(4EpIp)]1/4,β1=cos λjξcosh λjξ,β2=sin λjξcosh λjξ+cos λjξsinh λjξ,β3=sin λjξsinh λjξ,β4=sin λjξcosh λjξ-cos λjξsinh λjξ.又有Vj=(wj,θj,Mj,Qj,1)T，式中：θj为转角，以逆时针为正；Mj为弯矩，以使梁下侧受拉为正；Qj为剪力，以使截面右侧梁段顺时针旋转为正．根据图2，管线左端状态向量VL可通过各段梁传递矩阵U1，U2，…，Uj，…，Un传递至管线右端，得到右端状态向量VR，即VR=UnUn-1Un-2⋯U2U1VL=UTVL，(6)式中UT为管线的总体传递矩阵．管线两端受隧道开挖影响较小，管线两端的边界条件为：转角为0°，剪力为0 N．管土脱空须要迭代求解，迭代步骤如下．a. 当第一次迭代时，假设管土不脱空，通过式(6)将V0传递至管线右端，利用管线右端转角为0°、剪力为0 N的边界条件，得到关于管线左端未知状态变量(wL，ML)的方程组：f1(wL,ML)=0;f2(wL,ML)=0. (7)求解式(7)得到管线左端的未知状态变量，再通过各分段梁传递矩阵求得各点的变形和内力．b. 判断管土是否脱空，确定脱空范围，按式(2)~(4)对传递矩阵中ηe，j和qe，j进行相应替换，按步骤a进行第二次迭代．c. 重复以上步骤，直至前后两次迭代管土脱空范围的变化小于既定精度为止．由以上计算过程可见：与有限元法和有限差分法等数值方法相比，传递矩阵法涉及到的未知量只须控制微分方程的未知边值，不必求解大规模的线性方程组，也无须存储大量的刚度矩阵，同时具有较高的求解精度，提高了计算效率．3 算例及参数分析3.1　算例和有限元验证利用Matlab编写本文方法的计算程序，进行案例计算．利用ANSYS建立有限元模型，对本文方法进行验证，其中，管线采用梁单元BEAM4划分，单元长度与本文方法计算步长一致，地基弹簧采用不抗拉杆单元LINK180模拟，不允许管线与地基之间存在拉力，以此模拟管土脱空．计算隧道垂直下穿管线引起的管线沉降和弯矩．地铁隧道开挖引起的管轴线处地层沉降由Peck公式表示，即S(x)=-Smaxexp[-x2/(2i2)]，(8)式中：Smax为土体最大沉降；i为沉降槽宽度系数，即沉降曲线反弯点至隧道中线的距离．Smax=-136 mm，沉降槽宽度系数为i=2.6m．管线埋深为5 m，上覆土重度为18 kN/m3，则上覆土压力p=90 kPa．D=0.5 m，厚度e=18 mm．管线计算长度L=40 m，计算步长ξ=0.02 m．改变管线弹性模量Ep，ku和kd，计算5个工况．工况1取Ep=30 GPa，ku=1.19×107 Pa/m，kd=2.38×107 Pa/m；工况2取Ep=210 GPa，ku=1.19×106 Pa/m，kd=2.38×106 Pa/m；工况3取Ep=210 GPa，ku=5.95×108 Pa/m，kd=2.38×107 Pa/m；工况4取Ep=70 GPa，ku=1.19×107 Pa/m，kd=2.38×107 Pa/m；工况5取Ep=210 GPa，ku=5.95×107 Pa/m，kd=2.38×108 Pa/m．表1列出了本文方法及有限元法计算结果,表中：wmax为管线最大沉降；Mmax为管线最大弯矩；vu和vd分别为管线与上部土层和下部土层的脱空范围长度．两种方法所得结果基本一致，证明了本文方法的正确性．10.13245/j.hust.220405.T001表1案例计算结果工况wmax/mmMmax/(kN•m)vu/mvd/m本文方法有限元方法本文方法有限元方法本文方法有限元方法本文方法有限元方法1126.115126.114353.280353.281——3.363.36279.12579.112601.357601.344——3.203.203132.194132.1972 141.9902 141.67717.0016.96——4116.368116.301650.871652.7011.721.564.524.525118.378118.3632 088.8782 089.4067.487.4811.3211.28由表1的计算结果可以看出：当土体刚度较小或者管线刚度较小时，管线只与下部土体产生脱空；当土体刚度较大且管线刚度较大时，管线与下部土体产生脱空的同时，也可与上部土体产生脱空；当上部土体刚度明显大于下部土体时，将发生管线只与上部土体脱空的情况．3.2　管土相对刚度系数管线沉降和弯矩受多个参数的影响，难以得到单一参数与管线沉降和弯矩的对应关系[9]．在既有研究基础上[2，12]，提出一个无量纲的管土相对刚度系数R，表达式为R=EpIp(kukdkdku)1/(kd+ku)i4DSmaxi．(9)采用本文方法，计算不同参数组合下管线沉降和弯矩．ku的取值范围为1~50 MPa/m，kd的取值范围为1.38~138 MPa/m，i的取值范围为2~9 m，Smax的取值范围为13.6~272 mm，Ep的取值范围为4~210 GPa．按以上参数进行3 000次计算，对所得结果进行拟合，得到R与管线归一化弯矩(Mmaxi2/(EpIpSmax))和归一化沉降(wmax/Smax)的关系，如图3所示，同时也给出了不考虑管土脱空情况下计算结果的拟合曲线．10.13245/j.hust.220405.F003图3系数R对管线弯矩和沉降的影响图3中考虑在管土脱空情况下，管线归一化沉降与归一化弯矩与系数R存在以下关系：wmaxSmax=11+1.07R0.52；(10)Mmaxi2EpIpSmax=11+6.03R0.65．(11)由图3可见管线归一化沉降和弯矩均随R的增大而减小．当R1×10-4时，管线变形呈柔性，不论是否考虑管土脱空，管线归一化弯矩和沉降均趋近于1，说明此时管土沉降变形接近相等，不产生脱空，管线最大弯矩可通过土层沉降推算，即Mmax=-EpIpSmax″=-EpIpSmax/i2．(12)当R1×10-4时，两种算法管线归一化弯矩和沉降存在差异，说明此时产生了管土脱空，故R=1×10-4可作为管土产生脱空的分界线．考虑管土脱空情况下管线弯矩和沉降比未考虑管土脱空时小，这是因为管土脱空后管线承受的土体荷载与不考虑管土脱空相比较小，因而管线的变形和内力较小，由此也可以看出，不考虑管土脱空情况下，管线变形和内力的计算结果是偏保守的．当管土相对刚度系数大于1 000时，管线弯矩与沉降趋近于0，管线变形呈刚性，与下部土层脱空量达到最大，并且几乎不产生弯曲变形．此时管线弯矩可通过静力平衡计算[3]，即Mmax=kuD∫-∞∞S(x)dx4kuD/(4EpIp)4．根据式(11)和式(12)，当管土相对R位于1×10-4和1 000之间时，可采用下式估算管线沉降和弯矩，即：wmax=Smax1+1.07R0.52；Mmax=EpIpSmaxi2(1+6.03R0.65)．4 离心模型试验试验采用交通运输部天津水运工程科学研究院的TK-C500型土工离心机．模型箱净尺寸为1.15 m×0.40 m×0.70 m，试验设计加速度为80g．试验模拟6 m直径盾构隧道垂直下穿管线，试验装置及传感器布置如图4所示．采用线性可变差动变压器(LVDT)位移传感器测量管线沉降及管轴线同一水平处土层沉降；设置9个弯矩测量断面，每个测量断面用4个应变片连接成全桥电路；于管线底部设置5个土压力传感器，用于测量管底与下部土层的压力，以判断管土脱空情况．10.13245/j.hust.220405.F004图4试验装置及传感器安装(mm)试验土层原型为北京某地铁区间隧道穿越管线区域，由于开挖影响范围内土层多为砂层，故试验采用丰浦砂模拟，试验用土参数按原型土层参数的加权平均值确定，内摩擦角为31°，压缩模量为38.41 MPa，干密度为1.59 g/cm3．试验土层制备过程如图5所示．10.13245/j.hust.220405.F005图5试验土层制备室内加载试验测得地基系数ku和kd分别为2.67 GPa/m和3.56 GPa/m．试验管线为铝合金管，弹性模量为69 GPa，长为1 140 mm，外径为17 mm，以抗弯刚度相似的原则模拟隧道穿越DN1200铸铁管．管线两端用直线导轨约束，使之仅能够沿竖向自由滑动．试验隧道为外套钢套筒的液压油缸，推出套筒将引起地层损失从而使地层及管线产生变形．进行三组试验，试验条件分别为如下．a. 试验1：管线壁厚为1 mm，套筒厚度为1 mm(对应地层损失率为5%)；b. 试验2：管线壁厚为1 mm，套筒厚度为2 mm(对应地层损失率为10%)；c. 试验3：管线壁厚为3 mm，套筒厚度为2 mm(对应地层损失率为10%)．当试验时，逐级增大离心机加速度至80g，待管线与土层沉降稳定，将钢套筒按20 mm/min的速度匀速推出，模拟隧道开挖过程，开挖总耗时20 min，管土沉降稳定后逐渐降低离心机加速度至停止．将管轴线同水平处土体自由沉降数据按式(8)拟合，得到土体沉降参数(已换算为原型)，见表2．10.13245/j.hust.220405.T002表2土体自由沉降参数试验Smax/mmi/m1127.806.242262.144.023268.683.95按表2中地层变形数据计算管线沉降和弯矩，与试验结果(已换算为原型)对比，如图6所示．由图6可见：试验1与试验2相比，由于试验2土层沉降更大，因此管线产生了更大的沉降；试验2与试验3相比，试验3管线的抗弯刚度更大，故试验3管线沉降较小．管线沉降和弯矩计算值与实测值基本一致，相互能够较好地印证．10.13245/j.hust.220405.F006图6试验结果与理论计算结果对比图7给出了管线与下部土体相互作用力pr随时间t变化的时程曲线，图中：dc为测点与隧道开挖中线的距离(已换算为原型)；虚线所示处为隧道开挖完成．由图7(试验1中dc=12 m位置处土压力传感器损坏)可见：各组试验中管线与下部土体相互作用力的最小值均位于dc=0 m处(隧道开挖中线上方)，其次为dc=4 m处，说明隧道开挖中线附近管线与下部土体存在较大的沉降差异，下部土体对管线的支承作用减弱．在dc=20 m处，管线与下部土体相互作用力随隧道开挖过程较为稳定，压力值与管轴线同水平处土体自重应力相差不大，说明距离隧道开挖中线较远处，隧道开挖对管线和土体变形的影响较小．10.13245/j.hust.220405.F007图7管线与下部土体相互作用力时程曲线在dc=0 m和dc=4 m处，管线受到的下部土体作用力随隧道开挖过程呈下降趋势，表明隧道开挖中线附近管线与下部土体存在脱空趋势．在试验1中，由dc=0 m处管土相互作用力稳定数值来看，管线与下部土体未发生明显脱空．在试验2和试验3中，在dc=0 m处和dc=4 m处管土相互作用力稳定数值接近0，说明管线与下部土体在隧道开挖中线两侧产生了脱空，脱空范围长度约为8 m．由此可见：当计算管线变形和内力时，有必要对管土可能发生的脱空情况予以考虑．在dc=8 m和dc=12 m处，管线受到的下部土体作用力随隧道开挖过程呈增大趋势，压力值较隧道开挖前明显增大，这是因为脱空区内管线失去了下部土体的支承力，此时脱空区附近区域土体对管线的支承力增大，以平衡脱空区上部的土体荷载．5 结论a. 将管线视为弹性地基梁，管土相互作用简化为不抗拉的地基弹簧，建立了考虑地下管线与土层脱空的双层弹簧弹性地基梁模型．给出了管线变形的控制微分方程，采用传递矩阵法对控制微分方程进行求解，给出了管土脱空的迭代步骤．本文方法计算过程涉及参数较少，具有较高的求解精度，与有限元法相比，无须求解大规模的线性方程组，提高了计算效率．从计算假设可知：本文方法只适用于土质较好的情况，对于软弱土体，或者管线渗漏导致周围土体承载力下降的情况，是不适用的．b. 基于本文方法研究了管土相对刚度系数对管线沉降和弯矩的影响，与不考虑管土脱空相比，考虑管土脱空时管线沉降变形和弯矩均有所减小．R=1×10-4可作为判断管土能否脱空的界限．c. 管线与下部土体相互作用力的试验数据表明：在管土沉降差异较大的区域，管土之间是能够产生脱空的．不考虑管土脱空的传统计算方法偏于保守，不能正确反映管线与土体的沉降差异，在管线变形和内力计算中考虑管土脱空是有必要的．