海底管道在弯矩和水压作用下的屈曲压溃研究

余建星 ，李牧之 ，余 杨 ，韩梦雪 ，李 杨 ，于佳晖

(1. 天津大学水利工程仿真与安全国家重点实验室，天津 300072；2. 高新船舶与深海开发装备协同创新中心，上海 200240)

海底管道在安装铺设和服役过程中所处的环境复杂，“J”型、“S”型和卷管铺设法[1]都会使管道受到弯矩和外压的联合作用，与单独载荷的作用情况相比，联合载荷的作用会改变管道的屈曲破坏机理，导致其极限承载力的降低，因而对海底管道在联合载荷作用下的屈曲机理进行深入研究有助于海底管道的抗屈曲设计．Johns等[2]通过试验探究了管道在弯曲和外部静水压力组合作用下的破坏．Ozkan等[3]基于Mises屈曲准则和理想弹塑性假设，在管道的极限应力状态下，推导了管道在内压、轴力和弯矩组合作用下的相互作用方程．王慧平等[4]在此基础上采用双线性理想弹塑性模型，考虑了管道截面塑性区的椭圆化变形以及管道材料的各向异性，推导了管道在内压、轴力和弯矩联合作用下的极限弯矩承载力解析解．除了上述方法，能量法在求解管道屈曲破坏的理论模型中也有广泛运用．Corona等[5-6]通过研究发现管道初始曲率所产生的残余应力和管道的初始椭圆度会造成非对称屈曲的发生，并以虚功原理为基础提出了一套理论公式．张子骞等[7]从旋转壳体的一般几何方程出发，基于形变理论和能量理论，运用里兹法求解了具有初始曲率的圆柱壳体在弯曲作用下的塑性失稳临界弯矩．陈飞宇等[8]针对含有初始缺陷的海底管道，在不考虑加载路径对结构失效的影响下，建立了管道在弯矩、轴力和外部静水压力的联合作用下的屈曲压溃理论模型，研究表明，在复杂载荷的作用下，各因素对管道压溃压力影响的非线性特征明显．Kyriakides等[9-10]通过试验和理论分析的方法，将弯矩和外压作用下的无初始缺陷管道简化为二维模型，从而运用非线性环理论得出先施加外压再施加弯矩的加载路径更易使管道进入屈曲状态．

理论计算在研究复杂载荷下管道的屈曲破坏形式及极限承载力的计算中起到了重要作用，但经过试验结果验证了的有限元分析方法也是解决海底管道复杂载荷问题的重要手段．Bai等[11]采用数值模拟方法，研究了径厚比在 10～40之间的管道在外压和弯矩共同作用下的屈曲响应，得到了在外压和弯矩联合作用下的管道极限承载性能．C-FER公司采用深水试验舱进行了管道全尺寸试验[12]，得到了管道在固定外压下的极限弯矩承载力，并运用有限元方法进行建模，得到的管道变形形式和应力状态与试验结果基本吻合，证实了有限元方法计算复杂载荷作用下管道极限承载力结果的正确性，但为了模拟试验装置设置的模型边界条件无法实现自由弯曲．Ghazijahani等[13]运用非线性有限元方法分析单一加载路径下管道在弯矩和外压作用下的极限承载力和变形形式，并针对大径厚比管道进行了全尺寸试验．Gong等[14-16]针对纯外压和复杂载荷作用下的管道，运用理论和有限元方法对由于初始缺陷形式造成的管道非对称屈曲进行了研究，并得出几何和材料因素对管道承载能力的影响．杨政龙等[17]采用有限元软件建立了管道与海床的数值模型，模拟了外压与冲击载荷联合作用下结构动态响应过程．余建星等[18]针对含内部冲刷腐蚀的弯管，采用有限元方法模拟了弯管的 3种屈曲模态，并针对冲刷腐蚀损伤对弯管屈曲压力的影响因素作了敏感性分析．何璇等[19]选取含球形凹陷和椭球形凹陷的圆管作为研究对象，在有限元软件中建立管道的对称模型，分析了在外压和弯矩的组合作用下屈曲载荷临界值相对于纯载荷作用下的偏离程度，但在管道轴向上施加的约束也导致模型无法实现自由弯曲．

