管道热屈曲动力过程数值模拟

王立忠，王宽君，施若苇

(浙江大学 建筑工程学院，浙江 杭州 310058)

管道热屈曲动力过程数值模拟

王立忠，王宽君，施若苇

(浙江大学 建筑工程学院，浙江 杭州 310058)

运输高温高压油气的海底管道会发生整体热屈曲现象。管道热屈曲过程中可能会产生平衡状态的跃迁(snap-through)，且这样的跃迁过程必然会伴随着动力响应。管道热屈曲动力过程中侧向弹出的速度以及轴向缩进的速度对管土相互作用参数的取值有很大影响，然而关于管道热屈曲动力过程的研究却很少。本文给出了数值模拟过程中管道系统阻尼值和升温速率的确定方法，研究了管道初始几何缺陷以及海床参数对管道热屈曲动力过程的影响。

管道热屈曲；动力过程；阻尼值；升温速率；初始几何缺陷；海床参数

高温高压管道在内部温度、压力荷载以及海床轴向约束的作用下会产生巨大的轴力，当轴力达到一定值时管道就会发生类似于欧拉梁失稳的现象，如图1所示，管道工程中把这样的失稳现象称为管道整体热屈曲(Global pipeline buckling)。如图2所示，对于理想的平直管道，一定的温度荷载对应着两个平衡状态(曲线1)：①B点所示平衡状态，其变形幅值较小且轴力没有充分释放，属于不稳定平衡状态，②C点所示平衡状态，其变形幅值很大且轴力得到充分释放，属于稳定的平衡状态[1]。在微小扰动下，处于B点平衡状态的管道会发生失稳现象，突然形变(snap- through)至C点的平衡状态。对于初始几何缺陷较小的情况，管道的失稳过程也存在突然形变的现象，如曲线2中A点至D点；当初始几何缺陷增大后，突变现象消失，温度荷载与屈曲状态呈如曲线3所示的单调函数。借鉴钢结构稳定理论[2]中对轴心受压构件失稳的定义，管道热屈曲模式可以分为如下两种：①跳跃型(曲线2)：管道在失稳过程中会发生形变上不连续的跃迁，与轴心受压构件的跃越失稳模式[2]相似；②分岔型(曲线3)：管道的荷载形变曲线没有极值点，但平衡状态会发生明显改变(平衡状态的分岔)，与轴心受压构件的稳定分岔失稳模式[2]类似。需要注意的是稳定分岔失稳的定义是针对无缺陷轴心受压构件，而本文中定义的分岔型失稳是针对有缺陷的管道。在实际工程设计中需要避免管道发生跳跃型屈曲失稳，因为跳跃型屈曲失稳会使管道受到额外的动力荷载而更容易发生破坏。

图2 管道热屈曲温度幅值曲线Fig. 2 Temperature rise versus bucking amplitude

目前关于管道热屈曲的数值模拟大多采用静力方法(Static- implicit)，这对于分岔型屈曲是适用的，但对于跳跃型屈曲，静力算法忽略掉了管道失稳过程中的有关状态量，譬如管道弹出的速度。同时，采用静力算法计算跳跃型屈曲失稳时，收敛性很差，特别当管道初始缺陷较小、管道变形的突变量较大时，静力算法基本上无法收敛。赵天奉[3]采用了改进的Riks算法，可以得到管道后屈曲状态，但也无法模拟管道热屈曲过程中的动力响应。Liu等[4]基于ABAQUS采用显式动力算法(Dynamic- explicit)分析了管道整体热屈曲过程，并指出显式动力算法可以模拟管道热屈曲过程中的动力响应。然而，显式算法需要很小的时间增量步，其时间增量步取决于模型的最高固有频率而与荷载的类型和持续时间无关，通常模拟需要取一万到一百万个增量步。

管道升温至临界屈曲温度的过程显然不可能在几秒钟完成，根据管道离热源(油井)的距离通常需要几分钟到几个小时不等[5]。管道的屈曲失稳过程相对于管道整个升温过程是短暂的，大多数时间内管道处于静力状态，只有在屈曲失稳的瞬间存在动力响应。如果采用显式算法，由于受最小时间增量步的限制，需要的时间增量步数量巨大。然而，采用隐式动力算法(Dynamic- implicit)没有最小时间增量步的限制，对于动力响应以外的时间段可以提高时间增量步大小，提高了计算效率。因此，本文将基于ABAQUS采用隐式动力算法模拟管道热屈曲过程中的动力响应。

1 数值模型建立

1.1模型基本参数

