基于有限元模拟的跨断层输水管道失效风险评价——以阿红干渠管道为例

顾世祥，梅 伟，唐 畅，杨 帆，霍玉国，曹子君

（1.云南省水利水电勘测设计研究院，昆明650021；2.武汉大学水资源与水电工程科学国家重点实验室工程风险与防灾研究所，武汉430072）

0 引 言

滇中引水工程跨越金沙江、澜沧江、红河、南盘江四大水系，穿越横断山系高中山地貌区及滇中、滇东盆地山原区，沿线涉及多条地震带。活动断裂带突发地震带来的断层错动会导致埋地管道受力发生“S”型变形，管道可能因发生局部拉伸破坏或压缩破坏而失效，不利于整个输水管道工程的安全运行，甚至会带来严重的次生灾害［1］。研究典型跨断层埋地管道响应能够为保障活动断裂带输水管道安全提供参考依据，对输水管道工程防震减灾具有十分重要的意义。

跨断层埋地管道分析方法包括理论解析法和数值模拟法。解析方法通常以索模型或梁模型为基础进行研究，该方法难以分析管道受压和管截面大变形的情况［2］。数值模拟方法是研究埋地管道地震响应的重要手段，管-土相互作用建模方式是数值模拟方法中的重要环节，主要包括两种方法：土弹簧法和接触面法。Vazouras 等［3，4］采用非线性接触面，用壳单元和实体单元分别模拟管道和土体，对管道在走滑断层作用下进行了精细化数值模拟。目前，国内外埋地管道抗震设计指南（例如GB/T50470-2017［5］；ALA，2005［6］；PRCI，2004［7］）推荐采用土弹簧法进行跨断层埋地管道响应分析。土弹簧法相比于接触面法具有计算时间短、收敛性好及应用广泛的特点［8］。ABAQUS 软件中基于土弹簧法开发的PSI 接触单元能够合理地表征管-土相互作用的非线性行为，在埋地管道和周围土体的相互作用分析中广泛应用［9-11］。

实际工程中，断层错动引起的管道响应受到多种不确定性因素的影响［12］，比如管材力学参数的不确定性、土体参数的不确定性以及荷载不确定性等。基于蒙特卡洛的随机模拟方法能够合理地考虑上述不确定性对管道结构响应的影响。例如，樊恒等［13］研究了考虑地震水平加速度、管壁厚度、屈服强度、弹性模量和泊松比的不确定性对管道动力响应的影响。陶凯尔等［14］采用解析法分析跨断层埋地管道响应时考虑了地震烈度和管壁厚度的不确定性，采用数值模拟分析跨断层埋地管道响应时仅考虑了断层错动位移的不确定性。目前，国内外学者对跨断层埋地管道的不确定性分析研究相对较少，且一般根据跨断层埋地管道解析法进行分析，主要考虑管材参数和荷载的不确定性，没有考虑场地土体参数的不确定性，研究算例常为走滑断层。

本文以滇中引水二期工程玉溪段阿红干渠管道五段为例，研究了活动断层作用下管道的失效风险。利用ABAQUS 软件中的PSI单元模拟管-土相互作用的非线性行为，通过等效样本法［15］获得场地的土体力学参数样本值，考虑了地震位移荷载、管道材料强度以及土体不排水剪切强度的不确定性对管道动力响应的影响，运用蒙特卡洛模拟方法进行随机模拟分析，得到所研究管道的失效概率用以评估管道的安全性和可靠性。在此基础上，研究了地震震级、管壁厚度和填土深度对于管道失效概率的影响规律，为该工程及类似工程提供参考依据。

1 场地概况

滇中引水工程二期工程玉溪段输水管道布设不可避免地穿越地震活动带，管道抗震问题突出。断层活动引起的永久地面变形的发生场地相对固定，管道损坏率高，管道抗剪断问题成为管道工程中的重要问题之一。如图1所示，玉溪段阿红干渠管道同时跨越普渡河断裂带和曲江断裂带，其中阿红干渠管道五段跨越普渡河断裂带。

