吴 彬 刘晓宇

(1.中国科学院 力学研究所流固耦合系统力学重点实验室, 北京 100190;2.中国科学院大学 工程科学学院, 北京 100049)

长输油气管道沿线地理环境复杂,不可避免地会穿越地质滑坡、坍塌、疏松土壤沉降、湿陷性黄土、地震、水毁、煤矿采空区塌陷等地质灾害隐患区.特殊地段的埋地管道在地质灾害破坏性作用下,会形成大跨度悬空管段,导致管道发生变形,严重时会拉断管道,造成严重的事故[1-2].因此,开展悬空管道受力分析研究,准确评估悬空管道的可靠性,对于管道的安全可靠运行具有重要的意义.

目前,国内外有关埋地悬空管道受力分析的研究主要基于弹性地基梁模型[3-7]和有限元模型[8],虽然形式简单,参数少,使用方便,但因其忽略了地基土层的塑性变形,不能全面地反映地基梁的实际受力与变形性态.王同涛,闫相祯,杨秀娟,等[9]和Randolph等[10]发展了弹塑性Winker模型,但将地基土体的弹性变形段简化为有限长度或不考虑弹性变形区,在连续性上存在缺陷,仍然不符合管土相互作用的实际情况.此外,由于弹塑性地基模型引入物理非线性,致使求解过程十分复杂,需要发展高性能数值求解算法.

针对以上问题,本文基于弹塑性地基梁理论,根据埋地管道悬空段、塑性段和弹性段不同加载条件,分段建立管道挠曲微分方程,采用遗传算法和简化方法,分别给出了埋地悬空管道变形和内力的理论解及塑性区长度上限值,最后通过算例进行了验证.

1 埋地悬空管道力学模型

1.1 模型概化

如图1所示,将埋地悬空管道简化为无限长连续梁,其中未悬空部分管道置于弹塑性地基上.L为悬空长度,q为垂直作用于管道的竖向荷载(包括管道上覆土压力,管道自重以及管道内介质所产生的作用),p为未悬空管道受到的被动土体抗力.考虑到模型整体具有对称性,建立以悬空段管道中点为坐标原点的直角坐标系XOY,并假设管道在几何上及受力上都对称于Y轴.

图1 悬空管道的力学模型

如图2所示,未悬空部分管道管土作用采用非线性Winkler地基模型,即地基对管道的铅垂抗力p与管道下沉位移(挠度)ω服从理想弹塑性关系

式中：k为地基弹性抗力系数;pc为单位长度地基对管道的最大支撑反力;ωc为土体从弹性至塑性的临界位移.考虑未悬空管道可能受到轴向摩擦力f作用,本模型假定：当ω<ωc时,f=0;当ω≥ωc时,可通过公式计算获得.唐永进[11]给出k、pc和f的计算公式

式中：ρs为管道上覆土密度;H为管道上覆土厚度;D为管道外径;D0为管道平均直径;t为管道壁厚;φ为管道横向位移时的土体摩擦角;μ为管道与土体之间的摩擦系数;ρP为管道重量密度;ρf为管道内介质重量密度.由式(2)和(3),可导出ωc的表达式,即

图2 荷载位移关系

由于管道在几何上及受力上都对称于Y轴,可取管道的一半作为研究对象.如图3所示,这一半管道根据不同受力状况可被划分为3个部分：

1)悬空段(0≤x

2)塑性段(L/2≤x<Δ)：这部分管道范围为G1到G2段,Δ为待定常数,由管土作用决定.此部分土体处于塑性状态,土体抗力p=pc.考虑到管道轴向摩擦力与管道轴向拉力相比是个小量,故忽略土体轴向摩擦力.

3)弹性段(Δ≤x<∞)：这部分管道范围为G2到无穷远处.此部分土体处于弹性状态,土体抗力p=kω.考虑到在摩擦力作用下,管道受到轴向拉力逐渐减小至零,且管道处于小挠度变形状态,故忽略管道受到的轴向拉力.

图3 荷载位移关系

1.2 控制方程与基本解

