新建地铁隧道开挖打破了周围地层的应力状态，使地埋管线产生附加变形．管线变形过大将影响其正常使用，甚至引发严重的工程事故[1]．管线的纵向应变常作为管线的安全评价指标[2]．地埋管线纵向变形的理论计算方法主要基于弹性地基梁理论，将管线视为弹性地基梁，地层变形作为附加荷载．文献[3]采用修正Peck曲线作为隧道施工引起的地层沉降曲线，给出了管线变形的连续弹性地基梁法．文献[4]基于连续线弹性地基梁法，计算了管线在隧道施工影响下的弯曲应变．文献[5]给出了适用于管线变形计算的Winkler地基模量表达式．以上研究在管线荷载计算时只考虑了地层沉降，实际上，地层位移同时存在竖向与水平分量，除使管线产生挠度外，也会使管线产生轴向变形．另外，管线挠曲也会使管线伸长而产生轴力．目前，考虑管线轴力的管线变形计算方法多为小变形方法[2,6]，将挠曲和轴向变形单独计算，没有考虑挠度对轴向变形的影响．文献[7-8]考虑了管线挠曲引起的管线轴向伸长，但没有计算地层水平荷载引起的管线变形．为全面探讨地埋管线轴向变形行为，将管线视为弹性地基梁，考虑地层竖向和水平位移的共同作用，引入管线变形的几何非线性，建立管线挠度和轴向变形联立的控制微分方程组，采用优化方法进行求解．通过案例计算和离心模型试验，将有限元计算结果、试验结果同优化解进行对比，验证了优化方法的正确性．1 考虑轴向变形的控制微分方程1.1　平衡微分方程隧道垂直下穿既有地埋管线施工使管线受到竖向和水平向的地层位移荷载．将地埋管线视为弹性地基梁，计算模型如图1所示，图中：qv和qh分别为竖向和水平向分布荷载；mh为分布弯矩．10.13245/j.hust.210322.F001图1计算模型采取以下计算假设：a.管线和地层均为匀质线弹性材料，不考虑管线接头的影响；b.管线变形服从平截面假定；c.基底竖向变形服从Winkler假定，管土摩擦力与管土水平相对位移呈正比．管线微段变形后的受力状态见图2，Fn和Fs分别为轴力和剪力；M为弯矩；θ为转角；fv和fh分别为垂直和平行于管轴线的荷载．由管线微段的受力平衡可得：10.13245/j.hust.210322.F002图2管线微段受力示意图dFn+Fsdθ=-fhdx;dFs-Fndθ=fvdx;dM-Fsdx=mhdx. (1)管线变形后，轴线方向会因挠曲偏离水平方向，若管线变形后截面转角为θ，则有：fv=qvcos θ-qhsin θ；(2)fh=qhcos θ+qvsin θ．(3)1.2　管线荷载根据计算假设，管线受到竖向荷载为qv=KD(ws-wp)，(4)式中：K为地基系数；D为管线外径；ws为地层沉降；wp为管线挠度．根据计算假设，管周任意点受到水平荷载为qh,u=Kh(us-up)，(5)式中：Kh为管土摩擦系数；us为地层水平位移；up为管线任意点水平位移．根据平截面假定，管线水平位移沿截面高度呈线性变化，则管截面上任意点水平位移为up=u0-rθsin α，(6)式中：up为管线水平位移，其中管轴线水平位移记为u0；r为管截面任意点到截面中心的距离；α为管截面任意点与截面中心的连线同水平面的夹角．将式(6)代入式(5)，对管线外周各点管土摩擦力进行积分，得到管线水平荷载qh=∫02πqh,udC=KhC(us-u0)，(7)式中C为管线外圆周周长．将管线外周各点摩擦力对横截面中性轴取矩，得到管线分布弯矩mh=rsin α∫02πqh,udC=0.5KhCro2θ，(8)式中ro为管线外半径．1.3　几何方程和物理方程考虑管线挠度对管线水平位移影响时，管线任意点纵向应变为[9]εp=u0'-rsin αwp''+0.5(wp')2．(9)管线轴向应力与轴向应变的关系为σp=Eεp，其中E为管线弹性模量．则管线轴力和弯矩为Fn=∫AσpdA=EA[u0+0.5(wp')2]；(11)M=rsinα∫AσpdA=EIwp'''，(12)式中：I为管线截面惯性矩；A为管截面面积．将式(12)代入方程组(1)，得到管线剪力Fs=EIwp'''-0.5KhCro2wp'．(13)1.4　控制微分方程组转角与挠度的微分关系为θ=wp'．(14)将式(4)、(7)、(8)和(14)代入式(2)和(3)，再将式(2)、(3)、(8)、(11)、(12)和(13)代入方程组(1)，得到管线沉降与水平位移联立的控制微分方程组：       EIwp(4)=KDcoswp'(ws-wp)-KhCsinwp'(us-u0)+0.5KhCro2wp''+EA[u0'wp''+0.5(wp')2wp''];       EAu0''=-KDsinwp'(ws-wp)-KhCcoswp'(us-u0)+0.5KhCro2wp'wp''-EIwp''wp'''-wp'wp''. (15)2 控制微分方程的优化解法方程组(15)的非线性程度较高，寻找方程组的解析解是困难的，这里采用优化方法进行求解，通过构造目标函数求出未知边值的优化解[10-11]．令y1=wp，y2=wp'，y3=wp''，y4=wp'''，y5=u0，y6=u0'，将方程组(15)化为一阶常微分方程组为 y1'=y2; y2'=y3; y3'=y4;       y4'=KDcosy2EI(ws-y1)-KhCsiny2EI⋅ (us-y5)+KhCro22EIy3+AIy3y6+12y22y3; y5'=y6;       y6'=-KDsiny2EA(ws-y1)-KhCcosy2EA⋅ (us-y5)+KhCro22EAy2y3-IAy3y4-y2y3. (16)设计算区间为[-l，l]，因区间两端受隧道开挖影响较小，令两端水平位移、转角、剪力均为零，则方程组(16)的已知边值为y2-l=y2l=0，y4-l=y4l=0，y5-l=y5l=0．设x=-l端未知边值y1-l=x1，y3-l=x2，y6-l=x3，结合x=-l端已知边值采用四阶龙格-库塔法推算x=l端的边值y˜2l，y˜4l，y˜5l，将求得边值与x=l处的已知边值进行比较，得到目标函数为Z=(y˜2l-y2l)2+(y˜4l-y4l)2+(y˜5l-y5l)2．(17)目标函数Z的极小值点即为未知边值x1，x2，x3的最优解．由于龙格-库塔法有较多计算步，难以得到目标函数的显式表达，因此引入粒子群算法寻找目标函数极值[12]．粒子群算法的计算参数包括粒子个数n，学习参数c1和c2，惯性权重w，粒子空间维度d(未知边值数目)．粒子群算法将目标函数的解抽象为n个粒子，每个粒子在d维空间运动，位置矢量为xj(下标j为粒子编号)，速度矢量为vj．粒子根据自身经历的最好位置pj和粒子群总体经历的最好位置pg,j更新自身的位置和速度，具体为vj(t+1)=wvj(t)+c1r1[pj-xj(t)]+c2r2[pg,j-xj(t)]；(18)xj(t+1)=xj(t)+vj(t+1)，(19)式中：t为迭代次数；r1和r2为0~1均匀分布的随机数．求解控制微分方程组的计算流程为：a.生成n个随机粒子(也即生成n组由x=-l端未知边值组成的3阶向量)，赋予每个粒子初始位置和初始速度； b.结合x=-l端已知边值，采用龙格-库塔方法计算x=l端边值，利用式(17)计算目标函数值；c.根据步骤b所得目标函数值，记录粒子最好位置pj和群体最好位置pg,j，通过式(18)和(19)更新粒子速度vj和位置xj；d.重复步骤b和步骤c达到最大迭代次数，得到x=-l端未知边值最优解．3 案例计算及试验验证采用Matlab编写优化算法的计算程序．采用ANSYS建立有限元模型作为对比，管线由壳单元Shell63模拟，管土相互作用由弹簧单元Combin14模拟．管线与弹簧单元材料属性均为线弹性．弹簧单元一端与壳单元节点连接，另一端施加地层位移荷载，管线两端约束纵向位移．管线参数为：弹性模量70 GPa，外径0.5 m，壁厚0.018 m．土层参数为：最大沉降136 mm，沉降槽宽度系数2.6 m，地基系数2.38×107 Pa/m，管土摩擦系数8.34×106 Pa/m，管轴线与隧道中心的高差3.5 m．龙格-库塔方法计算步长取5×10-3 m，计算范围l=20 m．粒子群方法粒子数目n=40，学习系数c1=c2=2，惯性权重w=0.5，粒子空间维度d=3，最大迭代次数取1 000．隧道垂直下穿管线引起的管线周边地层沉降及水平位移为[13]ws(x)=ws,0exp[-x2/(2i2)]；(20)us(x)=(x/zr)ws(x)，(21)式中：ws，0为地层最大沉降；i为沉降槽宽度系数；zr为管线中心与隧道中心的埋深差．图3给出了优化方法和有限元法计算结果，试验装置及传感器布置如图4所示．可见两种方法得到的计算结果基本一致．经过寻优计算，目标函数值Z=2.58×10-28，说明优化方法有很高的计算精度．试验采用交通运输部天津水运工程科学研究院的TK-C500型土工离心机．试验设计加速度为80g(g为重力加速度)．根据π定理[14]，模型几何相似数为80．10.13245/j.hust.210322.F003图3管线变形和内力10.13245/j.hust.210322.F004图4试验装置及传感器布置(mm)1—LVDT(管线顶部)；2—LVDT(管轴线同水平处)；3—直线导轨．试验模拟隧道垂直下穿既有管线，采用线性可变差动变压器(LVDT)位移传感器测量管线沉降及管轴线同一水平处土层沉降；设置9个弯矩测量断面和9个轴力测量断面．试验管线为铝合金管，弹性模量69 GPa，长1 125 mm、外径17 mm、壁厚1 mm．试验管线两端用直线导轨约束，仅能够沿竖向自由滑动．采用外套钢套筒的液压油缸模拟隧道开挖，油缸外径74 mm，套筒外径76 mm、内径74 mm．套筒可沿油缸纵向滑动，推出套筒将引起地层损失从而使地层及管线产生变形．试验用土为丰浦砂，按0.8的相对密实度撒入模型箱．室内加载试验测得地基系数3.66×106 Pa/m，管土摩擦系数1.56×106 Pa/m．管线埋深62.5 mm，隧道埋深145.5 mm．如图5所示，采用沉降拟合数据计算管线沉降、弯矩、轴力，将管轴线同一水平处土体沉降数据按式(20)拟合，得到ws,0=126.8 mm，i=6.28 m(已换算为原型沉降)．10.13245/j.hust.210322.F005图5优化解与试验数据的对比由图5可见：管线沉降、弯矩、轴力计算值与试验数据趋势一致．离心模型试验测得管土最大相对位移δru*是36.7 mm，发生在隧道上方，为管线相对土体向上位移．管土相互作用力与管土相对位移的关系可简化为理想弹塑性模型[15]，当管土向上、向下相对位移分别达到δru和δrd时，管土相互作用力达到限值．根据δru的取值范围为0.005~0.015H(H为管轴线埋深)可得δru取值范围为28.4~85.2 mm．由于δru*超过δru取值下限但远未达到上限，因此将本试验中管土相互作用近似为线弹性，则本文优化解与试验结果相互能够进行较好的印证．4 影响因素分析本文方法与既有方法的不同之处在于，管线荷载计算中考虑了地层变形的水平分量，几何方程中考虑了管线挠度对轴向变形的贡献．下面结合案例计算参数分析地层水平荷载和几何非线性的影响．4.1　地层水平荷载的影响图6给出了管线顶部与底部纵向应变，与不考虑地层水平荷载得到的管线应变进行对比．由图6可见：管线在穿越中心纵向应变最大，最大压应变位于管线顶部，最大拉应变位于管线底部．与仅考虑竖向荷载的结果相比，本文方法计算得到的最大压应变明显增大，最大拉应变明显减小．说明水平荷载引起的轴向应变是管线纵向应变中不可忽略的一部分，仅考虑竖向荷载计算得到的管线纵向应变将高估管线的最大拉应变，低估管线的最大压应变，因此在管线荷载计算中考虑地层水平荷载是有必要的．10.13245/j.hust.210322.F006图6管线顶部及底部纵向应变图7给出了不同管土摩擦系数计算得到的管线轴力．由图7可见：不同管土摩擦系数下管线的受力状态不同．管土摩擦系数为0时，管线全长受拉，这是由于管线挠度使管线轴向伸长所致．随管土摩擦系数增大，管线两端拉力增大，中部拉力减小，由受拉状态转变为受压状态，且压力随管土摩擦系数的增大而增大．这是因为管土摩擦系数越大，地层对管线的水平作用力越大，使管线中部由受拉状态转变为受压状态．由此可见，管土摩擦系数较小时，通过控制管线拉应变以保证管线的安全；管土摩擦系数较大时，在控制管线拉应变的同时，也须要对管线压应变进行控制．10.13245/j.hust.210322.F007图7不同管土摩擦系数计算的管线轴力4.2　几何非线性的影响令地层最大沉降为300 mm，地层水平位移为0 mm，计算仅由地层沉降引起的管线轴力，如图8所示．图8中地层沉降引起的管线拉力达526.35 kN，可见当地层沉降较大时，管线挠度对管线轴向变形的影响是不能忽略的．10.13245/j.hust.210322.F008图8地层沉降引起的管线轴力令地层最大沉降为30和300 mm，分别计算管线轴力和弯矩，进行以下两组对比：a.与仅考虑地层水平荷载下的管线轴力进行对比，了解地层竖向荷载对管线轴力的影响；b.与仅考虑地层竖向位移荷载下的管线挠度和弯矩进行对比，了解地层水平荷载对管线挠度和弯矩的影响．两组对比结果分别列于表1和表2，其中：Ft和Fc分别为管线拉力和压力；Ms和Mh分别为管线正弯矩和负弯矩．10.13245/j.hust.210322.T001表1管线轴力ws，0/mm地层荷载条件Ft，max/kNFc，max/kN30竖向与水平向126.10-562.44仅水平向120.84-569.26300竖向与水平向1 754.83-4 990.05仅水平向1 224.73-5 769.3610.13245/j.hust.210322.T002表2管线弯矩ws，0/mm地层荷载条件Ms，max/(kN·m)Mh，max/(kN·m)30竖向与水平向150.20-69.83仅竖向149.76-70.20300竖向与水平向1 608.73-705.18仅竖向1 487.84-695.90由表1和表2可见：地层最大沉降为30 mm时，管线轴力和弯矩差别不大，这是因为地层变形较小时，管线几何非线性变形不明显，此时管线挠曲和轴向变形的相互影响较小．地层最大沉降为300 mm时，管线产生了明显的几何非线性变形．由表1可见：同时考虑竖向和水平荷载情况下，管线最大拉力比仅考虑水平荷载时增大43.3%，管线最大压力减小13.5%，这是因为管线挠度较大，使管线产生了较大的轴向伸长和轴力．由表2可见：同时考虑竖向和水平荷载的情况下，与仅考虑竖向荷载相比，管线最大正弯矩增大8.1%，这是因为水平荷载增大了管线挠曲线在穿越中心线附近的曲率，因而管线弯矩有所增大．可见管线挠度较大时，管线弯曲变形和轴向变形是相互影响的．传统管线变形计算方法不考虑管线变形的几何非线性，管线挠度和水平位移各自独立，仅与位移方向上的荷载有关[6]，不能计算管线挠度引起的管线长度改变，也无法反映弯曲变形与轴向变形的相互影响，因此，相比传统计算方法，本文方法更为合理．5 结论a.建立了管线挠度与轴向变形联立的管线变形控制微分方程组，采用优化方法求解．进行了案例计算及离心模型试验，将有限元计算结果、试验结果同优化解进行对比，验证了优化方法的正确性．b.地层水平荷载引起显著的管线拉压应变，在管线变形计算中是有必要考虑的．管土摩擦系数增大导致管线受力状态由全长受拉转变为部分受拉部分受压，在管土摩擦系数较大的地层中须对管线拉压应变同时进行控制．c.地层变形较小时，管线弯曲变形和轴向变形相对独立．地层变形较大时，管线的几何非线性变形比较明显，表现为管线挠度引起较大的轴力和轴向变形，且管线弯曲和轴向变形存在相互影响．