普渡河断裂带是滇中“南北向构造带”的主要成分之一。总体走向近南北向，倾向向东，倾角约75°～80°［16］，区域长度大于100 km，断裂性质为逆断层，左倾滑。其主干断裂大致沿普渡河延伸，总体走向近南北，全长约300 km。研究表明，普渡河断裂第四纪以来的平均水平运动速率为0.9～2.0 mm/a，平均垂直运动速率为0.35～0.47 mm/a［17］。该断裂强烈活动的最新时代为早中更新世时期，晚更新世以来，该断裂带仍有活动，但活动性明显减弱。普渡河断裂带自公元1600年记录以来，5 级以上地震17 次，包括1985年4月18日在云南禄劝-寻甸交界地区发生的Ms6.3级地震。

根据玉溪盆地钻孔探测成果图［18］，土层深度5.5 m 以上为黏土层，土层深度5.5～13.3 m 为粉质黏土。根据阿红干渠地质资料，黏土层的内摩擦角为12°～15°。阿红干渠管道五段与普渡河断裂带相交，夹角约为45°，该管段分配水头15.298 m，设计流量为5.9 m3/s，设计管径1.8 m，距离管轴的填土埋深为3.2 m，输水管道拟采用Q235C钢材。

2 确定性分析模型

2.1 位移分量

在研究跨断层埋地管道响应时，通常将其近似简化为以位移荷载控制的静力问题来分析，因此合理地估计断层位移，对研究特定场地处管道响应具有重要作用。本文仅考虑普渡河断裂带为逆断层的情况，根据地震震级，采用逆断层地震位移估算经验公式［19］，可以得出断层位移估计值。

式中：M为地震震级；MD为最大地表断层位移，m；AD为平均地表断层位移，m；σlog（MD）和σlog（AD）分别为两个回归方程的标准差。

需要注意的是，公式（1）和（2）仅适用于平行于断层走向的竖直平面上的位移，此外其只适用于M为5.4～7.4的情况。在本文有限元模拟中，仅考虑断层最大位移MD的工况，即：

式中：DFS为平行于断层走向的竖直平面上的断层位移，m。

根据位移与夹角的几何关系，可由平行于断层走向竖直平面上的断层位移推得断层的总位移［7］：

式中：Δ为断层总位移；α为断层倾角；β为断层走向与区域应力方位角之间的水平夹角（对于以逆断层为主的断层β=90°）。

图2为典型逆断层位移分量示意图。如图2（a）表示逆断层总位移与竖直方向位移的关系，图2（b）为水平面上平行于管轴方向位移与垂直于管轴方向位移的关系。已知逆断层倾角、断层与管道的夹角，可将断层总位移分解为三个相互垂直方向的位移［7，20］：

式中：κ为管道与断层之间的夹角；ΔX为平行于管道轴线方向的位移分量；ΔY为水平面上垂直于管道轴向方向的位移分量；ΔZ为竖直方向的位移分量。

2.2 管-土相互作用

管-土相互作用是影响埋地管道响应的关键因素之一。土弹簧法以弹性地基梁假设为基础，将土体对管道的作用简化为3 个方向上的土弹簧，通过土弹簧的特性曲线（图3）来表征管-土相互作用的本构关系［21，22］。

管-土相互作用的本构关系可以通过室内试验获得，也可以采用文献中的简化模型，其中美国ASCE 生命线工程技术规程中提出的管－土相互作用模型被广泛应用。本文中管道所埋地层为黏土层，表征土弹簧特性的特征值可计算如下［23］。