1)悬空段(0≤x

由力系平衡可知,管道挠曲微分方程为：

式中：E为管道的弹性模量;I为管道的惯性矩;T为管道所受的轴向荷载,是待定常数,由管土作用决定.求解微分方程(6)导出

2)塑性段(L/2≤x<Δ)

管道的挠曲微分方程为：

式中：c5～c8为待定系数.

3)弹性段(Δ≤x<∞)

管道的挠曲微分方程为：

求解微分方程(11)导出

1.3 定解条件

由对称性可知,管道在x=0处的转角和剪力为零,即

由无穷远处边界条件可知,管道在x=∞处的转角和弯矩为零,即

由连续性条件可知,管道在x=L/2,Δ处的挠度、转角、弯矩和剪力都相等,即

由于Δ为土体从弹性段发展至塑性段的转换点,因此Δ位置处管道的挠度ω2(Δ)满足如下关系

此外,整个管道在轴力T作用下的管道伸长量应该与管道的几何变形量相等,即

式中：A为管道的横截面积.

式(14)～(18)给出的14个定解条件,可求解式(7)、(10)、(12)中的14个未知量,即c1～c12、Δ和T,进而根据管道内力与挠度之间的关系

获得截面应力σx、转角θ、弯矩M与剪力Q的解析表达式.

2 模型求解

联立 式(7)、(10)、(12)、(14)、(15)、(16a)和(16b),可以导出由Δ和T表示的c1～c12表达式：

于是,悬空管道受力分析转变为求解由Δ和T为未知变量的非线性方程组：

其中：悬空段和塑性段管道挠度ω1～ω3由式(21)中c1～c12表达式代入式(7)、(10)、(12)导出.

通常,非线性方程组采用牛顿型迭代法进行求解,但这些方法每次迭代需要给定函数导数值及迭代初始值,且初值选取会影响算法的收敛性.为此,本文采用以遗传算法为代表的智能优化算法.遗传算法研究兴起于20世纪80年代末,是模仿自然选择和遗传机制的一种智能优化算法,隐含并行性和群体全局搜索是其两个显著特征,对于复杂非线性系统求解具有独特的优越性能[12-15].

如图4所示,遗传算法基本流程如下：

步骤1：编码设计.

采用二进制编码,将待求变量{Δ,T}通过线性变换,映射到[0,1]区间中的值,然后转换为二进制表示.

步骤2：种群初始化.

生成n组随机数,均匀分布在[0,1]区间上.每组随机数设定m个个体},i=1,…,n,j=1,…,m,作为初始群体的父代个体.

图4 遗传算法流程

步骤3：适应度函数定义.

在遗传算法中,每个个体的优劣体现在适应度上.目标函数为

取其倒数作为适应度函数,个体的适应度值越高目标函数值越小.

步骤4：构造选择算子.

采用随机遍历采样方式,从初始群体的父代个体中随机选择,生成第1代子个体}.

步骤5：构造变异算子.

为防止群体进化陷入搜索空间的某个局部区域,采用离散变异方式,产生第2个子代群体}.

步骤6：构造交叉算子.

步骤7：适应度判别.

计算当前个体{Δ(i,j),T(i,j)}适应度函数值P*,若P*≤P,则结束遗传算法迭代,记录最优个体数值}=},反之,}=}*.

步骤8：演化迭代.

3 塑性区上限长度求解

当不考虑管道所受的轴向荷载T时,悬空段和塑性段管道的挠曲微分方程分别简化为

式中：c1～c8为待定系数.

联立式(25)、(12)、(14)、(15)、(16a)、(16b)和(17),略去中间推导过程,可导出不考虑轴向荷载条件下悬空管道塑性区长度Δu的解析表达式,即

鉴于管道轴向荷载能够有效延缓悬空管道下沉挠度,不考虑轴向荷载条件下悬空管道塑性区长度一定大于考虑轴向荷载条件下悬空管道塑性区长度,因此,Δu在物理上可看作悬空管道塑性区上限长度.此上限解一方面可用于遗传算法中Δ搜索域的确定,另一方面也可用于工程设计分析时塑性区的估算.

4 算例验证