可以看出，目前针对弯矩和外压作用下管道的理论模型并没有充分考虑不同加载路径以及管道初始缺陷的问题，有限元模型也因为约束的方式而未能模拟出更贴近实际工况和理论模型的自由弯曲变形，并且分析重点主要针对管道在各种敏感性因素影响下的极限承载力，对于弯矩和外压联合作用下的管道压溃破坏机理以及出现危险加载路径的原因还有待进一步探索．

本文在对管道模型做了一定简化的基础上，根据非线性环理论列出几何方程，选择流动理论作为本构关系，采用虚功原理建立能量平衡方程，并假设含初始缺陷管道的位移函数，从而在不同的加载路径下建立了弯矩和外压联合作用下含初始椭圆度缺陷的管道压溃破坏的理论模型．同时在有限元软件ABAQUS中，为了实现管道模型的自由弯曲变形和不同加载路径的计算，设置合理的边界条件和加载方式建立模型，在不同的分析步下实现对不同加载路径下管道承载能力的计算，并与理论模型进行对比，验证了理论模型的准确性．之后在有限元模型中对复杂载荷作用下的管道进行分析计算，以缺陷截面作为研究对象，追踪截面的椭圆度变化情况，进而探究弯矩和外压荷载对管道的联合作用及不同加载路径下管道的屈曲破坏机理．

1 理论模型

管道在外压的作用下，会在缺陷处发生局部压溃，管道压溃压力的确定，与失稳之后的后屈曲阶段无关[20]．因此，基于 Kirchhoff-Love假定和平断面假定条件建立二维圆环模型，忽略圆环上与环向应力相比较小的径向应力，并且认为在管道轴向上各截面变形一致．由于采用了薄壁假定，管道的径厚比越大越接近实际情况[21]，因此只考虑径厚比大于20的管道.

根据以上假定，建立管道在弯矩和外压联合作用下的屈曲压溃理论模型，如图 1所示．图中：p表示外压；M表示弯矩；κ表示轴向曲率；R表示管壁中性层半径；z表示中性层法向坐标；v和 w分别为中性层处质点的环向位移和径向位移；η和ζ分别为中性层处某质点的横纵坐标，该点与纵轴的夹角为θ．

图1 弯压联合作用下的管道理论模型Fig.1 Theoretical model of the pipe under combined action of bending moment and pressure

1.1 几何方程

根据薄壳理论，管道模型的环向应变θε由薄膜应变εθ,0和弯曲应变两个部分构成，由非线性环理论[10]可得

式中：κθ为环向曲率；为非圆管道初始缺陷引起的应变．位移v和w的一次导数和二次导数都表示对θ求导．

因为假设管道轴向变形一致，所以管道截面任意一点的轴向应变xε为

式中εx,0为管道弯曲时中性面处的轴向应变．

1.2 本构关系

将加载过程划分为若干个增量步，采用流动理论建立应力增量dσij与应变增量dijε之间的关系[22]．同时，为了确定强化模量h的值，采用Ramberg-Osgood模型进行拟合．应变增量dijε可以划分为弹性应变增量dεij,e和塑性应变增量dεij,p两部分之和，即

式中：E为弹性模量；μ为泊松比；sij为应力偏张量；σy为名义屈服应力；n为硬化系数；σe为等效应力．

近似地，第i+1个增量步下的应力、应变全量即可由第i个增量步下的应力、应变全量以及当前增量步下的应力、应变增量之和确定，即

式中i为非负整数．

1.3 虚功方程

管道在弯曲和外部静水压力联合作用下，由虚功原理可得平衡方程为

对于弯矩载荷，曲率在加载过程中是预先给定的，因此弯矩对管道所做虚功为零[8]，所以有

1.4 求解过程

在MATLAB计算软件中用程序化的语言表示离散后的几何方程、本构关系和由能量法列出的非线性方程组．采用 Gauss-Legrendre法对圆环截面划分的积分单元进行求积计算，采用 Newton-Raphson法求解非线性方程组．

为了满足位移边界条件，选择三角级数对环向和径向位移v、w进行离散，即

式中vn和 wn为假设的未知参数．同时假定管道结构的初始缺陷w0是沿轴向一致的椭圆度缺陷，

在程序中将各偏导数表示为位移函数w和v对θ的导数的形式，进而可以根据总位能增量的变分δΠ˙对假设为三角级数形式的位移函数的参数求偏导数，得到了 2N+2个方程，对应εx,0、wn(n =0,1,… , N)、vn(n = 1 ,2,… , N )，共2N+2个未知数．