式中：D为管道外径；Su为土体不排水剪切强度；H为距离管道轴线的土体深度；α为与Su有关的系数；Nch和Ncv为与H/D有关的系数；Nc为与土体内摩擦角φ有关的系数；tu为管轴方向单位长度土体对管道的最大抗力；xu为yu对应的最小轴向位移；pu为水平方向单位长度土体对管道的最大抗力；yu为pu对应的最小水平方向位移；qu1和qu2分别为竖直向下和竖直向上方向单位长度土体对管道的最大抗力，zu1为qu1对应的最小竖向位移，zu2为qu2对应的最小竖向位移。

因此，当已知Su、φ、H和D时，即可获得土弹簧特性曲线。ABAQUS 中内置的PSI 单元基于土弹簧模型表征管-土相互作用。当采用PSI 单元模拟埋地管道和土体的相互作用时，可采用梁单元、管单元或弯接头单元来模拟管道，PSI 的一侧与管道共节点，另一侧代表土体表面，并用于施加边界条件，以描述地面的运动［9，10］。

2.3 建模过程

玉溪段阿红干渠管道五段管道壁厚为12 mm，外径为1.824 m，埋深H为3.2 m，φ取14°，初设土体不排水剪切强度Su为66 kPa。可估算由水压引起的管道轴向应力［24］：

式中：p为管道内部水压；δ为管壁厚度；σn为由内部水压引起的管道轴向拉应力。

初步计算可得σn=5.7 MPa，该值远远小于管材的屈服强度（235 MPa）。因此，建模时可忽略内部水压的影响。采用简化的折线弹塑性模型表征管材特性，图4所示为Q235C 钢材的理想应力-应变曲线，钢材的屈服强度为235 MPa，弹性模量为206 GPa，泊松比为0.3，拉伸强度为395 MPa。本文所建立的三维模型如图5所示，选取管道长度600 m，右侧为断层上盘，断层错动时，右盘上移。建立模型时，管道、地面和管-土相互作用分别采用PIPE31 单元、RB3D2 单元和PSI34 单元进行模拟，靠近断层的位置管单元网格划分较密，远离断层的位置管单元网格划分的较稀疏，以便获得良好的精度和计算效率。根据2.1节计算得到地震震级为6.3时的位移荷载分量，如表1所示。

表1 地震工况下的位移分量（M=6.3）Tab.1 Displacement components under earthquake of M=6.3

2.4 失效准则

基于应变的失效准则相比应力失效准则更能够发挥材料的塑性和延展性能，本文基于应变设计方法对通过活动断层的埋地管道进行抗拉伸和抗压缩验算。管道轴向容许拉伸应变和管道轴向容许压缩应变应分别按式（13）和式（14）计算［5］：

式中：［εt］F和［εc］F分别为埋地管道抗断的轴向容许拉伸应变和容许压缩应变；εtcrit和εccrit分别为管段的极限拉伸应变和压缩应变；φεt为拉伸应变承载系数，当环向应力小于屈服强度的40%时，取0.9；φεc为压缩应变承载系数，取0.6。本文中由GB/T50470-2017［5］附录D 估算Q235C 设计容许拉伸应变取值为0.91%，设计容许压缩应变取值为0.55%。

2.5 确定性结果分析

管道跟随土体发生3 个方向的位移错动，但以竖直方向的位移为主，右侧管道向上运动时，管道上部受拉，下侧受压，由于逆断层在水平方向的位移分量使得管道轴向受压，因此该管道可能发生局部严重受压。管道轴向应力云图如图6所示，断层错动发生后，管道发生“S型”形变，且大部分区域发生轴向压缩应变。图7中红色区域表示管道发生了塑性变形，可见位于断层附近的管道大部分区域已进入塑性阶段。图8为管道轴向应变包络图。如图所示管道最大拉应变εtmax为0.3%，小于［εt］F，管道最大压应变εcmax为2.2%，大于［εc］F，管道因而压缩失效。

3 不确定性分析

3.1 随机变量

