洪水作用下穿河悬空输气管道的动力响应研究

符炜东 田宇航 谭亮

(1.四川大学水力学及山区河流开发保护国家重点实验室，四川 成都 610000；2. 四川省广元水文水资源勘测局，四川 广元 628000）

穿河管段的安全性问题是能源领域非常典型的工程问题。针对穿河管道，人们起初关注的是管道的合理埋深，在保证管道安全的基础上力求减小工程量，从而达到节省工程投资的目的。为此，国内外学者总结了一些在一定条件下合适的计算冲刷深度的公式：如1981年，俞乐群借用桥渡建筑物附近一般冲刷深度的公式来计算管道穿越河流处的河床冲刷深度[1]；2011年Azamathulla 将人工神经网络方法带入到管道冲刷计算中[2]。有学者针对一些实际工程作出安全性分析：1998 年，黄金池提出了一些穿河管段保护措施[3]；2002年，先智伟分析了水下穿越事故的主要特征，指出了穿越事故主要原因，并提出相应防护措施[4]。随着计算机技术的发展，人们开始利用各类软件对管道的荷载进行研究：2012 年，徐涛龙等对穿越管道临界悬空高度进行动态数值模拟[5]；2013 年，姚安林等使用数值模拟的方法，对不同壁厚穿河管段的临界悬空高度进行了研究[6]；

这些学者大多关注管道的埋深、管道悬空壁厚、管道的悬空长度、管道的应力特性等单一参数，而很少关注洪水与管道之间产生流固耦合致使管道发生形变以及管道的应力特征改变。对跨河管道进行数值模拟时，通常将河道简化为矩形水槽，忽略了河床的地貌变化及河流自由水面的模拟，致使水流的流态与实际间存在差异，因而对于实际工程问题，数值模拟得到结果的准确性难以保证。为此，本文使用管道所在河道横断面数据建模，基于Workbench 软件对某山区穿河管工程进行流固耦合仿真分析，得出在一定的冲刷深度下，不同流速的洪水及管道壁厚所对应的悬空管段的应力特性及管道破坏模式。

1 计算模型

1.1 控制方程

对于一般的可压缩牛顿流体，守恒定律通过如下控制方程描述控制方程包括质量守恒方程、动量守恒方程、能量守恒方程，流体运动方程的具体表达式如下：

质量守恒方程

式中:t为时间,ρf为流体密度,υ为流体速度矢量。

动量守恒方程

式中：t为时间,ρf为流体密度,υ为流体速度矢量，ff为体积力矢量。

能量守恒方程

式中：t为时间,ρf为流体密度,υ为流体速度矢量，ff为体积力矢量，qf为单位体积热量损失，E为单位质量内能。

固体运动控制方程

式中：ρs为固体密度̈为固体域当地加速度矢量，∇σs为柯西应力张量，fs为体积力矢量。

流固耦合系统

式中：τf、τs为流体、固体应力，nf、ns为流体、固体单位方向向量，df、ds为流体、固体位移。

1.2 基本原理

流固耦合问题分为有单向流固耦合和双向流固耦合。由于悬空管道形变尺度相较于流场很小，可以忽略不计，故采用单向耦合法[8]。本文依据某实际穿河工程管道所在横断面数据建立流体三维模型与管道三维模型，将流体网格导入Fluent 进行不同流速下的仿真计算，提取管道所受阻力与升力，再将Fluent的计算结果作为管道荷载导入workbench里的单项流固耦合模块(FSI)得到管道在冲击下的形变及应力特征。

1.3 三维模型的建立

1.3.1 工程概况

某输气管道全长136.5公里，管径1000毫米，管道沿线地形地貌、地质构造复杂，全线均为中低山和丘陵地貌区，山峦起伏、沟壑纵横，管道多次穿越大小河流[7]。该管道穿越一河流时，管道顺斜坡敷设，管道直径1000mm，埋深1m，采用大开挖方式，回填土较松散，地表水土流失严重，在地表形成冲沟，少数冲沟已深达1-2米。

1.3.2 流体三维模型的建立

本文采用某实际穿河工程管道所在河道横断面数据建模，管道管径1m，壁厚14mm，距进口断面50 米，距出口断面150米，管道中心点高程177.48m，埋深1m，采用穿越的方式跨过河流，与河流成90°角。假设若干年后，河床整体由于水流冲刷下切2.5m变为悬空管道。

断面数据如表1.1，冲刷后的管道与河床示意图如图1.2。研究的洪水高程为178.98m，考虑自由液面波动，使用水平线181m与冲刷后的河床断面数据建立流体三维模型，如图1.3。

表1.1 断面冲刷前后数据表Tab.2.1 Data sheet before and after cross section scouring

图1.2相对位置示意图Fig.1.2Relative position diagram

图1.3流体三维模型图Fig.1.3 3D model of fluid

1.3.3 管道三维模型建立

管道模型根据5 种不同的工况的壁厚分别进行建模，管道两端截止于河床边壁，悬空长度为40m左右。管道三维模型见图1.4

图1.4管道三维模型图Fig.1.4 3D model of pipeline

2 网格划分与边界条件

2.1 网格划分

2.1.1 流体网格划分

将流体三维模型导入Icem划分网格，对管道中间网格部分进行加密，进口段和出口段采用结构化网格划分，管道中间段由于管道距离河床底太近，且河床底部断面形状不规则，故采用非结构化网格划分。流体网格见图2.1、图2.2。

图2.1 流体网格图Fig.2.1 fluid mesh

图2.2 流体网格局部放大图Fig.2.2 Fluid partially enlarged view of the grid

2.1.2 固体网格划分

固体部分根据5种不同工况的壁厚的管道进行建模并导入Workbench中进行网格划分，管道网格划分见图2.3、图2.4。

图2.3 管道网格图Fig.2.3Pipelinemesh

图2.4 管道网格局部放大图Fig.2.4 Pipeline grid enlarged view

2.2 边界条件

2.2.1 流体边界条件

流体求解器类型选择pressure-based，速度选择absolute，时间选择transient，考虑重力，由于建立的模型中流体包括水和空气，计算模型选择多相流（Volumeof Fluid），k-epsilon（2eqn），迭代中采用SIMPLE 算法。流体模型入口边界条件为速度入口，上边界为空气出口，下边界为wall，出口边界为压力出口，耦合面为管道外壁面（wall）。时间步长采用0.01s，分别模拟流速4m/s、5m/s、6m/s、7m/s、8m/s的流场。

2.2.2 固体边界条件

大多数情况下由于土壤的刚度等因素影响，穿河管道两端既不是两端固定也不是纯铰支的情况，一般鉴于两者之间[9]，本例采用一端固定另一端简支的约束条件使结果更符合实际情况。

穿河管道采用西气东输工程常用的L485 管线钢，密度7900kg/m3，弹性模量210Gpa，切线模量13.5Gpa，泊松比0.3，屈服应力485Mpa[10]。外径1m，分别模拟11mm、14mm、17mm、20mm、23mm壁厚下管道的应力特征及形变。

3 结果分析

3.1 流场分析

从下图3.1 组分分布图可以看出当流速为4m/s 时，由于管道阻隔，洪水自入口断面到管道附近水位一直在壅高，到管道处洪水水位高于入口断面设置水位，在管道后发生明显的跌水现象，水位迅速降低，在管道下游，水位略微抬高后又小幅下降，最终流态稳定。

图3.1 组分分布图Fig.3.1 Component distribution diagram

从下图3.2流速矢量图可以看出，在洪水行进过程中，洪水到达管道附近时，由于管道阻隔，水流流态比较紊乱。从进口到接近管道过程中，水流速度从4m/s 在逐渐减小，洪水越过管道后流速增加，管道上下两端附近流速较高，最大为7.97m/s，而后流速逐渐减小后趋于稳定，比较图3.1 与3.2 可以发现，水深变化规律与流速变化规律是协调统一的，即流速增加，水深减小，流速降低，水深增加。

图3.2 流速矢量图Fig.3.2 Velocity vector

从下图3.3压力分布图可以看出管道迎水面中点受到的压强最大为31KPa，管道上端压力最小，产生负压，管道底部压力较上端大。迎水面中点受到的水流冲击最大，故所受压强最大。顶部水深较小而流速较大，根据伯努利能量方程可知，管道顶部出现负压是合理的[11]。

图3.3 压力分布图Fig.3.3 Pressure profile

流体冲击管道时，会与管道产生相互作用，作用力在来流方向的分力为阻力（FD），与来流方向垂直方向的分力为升力（FL）。阻力可分为摩擦阻力和压差阻力，其中摩擦阻力是黏性力作用的直接结果，压差阻力是边界层分离的结果，管道属于钝体，压差阻力占阻力的绝大部分[12]。升力是与来流方向垂直方向上产生的压力差，本例中由于河床与管壁距离太近，全断面平均间隙比小于0.5（管道底部与河床表面的间隙，该间隙与管道直径的比值定义为间隙比），属于高雷诺数下的近壁圆柱绕流，根据相关文献，壁面阻碍了下部剪切层与外侧流动间的相互作用，使得下侧旋涡得不到充分发展，并且间隙比越小，与壁面接触的漩涡越靠前，壁面对漩涡脱落的抑制作用越大，致使升力的时均值为正数[13-14]。流场稳定后，升力大于阻力保持在257879N 左右，与压强分布图及文献描述相同，验证了数值模拟的可靠性。

以上为来流流速为4m/s 时的流场分析，通过改变流速，得到不同流速下的管道所受升力及阻力如下表3.4。从表中可以得知，管道所受的升力FL大于所受到的阻力FD，两者均为威胁管道安全的主要因素。

表3.4 速度、升力及阻力关系表Tab.3.4 Speed,lift and resistance relationship table

从表3.4可以得出水流流速与升力、阻力关系图（如图4.5）

图3.5 流速与升力、阻力关系图Fig.3.5 Relationship between flow velocity,lift and resistance

从图4.5可以得知，水流速度增大，管道所受的升力与阻力也增大，且随着水流速度的增加，管道所受的升力与阻力增大的越来越快，管道失稳的风险也越来越大。

3.2 管道应力分析

分别将4m/s、5m/s、6m/s、7m/s、8m/s入口速度生成的荷载结果导入Workbench 中，再分别模拟11mm、14mm、17mm、20mm、23mm 壁厚悬空管道的应力及形变。以4m/s 时，工程中所使用的14mm壁厚L485管道为例。

如图3.6 所示，管道的最大有效应力产生在管道左端底部为487.88MPa，管道使用L485管线钢，管线钢为塑性材料，适用第三强度理论，即不论材料处于何种应力状态，只要最大有效应力超过材料屈服应力材料就屈服[15]。L485 管线钢屈服应力为485Mpa，可知该实际穿河管道工程采用的14mm壁厚的L485管道在河床冲刷下切2.5m后，流速4m/s时难以抵挡洪水的冲击而失稳。比较各工况下的管道最大有效应力产生处均为管道左端底部，说明失稳破坏最先发生在管道左端底部。

图3.6 管道有效应力分布图Fig.3.6 Distribution of effective stress in pipelines

如下图3.7 所示，管道的最大形变发生在距左侧管道25m处。

图3.7 管道形变分布图Fig.3.7Pipeline deformation distribution

以上为流速4m/s、管壁14mm时的工况，通过改变流速大小及管壁壁厚可以得出，该工程悬空管道的有效最大应力总是发生在管道底部左端与土壤相交处附近，这与河床断面形状以及流体形态有关。

不同流速下悬空管道的应力和形变如下表3.8

表3.8 不同水流速度及壁厚下管道最大有效应力及最大形变表Tab.3.8 maximum effective stress and maximum displacement under different water velocity and wall thickness

从表3.8 数据可的管道最大有效应力与管道壁厚的关系（如图3.9所示），管道最大形变与管道壁厚的关系（如图3.10所示）水流速度与管道最大有效应力的关系（如图3.11），水流速度与管道最大形变的关系（如图3.12所示）

图3.9壁厚与最大有效应力关系图Fig.3.9 Wall thickness and maximum effective stress diagram

图3.10壁厚与最大形变关系图Fig.3.10 Wall thickness and maximum deformation diagram

由4.9 及4.10 可以看出，管道的最大有效应力及最大形变随着管道壁厚的减小而增大，且随着管道壁厚减小，最大有效应力增大的越来越快；

图3.11流速与最大有效应力关系图Fig.3.11 The maximum effective stress and flow diagram

图3.12流速与最大形变关系图Fig.3.12 Flow rate and the maximum deformation diagram

由3.11及3.12可以看出，管道的最大有效应力及最大形变随着水流速度的增大而增大，且随着管道流速增大，最大有效应力增大的越来越快。

4 结语

本文依托某实际穿河管道工程，基于workbench平台，采用流固耦合的数值计算方法，分析了不同流速下管道周围的流场特性及管道本身的动力响应特征，研究结果表明：管道悬空后，上游段出现壅水现象，水位升高，管道处形成跌水，水位迅速下降；管道迎水面中线所受压力最大，背水面次之；本文模拟工况下的升力FL均大于阻力FD。不同工况下的最大有效应力位置均为管道左端底部与岸边交接处，根据不同工况下，壁厚、最大有效应力、最大形变及流速的关系可以得出，流速越大，壁厚越小，最大有效应力及形变越大且增长越快，管道越容易屈服破坏。