管道发生屈曲时管道与海床相互作用的长度较小，摩擦力对两端几乎没有影响，两端采用固支约束不会对管道屈曲产生影响。所以本文数值模拟采用的管道长度为1 000 m，管道末端为固支约束。管道中部设有正弦函数形状的侧向初始几何缺陷，缺陷长度取100 m，缺陷幅值为vom，具体管道的有限元模型如图3所示。管道的几何和力学参数如表1所示，其中A是管道截面积，I是截面转动惯量，E是管道材料的弹性模量，α是管道材料的膨胀系数，w是管道的有效重度，D是管道直径。

图3 有限元模型Fig. 3 The finite element model

表1 管道参数Tab. 1 Pipeline parameters

管道采用B21H单元，即线性的二维梁单元，管道单元长度均为1 m，这样设置可以统一管土相互作用参数。采用非线性土弹簧模拟管土侧向和轴向相互作用，管土相互作用模型如图4所示。轴向管土相互作用采用双线性模型，滑移距离s取0.003 m。侧向管土相互作用分别采用双线性和三线性模型并进行对比分析。根据White和Cheuk[6]的建议，滑移距离s(s1)和s2分别为0.1D和0.25D左右，取s=s1=0.06 m，s2=0.24 m。

图4 管土相互作用模型Fig. 4 Soil- pipe interaction model

1.2阻尼值和升温速率的选取

(1)阻尼值的选取

在ABQUS隐式动力算法中，采用瑞利阻尼：

式中:α和β是由用户自定义的参数，M为系统质量，K为刚度矩阵。瑞利阻尼属于数值阻尼，没有物理意义。在数值试算中，取β=0；初始几何缺陷为vom=1m；管土相互作用为双线性模型，侧向摩擦系数取μy=1，轴向摩擦系数取μx=0.5；升温速率为0.01°C/s。

由于目前尚无管道热屈曲过程中阻尼值选取的相关资料，DNV- RP- F105规范给出了管道涡激振动时阻尼比的取值方法[7]：

图5 不同阻尼值下的管道中点侧向位移与温度的关系Fig. 5 Lateral displacement versus temperature of pipeline midpoint for different damping values

阻尼比是当前阻尼值与临界阻尼值的比值。任何一个振动系统，当阻尼增加到一定程度时，物体的运动是非周期性的，当阻力使振动物体刚能不作周期性振动而又能最快地回到平衡位置的情况，物体振动连一次都不能完成时的阻尼值称为“临界阻尼”。方程(2)中，ζT是系统总的的阻尼比，ζstr是结构阻尼比，ζsoil为土的阻尼比，ζh为水的阻尼比。ζstr根据管道配重的不同，取值范围约为为0.005～0.02，ζsoil根据土的性质不同，取值范围约为0.01～0.03，通常ζh取0。对于涡激振动，管道悬跨段不与土体接触，而管道侧向屈曲过程中整个管道都与土体接触，所以土体阻尼比需要相应增大。另一方面，瑞利阻尼模型对于大阻尼系统(阻尼大约超过临界阻尼的10%)是不可靠的[8]。因此，本文取阻力比为ζT=0.05。

为了确定自定义参数α的值，图5给出了不同α情况下管道中点位移vo与温度T的关系曲线。当α=0.5和5时，管道失稳时会发生动力抖动，之后逐渐接近静力计算结果；当α=20时，管道不存在动力抖动且变形值比静力计算结果小，属于过阻尼；当α=10时，管道动力计算结果与静力计算结果十分接近，此时的阻尼值为临界阻尼。根据阻尼比的定义可以得到本文模型的阻尼值α=0.05×10 =0.5。

(2)升温速率

管道系统的热传递受很多因素影响，譬如石油的温度，油气与管道之间的热传导系数，管道绝热层厚度，管道所处的环境温度等。根据Alves 等[5]的算例，热源处管道从0°C到50°C的平均升温速率约为0.12°C/s，典型的管道升温速率在0.01°C/s的数量级。

在数值试算中，管土侧向相互作用取双线性模型，侧向摩擦系数μy=1，轴向摩擦系数μx=0.5；阻尼值取α=0.5；初始缺陷分别取vom=1 m；管道的升温速率分别取0.01，0.1，1和10°C/s。

图7 升温速率为0.1和0.01°C/s时管道中点侧向速度对比Fig. 7 Lateral velocity versus temperature of pipeline midpoint for 0.1 and 0.01°C/s