跨断层埋地管道响应受到多种不确定性因素的影响。本文中考虑了管材力学参数、位移荷载以及土体参数的不确定性对管道响应的影响。钢材由于其在锻造过程中存在缺陷或者不均匀性等因素使得材料强度具有不确定性，已知Q235C 钢材屈服强度235 MPa。本文中假设钢材屈服强度服从对数正态分布，其均值为235 MPa，变异系数为10%，即标准差为23.5 MPa，其他强度参数（如极限抗拉强度与屈服强度的比值）保持不变。地震震级与最大地表断层位移之间并非满足映射关系，其关系往往由统计分析得到。因此，利用经验公式（1）计算最大地表断层位移时具有不确定性，其中log（MD）服从均值为（-1.84+0.29M），标准差为0.42的正态分布。

岩土工程勘探数据通常十分有限，以至于无法直接根据这些勘探数据获得有效的土体参数统计特征和概率分布。基于贝叶斯理论的等效样本法［15］能够将场地有限的勘探数据和先验信息结合，得到所需岩土体参数的后验分布，并采用马尔科夫蒙特卡洛模拟方法获得后验分布的大量等效样本。

本文采用玉溪盆地地区黏土层的剪切波速数据［18］计算管道所埋地层的不排水剪切强度。不排水剪切强度Su与剪切波速Vs具有以下经验关系［25］：

式中：σs为上述回归方程中不排水剪切强度Su的标准差。等效样本法已经内置于BEST 软件［15］，本文采用BEST 软件基于5 个不同深度处剪切波速Vs数据获得该土层30 000 个不排水剪切强度Su的等效样本，样本均值为66 kPa，标准差为14.5 kPa。由于BEST 采用马尔科夫链蒙特卡洛模拟产生等效样本，不同样本间存在一定相关性。为了降低样本相关性的影响，本文从30 000个样本中等间距选取1 000个样本进行管道响应分析。

在对跨断层埋地管道响应分析进行随机模拟时，根据屈服强度和最大地表断层位移的概率分布生成1 000 组随机样本，这些样本与Su的1 000 个等效样本共同作为管道分析的输入参数，进行1 000 次管道响应确定性分析，统计管道响应分析结果。

3.2 模型对比

由于材料非线性、大变形和管-土非线性接触等问题，基于有限元的随机模拟中不可避免地存在不收敛的情况。为避免或减少不收敛的情况出现，本文在随机模拟中进行了模型简化，由于本算例中基于轴向应变失效准则判断管道是否失效，轴向应变主要与轴向拉压和管道弯曲有关。如图2所示，轴向位移ΔX会直接引起管道轴向拉压变形，剪切位移ΔY和ΔZ会引起管道弯曲和轴向拉伸。管道在剪切作用下，弯曲变形使得管道两侧分别承受大小相等的拉压应变，轴向拉伸增大了受拉一侧的拉应变，而削弱了受压一侧的压应变。因此，代表剪切作用的ΔY和ΔZ对拉应变的影响较大，而对压应变的影响较小。

由公式（5）～（7）计算可得，轴向位移ΔX和水平向位移ΔY相等，且断层错动引起的竖直方向位移ΔZ为轴向位移ΔX的5.3 倍。因此，随机模拟中仅考虑施加轴向和竖直方向的位移荷载以简化模型，即将其简化为X-Z平面上的二维等效模型。在此，仅施加ΔX和ΔZ两个方向荷载的模型简称“2D”模型，而施加3个方向荷载的模型简称“3D”模型。

设管材屈服强度为235 MPa，管道壁厚为12 mm，外径为1.824 m，埋深H为3.2 m，φ取14°，Su=66 kPa。仅改变地震震级进行2D 和3D 模型的确定性分析结果对比。图9（a）和（b）所示分别为不同地震震级作用下两种模型所计算的管道最大拉应变和最大压应变。由图9可知，两种模型计算所得管道最大拉应变和最大压应变十分接近，说明了2D 模型的合理性。因此，在随机模拟中采用2D模型替代3D模型是可行的。