模型求解流程如图 2所示．当 i＝0时，εij、σij及h的初始值均为 0，即认为结构上不存在残余应变或残余应力．以弯矩—外压加载路径为例，为了得到弯矩作用下的管道承压能力，采用曲率(κ)-外压(p)的加载方式，即通过对管道施加曲率进而产生弯矩，在程序中先施加轴向曲率κ到预设值，再令外压从0开始逐步加载．

图2 理论模型求解流程Fig.2 Flow chart of the theoretical model solution procedure

鉴于含初始椭圆缺陷的管道在弯矩和外压作用下的对称性，采用半圆截面作为分析对象，其中位移级数 N取 9，径向和环向高斯积分阶数分别取 5和32．在每个加载步下，各未知数的初值列阵X设定为0，采用牛顿迭代法求解非线性方程组，即将 X代入上述几何方程和本构方程，得到由各个方程的值组成的列阵F，并对各个方程依次求关于各个未知数的偏导数，形成雅可比矩阵 J，进而得到更新的解列阵X′，即

根据|X-X′|的二范数是否小于某一个较小的容许误差δ来判定解是否满足收敛条件，即当|X-X′|的二范数小于 δ时，认为方程的解已收敛于 X′，否则令X＝X′，继续迭代过程．在程序调试中，为保证数值方法的计算结果能收敛于精确解，应当设置合理的收敛步数，保证结果的收敛性．若超过一定迭代步数仍未达到收敛，令 i＝i+1，并将载荷控制增量减半，重新进行加载计算，直到载荷增量小于预设值，则认为此时为极限加载状态，可以求得极限载荷．输出管道压溃破坏时半圆环截面构型及 Mises应力分布如图 3所示．

图3 管截面计算构型及应力分布Fig.3 Calculation configuration of the pipe section and stress distribution

2 有限元模型

2.1 建模过程

对于复杂载荷作用下的管道在ABAQUS中的模拟来说，最重要的是确定合理的边界约束，要求既不影响管道的变形，又能限制模型的刚体位移．因此采用 C3D8R单元建立 1/4管道模型，一端为含初始椭圆度缺陷的半椭圆面，一端为无缺陷的半圆截面．在椭圆截面和管道的轴向截面设置相应的对称约束来反映模型的对称情况并限制模型在这两个方向的刚体位移，此外，在椭圆截面处采用接地弹簧的方法限制第3个方向的刚体位移，管道模型如图4所示．模型的几何与材料属性如表 1所示，其中 D为管道外径，L为管道轴向长度．

图4 管道模型示意Fig.4 Diagram of the pipe model

表1 管道模型的材料与几何属性Tab.1 Material and geometric properties of pipe

为了与理论模型的加载方式保持一致，设置参考点与半圆端面相耦合，通过对参考点施加转角α从而控制管道的轴向曲率，转角与曲率的对应关系为κ=2α/ L ．为了模拟不同的加载路径，分别在两个分析步下设置相应的载荷进行加载．在 κ—p加载路径下，设置第 1个分析步 Step1为静力通用计算方法，从而对管道施加一定的转角载荷达到预设值，管道在弯矩载荷的作用下发生了弯曲变形，产生了一定的曲率．在第2个分析步Step2中选择弧长法进行分析，在保持弯矩载荷稳定的基础上，对管道外表面施加静水压力．管道模型在保持曲率不变的情况下，当外压载荷达到一定数值时，在缺陷处会发生局部屈曲现象．计算完成后可以得到管道的变形形式和极限承载力，在后处理模块中通过对称的方法得到完整管道模型的变形形式，如图 5所示．可以看出，在此约束和加载条件下，管道可以发生自由弯曲变形．

同理，可以通过在 Step1和 Step2中，分别施加外压和转角载荷实现对p—κ加载路径的有限元模拟.

图5 κ—p加载路径下管道的变形形式Fig.5 Deformation forms of pipe under κ—p loading path

2.2 模型验证

利用理论和有限元计算方法，对不同加载路径下管道的承载能力进行分析，在 κ—p加载路径下对理论模型和有限元模型施加相同的曲率并计算压溃压力，在 p—κ加载路径下对两种模型施加相同的外压并计算管道出现压溃破坏时的曲率，得到结果如图 6所示，其中横坐标为管道的轴向曲率 κ，纵坐标为在此弯矩载荷下管道达到压溃破坏的压力 p．结果表明：①无论在何种加载路径下，理论模型与有限元模型对于管道极限承载力的计算结果相差较小，可以证明此理论模型及采用的数值求解方法对管道模型的适用性；②p—κ加载路径较为危险，在相同的曲率下，此路径下管道达到压溃破坏时所需外压载荷较小.