从图6中可以看出，对于0.1°C/s和1°C/s的升温速率，管道位移曲线基本吻合，且管道屈曲后侧向弹出的速度会被系统阻尼消耗而逐渐减小至零。然而对于10°C/s的升温速率，管道位移与其它两种升温速率相比已有明显变化，且管道屈曲后侧向速度会一直存在，数值计算不收敛。1°C/s的升温速率可以得到收敛的计算结果，但其侧向速度峰值比0.1°C/s时偏大约40%。

图7对比了0.01°C/s和0.1°C/s升温速率情况下管道中点的侧向速度。升温速率为0.01°C/s时，管道的峰值速度约为1.61 m/s；而升温速率为0.1°C/s时，管道的峰值速度约为1.69 m/s，两者相差不超过5%。然而采用0.01°C/s升温速率所耗费的计算时间是0.1°C/s的10倍，因此，本文数值计算中采用升温速率为0.1°C/s是合适的。

2 参数分析

2.1初始几何缺陷的影响

管土侧向相互作用采用双线性模型，侧向摩擦系数取μy=1，轴向摩擦系数取μx=0.5。对于vom=0.5～10 m的初始几何缺陷，管道中点轴力po、侧向位移vyo与温度T的关系如图8所示。初始几何缺陷vom=0.5和1 m时，管道发生跳跃型屈曲失稳，其它初始缺陷情况下管道发生分岔型屈曲失稳。管道屈曲失稳前处于轴向受压状态，轴力随温度升高而增长，位移基本为零；管道屈曲失稳后处于压弯状态，位移vyo迅速增长，轴力得到释放。管道的临界屈曲温度Tcr随其初始缺陷vom的变大而逐渐减小；但初始缺陷持续增大时，管道临界屈曲温度Tcr不会持续减小而是稳定在23°C左右，如图9所示。

图8 不同初始缺陷下管道中点轴力和位移的变化Fig. 8 The influence of different initial geometric imperfection: (a) axial force, (b) displacement

图9 管道屈曲失稳临界温度与初始缺陷的关系 Fig. 9 The critical temperature versus initial geometric imperfection

管道中点侧向速度uyo在升温过程中的变化如图10所示。初始缺陷vom=2 m时，管道在整个升温过程中不发生动力响应，侧向速度基本为零；当vom=1 m时，管道发生失稳时的最大侧向速度约为1.8 m/s,vom=0.5 m时管道最大侧向速度约为5 m/s，可见管道热屈曲失稳时的侧向速度随初始缺陷的减小而增大。升温过程中，离管道中点50 m处的轴向运动速度ux50如图11所示。管道屈曲失稳时，轴向速度的数值比侧向速度小一个数量级，轴向速度也随初始缺陷的减小而增大。可以看出，管道初始缺陷越小管道发生跳跃型屈曲失稳时的动力响应越明显。

管道失稳时的侧向与轴向速度大小对管土相互作用参数的取值有着很大的影响。土体的不排水抗剪强度su跟应变率有关，且管道运动速度影响海床土体中孔隙水压力的产生和消散，从而影响土体抗力的大小[9]。通过基于动力算法的有限元计算能得到管道屈曲时的侧向和轴向速度，从而为管土相互作用试验和数值模拟提供参数，确定管土相互作用加载速率的选取范围。

图10 管道中点侧向速度Fig. 10 Lateral velocity of pipeline midpoint

图11 离管道中点50 m处的轴向速度Fig. 11 Axial velocity of the point 50 meters away from pipe midpoint

随着初始缺陷的减小，管道失稳时的突变量会逐渐增大。当vom<1 m后，静力算法通常只能计算到屈曲突变之前，随后即无法收敛，不能得到屈曲突变的过程和后屈曲状态。因为此时管道屈曲失稳前后所对应的两个平衡状态的形变量相差过大，数值计算无法从前屈曲的平衡状态迭代到后屈曲的平衡状态。然而，无论初始缺陷值取多小动力算法都不会产生不收敛的问题，且动力算法可以得到管道热屈曲过程中侧向弹出和轴向缩进的速度。因此，动力算法更适用于对管道热屈曲的数值模拟，特别是初始缺陷较小且管道发生跳跃型失稳的情况。

2.2管土相互作用参数影响