3.3 结果分析

图10（a）所示为根据1 000 组随机样本所得的管道最大拉伸应变分布图，其中横轴表示管道最大拉伸应变值，左侧纵轴为每一个直方区间对应的样本频率值，右侧纵轴表示样本累计频率。采用2.4 节的应变失效准则判别，可得地震震级M 为6.3时，共有744 个失效样本。图10中随机样本的轴向最大拉应变约有90%落在了0 到［εt］F之间，随机样本的轴向最大压应变有26.6%落在0 到［εc］F之间，即所有失效样本均因局部受压失效，且其中有10%的样本既发生拉伸失效又发生压缩失效。此外，在本文算例中，随机变量均值的响应（见2.5 节中的确定性分析）约为随机模拟中响应的中位数。当发生6.3 级地震时，该管道的失效概率为74.4%，其主要是由于压缩应变较大引起，即发生压缩破坏。这与确定性分析时的结果基本一致。在这种工况下，该管道不安全。

3.4 敏感性分析

除了考虑管道响应分析中的各种不确定因素之外，本节探讨地震震级、管壁厚度以及管道埋深对管道动力响应的影响。如图11所示，当管壁厚度为12 mm，埋深H为3.2 m 时，随着地震震级M的增大，管道的失效概率也急剧增大，管道的动力响应对于地震震级十分敏感。

本文研究了地震震级M为6.3，埋深H为3.2 m，管壁厚度从10 mm 到20 mm 对管道动力响应的影响。如图12所示，随着管壁厚度的增加，管道失效概率呈下降趋势，在所选管道壁厚范围内管道的可靠性对于管道壁厚较敏感。

此外，本文研究了地震震级M为6.3 时，管道填土深度从1 m 到5 m 变化对管道动力响应的影响。如图13所示，随着管道埋深的增加，未发现显著的变化趋势，与相关文献［26-28］结论基本一致。由此可见，所研究的管道的安全性能对于管道埋深不敏感。

4 结果与讨论

本文以滇中引水二期工程玉溪段管道五段为例，研究了跨断层埋地管道的安全性和可靠性问题。采用ABAQUS 内置管-土相互作用单元（PSI）建立了由地震位移控制的有限元模型，开展了地震作用下施加永久位移荷载时管道的响应分析。考虑了地震位移荷载、管材强度和土体不排水剪切强度的不确定性，分析了管道动力响应（轴向应变）的分布，基于蒙特卡洛模拟计算了特定工况下管道的失效概率。此外，本文探讨了地震震级、管壁厚度和填土深度对于所研究管道的失效概率的影响规律。主要结论如下。

（1）在本文算例中，忽略水平向位移ΔY（即将三维模型简化为二维模型）对于跨断层埋地管道响应分析影响较小，主要是由于ΔY相对于ΔZ较小，且ΔY主要引起管道弯曲，对管道的最大轴向拉应变和压应变的影响较小。因此，随机模拟可采用2D模型，以提高计算效率。

（2）敏感性分析的结果表明，跨断层埋地管道的动力响应对于地震震级和管壁厚度较为敏感，管道的失效概率随着地震震级的增加迅速增加。因此，在设计中应合理确定地震震级。管道的失效概率随着管壁厚度的增加逐渐下降，增加管壁厚度可以提高跨断层管道的安全性和可靠性。另一方面，填土埋深对管道的失效概率的影响较小。

（3）在本文算例中，管道厚度增加到20 mm 时，管道的失效概率仍大于50%，管道的安全性和可靠性较低，仅依赖增加管壁厚度并不能有效提高管道的安全性和可靠性。本文算例中可以考虑采用强度更高的管材（比如Q345 或管线钢等）或者减小断层与管道的夹角以降低轴向位移荷载从而来提高工程的安全性。 □