图6 不同加载路径下理论和有限元模型计算结果Fig.6 Results of the theoretical and finite element model under different loading paths

3 复杂载荷的屈曲破坏机理分析

3.1 弯矩对压溃压力的影响

管道在外压作用下会在缺陷最大处发生局部压溃，为了探究弯矩对压溃压力的影响，追踪缺陷截面的椭圆度在弯矩作用下的变化情况．在 ABAQUS中对 κ—p加载路径下模型进行计算，在 Step1计算完成后，提取纯曲率 κ作用下管道缺陷截面的椭圆度Δ，之后在 Step2中计算在此曲率载荷作用下管道的压溃压力 pc．同时，以椭圆度Δ作为管道的初始椭圆度缺陷，用弧长法计算纯外压作用下管道在此椭圆度缺陷下的压溃压力 pco，与弯矩载荷作用下的压溃压力pc相对比，结果如表2所示．

表2 复杂载荷和纯外压作用下的管道压溃压力Tab.2 Collapse pressure of the pipe under complex load and pure pressure

由表2可以看出：①弯矩载荷会加大含初始椭圆度缺陷管道的截面椭圆度，又因为压溃压力会随着椭圆度的增大而减小，所以弯矩会降低管道的承压能力；②弯矩作用下管道的压溃压力 pc小于纯外压作用下直管在此弯矩载荷产生的椭圆度缺陷Δ下的压溃压力 pco，因此认为弯矩载荷对管道承压能力的削弱作用不仅体现在椭圆度的变化上，尤其是曲率越大，造成的椭圆度缺陷越大，pc与 pco两者之间相差最大可达到 30%．此外，在应力分析中可以看出，虽然在管道接近压溃时外压会在管道厚度方向上产生相反的轴向应力，对弯矩作用产生的拉或压应力有一定程度的削弱，但随着载荷的增大，Mises应力始终增大，直到管道压溃．

综上可见，弯矩载荷对压溃压力的削弱体现在Mises应力和截面椭圆度增大两个方面，但在管道曲率载荷较小时，椭圆度的影响起主导作用．

3.2 不同加载路径变形截面分析

为了探究不同加载路径下管道承载能力变化的原因，运用有限元计算方法对管道模型在Step1中施加一定的外压 p1，在加载外压达到预设值后，在Step2对管道施加转角载荷直至破坏，记录管道破坏时的曲率 κ，得到管道缺陷截面压溃瞬间的截面椭圆度为 Δ1．之后，为了探究不同加载路径下的管道屈曲状态，建立κ—p路径下的管道模型，在Step1中对管道施加 p—κ路径下管道破坏时的曲率 κ，在管道弯曲达到曲率κ后，在Step2中对模型施加外压载荷直至破坏，记录管道在弯曲作用下的压溃压力为 p2，管道缺陷截面压溃瞬间的截面椭圆度为 Δ2，绘制出不同加载路径下的计算结果，如图7所示．

通过计算不同加载路径下管道在压溃时刻缺陷截面的椭圆度可以看出，管道在弯矩载荷作用达到相同的曲率时，p—κ加载路径下的管道在压溃瞬间的截面椭圆度较大，因此压溃瞬间的椭圆度是产生危险加载路径的原因，并且在复杂加载中，随着曲率的增大，不同加载路径下所产生的压溃瞬间的椭圆度差距也增大．

图7 管道压溃瞬间的截面椭圆度和压力Fig.7 Section ovality and pressure at the pipe collapsing instant

4 结 论

(1) 本文建立的理论模型和有限元模型可以准确计算出不同加载顺序下弯矩和外压对管道的极限承载力的影响，相比于弯矩—外压的加载路径，管道在外压—弯矩的加载路径下承载能力较差，此路径较为危险．

(2) 在有限元模型中，以压溃瞬间管道缺陷截面的椭圆度为分析指标，通过对压溃瞬间管道缺陷截面椭圆度的分析，得出压溃瞬间的截面椭圆度变化是不同加载路径下管道承载能力存在差异的原因，管道在外压—弯矩的加载路径下截面椭圆度较大．

(3) 弯矩载荷对压溃压力的削弱体现在 Mises应力和截面椭圆度增大两个方面，在管道弯矩载荷较小时，椭圆度的影响起主导作用．