采用王小龙,姚安林[5]中埋地钢管算例：管道材料为X70管材,管道弹性模量E=210GPa,外径D=1.016 m,壁厚t=0.017 5 m,截面惯性矩I=6.84×10-3m4;管道埋深H=1.5 m,自重4.282 k N/m,上覆土体对管道的垂直作用力为48 k N/m.管道悬空长度L=30m、100 m、200 m,地基弹性抗力系数k=4.064×104k N/m2,土体弹塑性临界位移ωc=0.037 7 m.

采用Matlab软件及其遗传算法工具箱(GATBX)[16-17]对式(22)中的Δ、T变量进行优化求解,计算管道的挠度和内力.选取的种群大小n=20,种群中个体的数目m=10,每个种群的长度取20,遗传代沟取0.9,最大遗传代数取100.选用随机遍历采样的选择函数,离散变异算子和单点交叉算子,交叉概率为0.7,变异概率为0.006,Δ搜索域为(L/2,Δu).

为验证理论解和近似解析解的准确性和适用性,采用通用有限元软件ABAQUS对悬空长度L分别为100 m、200 m和300 m的3种工况开展有限元数值模拟.参考文献[10,18-19]中建模方法,管道选用Timoshenko梁单元(B21),地基土层采用管-土相互作用单元(PSI24).其中,100 m算例中共划分420个PSI24单元,200 m算例中共划分800个PSI24单元,300 m算例中共划分1 200个PSI24单元.管道计算总长度为1 500 m,共划分1 500个B21单元,且左右端点固支.

图5～7分别为L=100 m、200 m和300 m 3种工况下,弹性地基梁模型(考虑轴向荷载、不考虑轴向荷载)、弹塑性地基梁模型(考虑轴向荷载、不考虑轴向荷载)、有限元模型等5类模型计算出的管道右半段挠度ω、弯矩M、最大应力分布曲线.具体计算数值见表1.

图5 L=100 m工况时各模型计算曲线图

图6 L=200 m工况时各模型计算曲线图

图7 L=300 m工况时各模型计算曲线图

表1 3种悬空长度下不同模型计算结果

可以看出,3种工况下弹塑性地基梁模型(考虑轴向荷载)始终与有限元结果吻合很好;弹性地基模型(考虑轴向荷载)在悬空长度小于200 m时,与有限元结果吻合也较好,但随着悬空长度的增加,最大应力逐渐偏离有限元结果,计算结果趋向保守;未考虑轴向荷载的力学模型计算结果与有限元结果差异巨大,不适于工程设计与分析;因此,在悬空管道受力计算时,必须考虑管道轴向荷载的影响,同时当管道出现大跨度悬空时,宜采用弹塑性地基梁模型进行计算分析,避免过高估计悬空管道的受力情况.

表2为不同悬空长度的情况下,分别采用遗传算法和解析式(26)计算出的悬空管道塑性区长度精确解和上限解.可以看出,悬空长度为100 m时,上限解约为精确解的2倍左右;随着悬空长度增加,上限解与精确解的差异性增大;悬空长度为300 m时,上限解约为精确解的4倍左右;因此,采用式(26)进行悬空管道塑性区估算时,实际值约为估算值的1/2至1/4.

表2 不同悬空长度下管道塑性区长度

5 结 论

1)基于弹塑性地基梁理论,将埋地悬空管道划分为悬空段、塑性段和弹性段,分段建立基本微分方程及定解条件.

2)提出一种半解析半数值的高效计算方法,先将悬空管道受力分析转变为求解由弹塑性临界点Δ和轴力T为优化变量的两个非线性方程组问题,然后采用遗传算法进行多目标全局优化搜索,计算结果与有限元数值模拟吻合很好.并且,在不考虑管道轴向荷载的条件下,推导出悬空管道塑性区长度的上限值,方便遗传算法计算与工程设计分析.

3)算例分析表明,当管道出现大跨度悬空时,弹性地基梁模型给出的管道内力显著高于弹塑性地基梁模型结果,分析结果偏于保守.因此,在地质灾害易发区,宜采用弹塑性模型开展管道强度设计和安全评估,以减少管道建设及运营成本.