管土侧向和轴向相互作用都采用双线性模型，轴向摩擦系数取定值μx=0.5，对侧向摩擦系数μy进行参数分析。从图12和13中可以看出，海床侧向抗力减小会导致管道临界屈曲温度减小，且管道屈曲失稳时的动力响应会减弱；相同温度情况下，海床侧向抗力减小会导致管道屈曲幅值增大。以初始缺陷vom=1 m为例，侧向摩擦系数为1.5和1时，管道会发生跳跃型失稳；但当侧向摩擦系数减小到0.5时，管道发生分岔型失稳。当vom=2 m，侧向摩擦系数为1和1.5时，管道发生分岔型失稳；但当侧向摩擦系数达到3时管道会发生跳跃型失稳。因此，管道屈曲失稳的类型不仅与初始缺陷有关还与海床侧向抗力有关，海床侧向抗力越大管道越容易发生跳跃型失稳。

图12 初始缺陷vom=1 m时管道中点侧向位移和速度Fig. 12 Lateral displacement and velocity of pipeline midpoint for vom=1 m

图13 初始缺陷vom=2 m时管道中点侧向位移和速度Fig. 13 Lateral displacement and velocity of pipeline midpoint for vom=2 m

管土侧向和轴向相互作用都采用双线性模型，侧向摩擦系数取定值μy=1，对轴向摩擦系数μx进行参数分析。从图14和图15中可以看出，轴向摩擦系数对管道的临界屈曲温度基本没有影响，但管道的后屈曲变形幅值随μx的提高而减小。管道发生热屈曲以后，管道平直段会有向屈曲段挤进的趋势(feed- in)，轴向摩擦力约束了管道平直段的轴向挤进从而抑制了屈曲变形。从图14中可以看出，对于跳跃型屈曲，轴向摩擦系数越大，管道发生屈曲失稳时的动力响应越明显，管道侧向弹出的速度越大。然而，提高轴向摩擦力并不会改变管道屈曲失稳的模式，对于vom=2 m的情况(图15)，轴向摩擦系数从0.5增长到1时，管道的屈曲失稳类型一直为分岔型失稳。

图14 初始缺陷vom=1 m时管道中点侧向位移和速度Fig. 14 Lateral displacement and velocity of pipeline midpoint for vom=1 m

图15 初始缺陷vom=2 m时管道中点侧向位移和速度Fig. 15 Lateral displacement and velocity of pipeline midpoint for vom=2 m

海床对管道的轴向摩擦力影响管道后屈曲状态，但对管道前屈曲状态和临界屈曲温度没有影响。这是因为管道屈曲前处于轴向压缩状态，轴向摩擦力不发挥作用；一旦管道发生屈曲失稳形成屈曲段，管道平直段会产生向屈曲段挤进的趋势，此时轴向摩擦力开始发挥作用。轴向摩擦力越大对管道轴向挤进的约束越大，从而管道屈曲段变形越小。

下面进行三线性管土相互作用模型的分析，相关参数见表2。在三线性管土相互作用模型中，Hbrk=μy1w表示管道破土时海床对管道的最大抗力(break out resistance)，Hres=μy2w表示管道破土后海床对管道的残余抗力(residual resistance)；其中μy1和μy2是海床与管道之间的摩擦系数，w是管道有效重度。三线性管土作用模型通常用来模拟“轻管”的侧向管土相互作用过程。所谓“轻管”就是Vmax/V> 2.5的管道[10]，其中Vmax指管道整个生命周期中海床对其最大抗力，V是目前海床对管道的竖向抗力。Vmax/V越大会导致μy2/μy1越大。

表2 三线性管土相互作用参数Tab. 2 Parameters of trilinear response for soil- pipe interaction

从图16中可以看出，μy1越大管道临界屈曲温度越高，管道屈曲失稳时的动力响应也会越明显，管道侧向弹出的速度越大。然而，管道屈曲后的幅值几乎不受μy1的影响。因为此时管道位移已经相对较大，管土相互作用已经进入三段线的最后一段，此时海床侧向抗力为Hres=μy2w，所以管道后屈曲变形是由μy2控制的。

图16 初始缺陷vom=1 m时管道中点侧向位移和速度Fig. 16 Lateral displacement and velocity of pipeline midpoint for vom=1 m

从图17可以看出，μy1影响管道的屈曲失稳类型。当μy1从1变为1.2时，管道从分岔型失稳转变成跳跃性失稳。随着μy1增大，管道的动力响应越来越明显，管道弹出的速度也越来越大。相对于双线性的管土相互作用模型，三线性管土相互作用模型更准确地模拟了“轻管”的管土相互作用模式。由于三线性管土相互作用存在峰值抗力和抗力跌落的过程，管道更容易发生跳跃型失稳。

3 结 语

本文基于ABAQUS有限元软件采用隐式动力算法(Dynamic- implicit)，对管道热屈曲过程进行了数值模拟，并给出了阻尼取值和升温速率的合理取值。详细结论如下：

1)隐式动力算法能模拟管道热屈曲过程中的动力响应，得到管道屈曲失稳过程中的侧向弹出和轴向缩进速度，从而为管土相互作用模型提供参数。

2)管道的热屈曲失稳过程分为两种：跳跃型失稳和分岔型失稳。管道发生跳跃型失稳时会产生动力响应，发生分岔型失稳时没有动力响应。

3)管道初始缺陷越小越容易发生跳跃型失稳，反之则发生分岔型失稳。管道的临界屈曲温度随其初始缺陷的增大而逐渐减小；但初始缺陷持续增大，管道临界屈曲温度不会持续减小而是稳定在某一温度。管道初始缺陷越小，其屈曲失稳过程中的动力响应越明显，管道侧向弹出的速度和轴向缩进的速度越快。

4)海床对管道的侧向抗力是除初始缺陷以外影响管道热屈曲类型的另一个重要因素。在双线性管土相互作用模型中，海床对管道侧向抗力越大，管道越可能发生跳跃型失稳，失稳时的动力响应也越大。在三线性管土相互作用模型中，管道破土时海床最大抗力Hbrk对管道热屈曲的动力响应影响很大，Hbrk越大管道越容易发生跳跃型失稳，且失稳时的动力响应也越大。而管道屈曲失稳后的变形由海床对管道的残余抗力Hres控制。

[1] HOBBS R E. In- service buckling of heated pipelines[J]. Journal of Transportation Engineering, 1984, 110(2): 175- 189.

[2] 陈骥. 钢结构稳定[M]. 理论与设计. 北京: 科学出版社, 2003. (CHEN Ji. Stability of steel structures- theory and design[J]. Beijing: Science Press, 2003. (in Chinese))

[3] 赵天奉. 高温海底管道温度应力计算与屈曲模拟研究[D]. 大连理工大学, 2008. (ZHAO Tianfeng. Research on thermal stress and buckling of HT submarine pipelines[D]. Dalian: Dalian University of Technology, 2008. (in Chinese))

[4] LIU R, XIONG H, WU X, et al. Numerical studies on global buckling of subsea pipelines[J]. Ocean Engineering, 2014, 78: 62- 72.

[5] ALVES L, SOUSA J R M, ELLWANGER G B. Transient thermal effects and walking in submarine pipelines[J]. International Journal of Modeling and Simulation for the Petroleum Industry, 2012, 6(1):37- 44.

[6] WHITE D J, CHEUK C Y. Modelling the soil resistance on seabed pipelines during large cycles of lateral movement[J]. Marine Structures, 2008, 21(1): 59- 79.

[7] VERITAS N. Free spanning pipelines[M]. Det Norske Veritas, 2002.

[8] SYSTèMES D. Abaqus analysis user’s manual[M]. ABAQUS Inc., Dassault Systèmes, France, 2007.

[9] RANDOLPH M F, WHITE D J, YAN Y. Modelling the axial soil resistance on deep- water pipelines[J]. Géotechnique, 2012, 62(9): 837- 846.

[10] RANDOLPH M, GOURVENEC M R S. Offshore geotechnical engineering[M]. CRC Press, 2011.

Numerical simulation for the dynamic process of pipeline buckling

WANG Lizhong, WANG Kuanjun, SHI Ruowei

(College of Civil Engineering and Architecture, Zhejiang University, Hangzhou 310058, China)

Subsea pipelines are used for transporting high temperature and high pressure oil and gas, which may lead to global buckling of the pipeline. In the process of thermal buckling, the snap- through of the pipeline always exists, which is regarded as a dynamic process. The velocities of the lateral and axial in the dynamic process have great effect on the evaluation of parameters in pipeline- soil interaction analysis. But the study of thermal buckling dynamic process is rarely seen. In this paper, the dynamic process of pipeline thermal buckling is simulated by ABAQUS, and a reasonable evaluation of damping value and heating rate is provided. The variation of the initial geometric imperfection and seabed parameters are also taken into consideration to analyse the dynamic process of pipeline thermal buckling.

subsea pipeline; thermal buckling dynamic process; damping value; heating rate; initial geometric imperfection; seabed parameters

TE832

A

10.16483/j.issn.1005- 9865.2015.05.010

1005- 9865(2015)05- 073- 08

2014- 07- 18

国家杰出青年科学基金资助项目(51325901)；国家自然科学基金面上资助项目(51279176)

王立忠(1969- )，男，教授,主要研究方向为海洋岩土工程。E- mail: wlzzju@163.com