双拱初始缺陷海底管线水平向整体屈曲数值模拟分析

刘 润, 刘文彬, 洪兆徽, 闫澍旺

(天津大学 水利工程仿真与安全国家重点实验室, 天津 300072)

随着全球经济的高速发展, 各国对石油天然气的需求量日益增加, 海底管线具有连续输送、运输量大、管理方便、稳定安全等一系列优点, 成为海上运输油气的主要手段, 据统计, 发达国家管线油气运输方式的输油量约占油气运输总量的2/3之多。管线运输在原油、天然气的生产、精炼、储存及到使用的全过程起到了重要的作用[1]。为了避免在运输过程中石蜡分馏产生的固化影响并降低运输难度, 油气必须在高温高压下运输。温度应力与材料泊松效应的共同作用使管线中产生较大的附加应力, 由于海床土体对管线的约束作用, 管线无法自由变形释放应力, 应力的不断累积使管线发生竖直向或水平向的整体屈曲变形, 对于不埋或浅埋的管线主要发生水平向屈曲[2]。由于管线内承载着大量油气, 一旦发生屈曲甚至屈服破坏, 就会导致油气泄露, 不仅造成资源的浪费, 而且对海洋生物及人类的生存环境构成威胁。因此研究海底管线在温度应力下的整体屈曲特性具有重要的现实意义。

由于受到研究对象尺度和工作条件的限制, 管线整体屈曲的现场实验和室内实验开展较少, 理论分析和数值模拟为主要的研究手段。1984年, Hobbs[3]首先将Ker[4]关于铁轨的研究成果应用于海底管线整体屈曲的研究中, 提出了理想海底管线整体屈曲的4种经典模态, 如图1所示。将图中距离管线中点较近的、水平向变形较大的变形段称为主弯段, 而将第三阶和第四阶中主弯段以外的较小变形段称为次弯段,图中v11、v21、v31、v41分别为第一阶、第二阶第三阶和第四阶屈曲模态下主弯段上点的水平向位移,v32、v42分别为第三阶和第四阶屈曲模态下次弯段上点的水平向位移,L1、L2、L3、L4分别为第一阶、第二阶、第三阶和第四阶模态下主弯段的屈曲长度,L31和L41分别为第三阶和第四阶模态下主弯段和次弯段长度之和,Ls为滑移段管线长度。

图1 整体屈曲模态Fig.1 Modes of global buckling

实际工程中海底管线在制造和铺设过程中, 会由于人为因素或海床的不平坦产生局部的弯曲变形,称为几何初始缺陷, 几何初始缺陷导致海底管线不能平直地裸置于海底。1986年, Taylor等[5]考虑了土抗力在管线变形过程中的变化, 首先提出了海底管线具有的单拱和双拱几何初始缺陷的缺陷形态, 推导了缺陷管线发生低阶模态水平向整体屈曲的解析解。1996年, Taylor等[6]等通过对海底管线竖直向的整体屈曲进行理论和实验研究, 提出了三种不同模式的缺陷形式, 并给出了不同缺陷形式下管线临界屈曲荷载计算公式。1997年, James[7]在考虑几何初始缺陷及残余应力影响的条件下, 采用简化计算模型研究了管线的隆起屈曲, 给出了管线隆起屈曲的临界荷载计算公式。2004年, Junes等[8]建立了埋设管线的非线性有限元模型, 对非线性地基土约束作用下含几何初始缺陷管线的竖直向整体屈曲变形进行了数值模拟。2004年, Edgard等[9]开展了关于管线整体屈曲的经典实验, 但在实验过程中没有考虑地基土的作用。2005年, David等[10]基于非线性有限元分析, 揭示了整体屈曲形成的机理并指出管线内部轴力、管线的初始几何缺陷和土体对管线的约束力是控制管线整体屈曲形态的主要因素。2007年, Yu等[11]运用2D和3D数值分析方法研究了海底管线水平向整体屈曲中管-土相互作用过程及管线直径和重量对屈曲形态的影响。近些年, 一些学者开始对含几何初始缺陷的海底管线进行数值模拟, 关于几何初始缺陷的引入方法, 一般都是在建立计算模型的过程中, 通过设定缺陷幅值与波长的大小, 将几何初始缺陷绘制于计算模型中, Yu等[11]、Newson等[12]、Liu等[13]都是通过该方法将海底管线的几何初始缺陷引入到计算模型中开展了对整体屈曲的数值模拟研究。2010年, Antunes等[14]基于Hobbs的研究成果推导了含浮力模块管线的整体屈曲解析解, 揭示了含浮力模块的海底管线发生整体屈曲的变形特性。2014年, Liu等[15]运用了4种有限元方法模拟海底管线水平向的整体屈曲。2014年, Liu等[16]运用能量法推导了理想管线发生高阶模态整体屈曲的解析解,分析了温度应力、管线截面刚度、地基土特性以及干扰力对管线整体屈曲的影响。以上学者运用解析解和数值模拟的方法研究了涉及海底管线的整体屈曲问题, 但对于管线整体屈曲过程的变形发展研究较少。

本文假定海底管线具有双拱反对称初始缺陷,采用模态分析法将最可能出现的缺陷形态引入数值分析模型中, 针对管线在高温高压作用下发生整体屈曲的动态变形特征, 运用显式动力数值模拟方法研究了管线整体屈曲过程中水平向变形与轴向变形随温度的变化规律, 分析了管线整体屈曲过程中屈曲段与滑动段轴力的变化过程。将数值模拟结果与经典解析解和 Edgard室内试验结果相对比, 验证了本文方法的可靠性。

1 整体屈曲的数值模拟

1.1 显式动力法

管线在铺设和制造过程中由于人为或自然因素, 不可避免地产生局部的弯曲现象, 称为几何初始缺陷。大量研究结果表明, 含几何初始缺陷的海底管线较理想平直管线更易发生整体屈曲。由于实际工程中海底条件极为复杂, 发生整体屈曲的管线较长, 现场试验和室内取样较难开展, 因此研究海底管线整体屈曲的手段主要以理论分析和数值模拟为主。

温度荷载作用下, 海底管线发生整体屈曲的过程实际上是细长材料的失稳过程, 图2为Taylor等[5]分析得到的管线整体屈曲过程中屈曲幅值(vm)随温度荷载(T)的变化曲线。由图2可知, 当管线初始缺陷幅值与波长的比值v0mL0较小时, 管线的失稳过程存在屈曲幅值动态跳转的现象, 即相同温度荷载下, 管线的屈曲幅值由曲线中的E点跳转至F点, 当初始缺陷幅值与波长的比值v0mL0增加后, 动态跳转现象消失。

图2 管线屈曲幅值与温度关系曲线Fig.2 The change curve of vm vs T

20世纪以前, 对温度荷载作用下海底管线整体屈曲的数值模拟大多选用静力算法, 这种算法对于无动态跳转的管线失稳是适用的。但当管线初始缺陷形态较小, 屈曲过程发生动态跳转时, 一方面, 静力算法的收敛性大大降低, 基本不能收敛, 虽然Risk算法能够得到管线后屈曲的状态, 但也无法得到管线屈曲过程的动态响应。另一方面, 对于有动态跳转的整体屈曲, 静力算法计算不到管线屈曲过程中的状态量, 如管线变形速率。因此本文基于ABAQUS有限元软件, 选用动力算法计算管线的整体屈曲过程。

ABAQUS中动力算法包括显式算法和隐式算法2种, 隐式算法采用隐式时间积分, 其时间步长由算法自动计算, 在屈曲发生瞬间投入大量计算资源,能够精确得到管线整体屈曲发生瞬间的应力状态,但进入后屈曲状态后其计算精度降低。显式算法的计算精度与计算总时间有关, 管线整体屈曲发生瞬间和后屈曲过程中都可以保持相同的计算精度。本文着重研究管线从整体屈曲发生时到后屈曲过程的应力状态, 因此, 需要保持相同的计算精度, 因此选用显式动力算法对海底管线的整体屈曲进行数值模拟。

1.2 模态分析法

管线几何初始缺陷的引入, 可采用以下 2种方法: 一是在建立管线的分析模型时, 直接按照缺陷幅值与波长的数值将缺陷绘入几何模型; 二是通过2次建模, 基于模态分析法引入缺陷。本文采用了后一种方法, 建立的管线几何缺陷模型如图3所示。

图3 模态分析法引入的几何初始缺陷Fig.3 Imperfection introduced by buckle analysis method

图3采用的模态分析法是一种基于概率理论的分析方法, 可以有效引入完全光滑的初始缺陷, 避免应力集中现象, 其模拟过程主要分为2个阶段: 首先, 根据管线的几何参数运用 Buckle分析法对管线进行屈曲模态分析, 得到管线最可能产生的几何缺陷形态, 通过“Nodefile”语句输出管线的节点位移;而后, 重新建立用于显示动力分析的计算模型, 运用“Imperfection, file=job-name, step=1”语句将管线节点位移以比例因子的形式导入到计算模型中, 完成管线几何初始缺陷的引入。

2 算例分析

2.1 算例概述

渤海某油气管线, 材料规格为 API 5L X65, 设计温差88℃, 设计内压4.65 MPa, 详细参数见表1。

表1 管线及地基土参数Tab.1 Parameters of pipeline and soil

为了便于计算, 管线运营过程中所受的内压一般转化为温差, 二者转化公式为:

公式(1)中,A为管线截面面积,p为管线所受的内压,r为管线的半径,E为管线钢材的弹性模量,α为管线钢材的线膨胀系数,t为管线的壁厚,ν为泊松比。根据上式可确定管线所受4.65 MPa的内压等效温差为5℃, 因此管线所受总温差为93℃。

2.2 计算模型

按照表1中的参数建立三维有限元模型, 管线的模型长度由屈曲段长度和滑动段长度组成[17], 根据本文选用的算例数据, 将管线的模型长度设为2 000 m, 地基土体的模型长度在管线模型的基础上两端各增加10 m, 为2 020 m, 地基土体模型的宽度为10 m, 近似等于30倍管径。按照前文所述的模态分析法将双拱反对称几何初始缺陷引入到计算模型中, 模型中运用接触单元模拟管土间的相互作用,其接触类型为面-节点接触, 接触法则为运动接触法,相对位移的跟踪遵循有限滑移法。运用罚函数模拟管土间轴向与水平向的摩擦特性, “hard”接触模拟管线与土体的法向接触。在整体屈曲的计算模型中分2个分析步模拟真实条件下含双拱反对称几何初始缺陷管线的整体屈曲特性: 第一分析步为缺陷管线在重力作用下产生自沉与地基土体接触; 第二分析步为缺陷管线受温度荷载作用发生整体屈曲。得到的管线有限元计算模型如图4所示。

图4为模态分析法引入双拱几何初始缺陷管线的计算模型, 其中管线选用pipe单元, 遵循Simpson积分规则。地基土体采用六面体 C3D8R实体单元,应力应变关系采用Mohr-Coulomb模型, 土体深度为1 m。管线和土体在非独立网格中进行装配。其中管线共2 000个网格; 地基土体共40 000个网格。由于管线较长, 摩擦影响长度延伸不到管线的两端, 所以管线两端边界条件设为自由, 土体四周边界条件为水平向及轴向变形固定, 底部边界为完全固定。

图4 数值模拟中的计算模型Fig.4 Calculation mode in numerical solution

2.3 分析结果

2.3.1 变形分析

数值模拟得到含双拱缺陷的海底管线在温度荷载作用下发生水平向整体屈曲的发展过程, 如图5所示。随着温度荷载增加程度不同, 管线水平屈曲变形呈现不同的变化特点, 当温度荷载由 0℃升高到9.3℃时, 管线的屈曲变形在几何初始缺陷的基础上进一步的发展, 此过程中管线的水平位移增大但屈曲波长减小, 发生水平“凸起”轴向“缩进”的变形现象, 其变形形态属于图1中的二阶屈曲模态。随着温度荷载的增加, 尤其是温差由9.3℃升高到93℃的过程中, 管线的水平位移及屈曲波长的增长速率增加,屈曲段两端产生反方向的小弯曲段, 其变形形态属于图1中的四阶模态, 此温度范围内管线发生较大的变形, 且其形态也发生了改变。2个阶段管线变形速率的变化为先慢后快。

图5中 A点为管线屈曲变形的转折点, 管线发生整体屈曲时该点的屈曲幅值最大, 应力集中, 最容易产生屈服破坏。绘制A点的Mises应力随屈曲幅值的变化曲线如图6所示, 并与材料的屈服应力相比较, 判断管线屈曲过程中是否发生屈服破坏。由图6可知, A点的Mises应力随管线屈曲幅值的增加逐渐增大, 其最大值小于材料的屈服应力448 MPa。说明海底管线在温度荷载作用下, 发生了屈曲大变形, 但没有发生屈服破坏。

图5 管线整体屈曲过程中水平位移的发展过程Fig.5 The developing processes of lateral deformation in lateral buckling

图6 A点Mises应力随屈曲幅值的变化曲线Fig.6 Stress change curve of turning point A

数值模拟得到含双拱反对称几何初始缺陷的管线在温度荷载作用下发生整体屈曲时轴向变形的发展过程如图7所示。对于理想的平直管线, 温度作用下管线不发生整体屈曲, 但由于热膨胀的作用, 管线会产生以中点为对称点的轴向扩张变形, 其变形特征为在管线中点处变形为0, 沿管线向两端逐渐增长, 在两端点处达到最大值, 称为最终扩张量。管线的最终扩张量随着温度荷载的增加而不断增大, 当温度达到设计温差93℃时, 管线最终扩张量为0.94 m。

图7 管线整体屈曲过程中轴向位移的发展过程Fig.7 The developing processes of axial deformation in lateral buckling

对于具有双拱几何初始缺陷的管线, 随温度荷载的增长管线在初始缺陷的基础上发生水平向变形,根据变形协调原理, 管线在发生水平向整体屈曲时,轴向会产生相应的滑动以补偿水平向的变形, 在数值模拟中表现为位于屈曲段两端的管段产生向中点滑动的轴向缩进变形, 其变形特点是在管线中点处改变运动方向, 沿管线两端变形量逐渐减小最终减小为0。对于轴向缩进变形量为0的点在工程上称为虚拟锚固点, 相邻 2个虚拟锚固点之间的距离称为虚拟锚长度如图8所示, 虚拟锚长度是研究管线整体屈曲的最小单元, 它包含了管线屈曲变形后的屈曲段长度和滑动段长度, 是数值模拟中衡量模型长度是否能模拟真实管线的重要标准。随着温度荷载的增加, 滑动段产生的轴向缩进变形逐渐增长, 当温度达到设计温差 93℃时, 管线的轴向缩进变形量为0.2 m, 虚拟锚长度以外的两端管线仍然发生轴向扩张变形, 扩张量为 0.76 m, 小于平直管线的最终扩张量。

数值模拟得到含双拱几何初始缺陷的管线在温度荷载作用下发生整体屈曲的最终变形图如图8所示。图8描述了数值模拟管线水平向整体屈曲后的最终变形形态, 包括水平向屈曲变形与轴向变形。水平向屈曲变形关于管线中点反对称, 其形态为图1中的第四阶屈曲模态, 中间有 2个较大的反对称弯曲变形两侧有反方向的小弯曲变形; 轴向变形主要分为 2部分, 一部分为虚拟锚段管线为了补偿水平屈曲变形而产生的以管线中点为中心的轴向缩进变形, 另一部分为虚拟锚两端的管线因受热膨胀而产生的轴向扩张变形, 轴向缩进变形与轴向扩张变形的交点 M、N是既不无水平变形也无轴向变形的屈曲不动点, 即前文所述的虚拟锚固点, 相邻2个虚拟锚固点之间的距离为虚拟锚, 虚拟锚长度主要由发生水平整体屈曲的屈曲段长度与两端发生轴向缩进变形的滑动段长度组成。

图8 数值法中管线整体屈曲的最终形态Fig.8 The final deformation of pipeline in numerical solution

2.3.2 轴力分析

初始缺陷是影响海底管线整体屈曲的重要因素,主要由缺陷幅值vm与波长L0两部分组成。由于初始缺陷幅值与波长的比值v0mL0能够综合反映两部分的特性, 先前的学者们(如 Taylor[5])以v0mL0的大小来反映初始缺陷的凹凸程度, 初始缺陷幅值vm越大,波长L0越小,v0mL0越大, 初始缺陷形态越凸出,温度荷载作用下管线越容易发生整体屈曲。本文基于 Taylor[5]的研究成果选取初始缺陷幅值与波长的比值大小为v0mL0= 0.007、v0mL0= 0.01、v0mL0=0.015、v0mL0= 0.02和v0mL0= 0.025, 分析不同几何初始缺陷形态下, 发生整体屈曲的海底管线所受的轴力。

在数值模拟中以产生最大水平屈曲变形的点A(参见图4)作为研究对象, 通过改变管线几何初始缺陷的v0mL0分析缺陷形态对屈曲管段轴力的影响,得到的计算结果如图9所示。

由图9可知, 屈曲段的轴力随屈曲幅值的增加呈现先增大后减小的变化规律, 说明屈曲段的轴力在温度荷载作用下先累积而后随屈曲幅值的增加而释放。当v0mL0由0.007增加到0.025时, 管线起始发生整体屈曲的屈曲段轴力Pm受几何初始缺陷形态v0mL0的影响较大, 由0.31 MN减小为0.18 MN。

图9 数值模拟v0 m L0对屈曲段轴力的影响Fig.9 The influence ofv0 m L0 on axial force in buckle section

以轴向缩进变形为0的虚拟锚固点M(参见图8)作为研究对象, 分析管线几何初始缺陷形态的v0mL0对产生轴向缩进变形的滑动段轴力的影响,得到的计算结果如图10所示。滑动管段的轴力随屈曲幅值的增加而增长最终趋于稳定, 这是由于滑动管段在产生轴向缩进变形的过程中受到地基土的摩阻力作用使得管段内的轴力发生累积, 其累积程度与滑动管段所受的摩阻力大小有关。管线发生轴向缩进变形的滑动段轴力Ph随几何初始缺陷形态v0mL0的增加逐渐减小。

图10 数值模拟v0 m L0对滑动段轴力的影响Fig.10 The influence ofv0 m L0 on axial force in slip section

计算得到含双拱几何初始缺陷的管线发生整体屈曲后管线内轴力的分布曲线, 如图11所示。理想管线在温度荷载作用下, 受热产生的轴力为377.88 kN。相同条件下, 含双拱几何初始缺陷的管线由于屈曲变形的发生使管线内的轴力得到一定的释放。通过图中不同v0mL0管线发生整体屈曲的结果对比可知,v0mL0越大管线发生整体屈曲时释放的轴力越大,残存于管线内壁上的轴力越小。

图11 整体屈曲后管线轴向压力的分布Fig.11 Distribution of axial force after lateral buckling

2.4 与经典解析解的对比

为了验证本文分析方法的准确性并说明与经典解析解方法的差异, 将含双拱缺陷管线发生水平向整体屈曲的四阶模态数值模拟结果同经典的二阶模态整体屈曲的解析解绘制于同一个图中, 如图12所示。解析解与数值模拟两种方法得到缺陷管线的屈曲段轴力均呈先增大后减小最终趋于稳定的变化规律, 但由于经典解析解的屈曲变形形态为二阶模态,数值模拟的屈曲变形形态为四阶模态, 二阶模态的屈曲段轴力一般小于四阶模态, 因此本文中的数值模拟结果大于经典解析解结果, 但两者的结果较相近, 最大相差不到32%。

图12 数值模拟结果同经典解析解的比较Fig.12 Comparison between numerical solution and classical analytical solution

2.5 与室内试验结果的对比

2004年Edgard开展了关于海底管线整体屈曲的室内模型实验, 其具体试验参数如表2所示。

表2 Edgard的试验参数Tab.2 Parameters of Edgard’s test

运用本文数值方法对Edgard的室内模型实验进行模拟, 将数值模拟的结果与试验结果绘制于图13中。由图13可知, 数值模拟方法得到的管线整体屈曲形态与 Edgard的实验结果相近, 管线中点处的位移基本相等, 数值模拟得到管线屈曲位移略小于实验观测值, 但二者最大相差不到 30%, 验证了本文采用的数值模拟方法的可靠性。

图13 数值模拟结果同室内试验结果的对比Fig.13 Comparison between numerical solution and test results

3 结论

本文将模态分析法与显示动力法相结合, 实现了ABAQUS中对含双拱反对称几何缺陷管线整体屈曲过程的数值模拟, 将模拟结果同经典解析解和Edgard的室内实验相对比。现得到结论如下:

1) 整体屈曲的发生是一个由低阶向高阶发展的过程, 具有双拱反对称几何初始缺陷的海底管线在温度荷载作用下首先发生二阶模态的整体屈曲进而发展为四阶模态, 此过程中管线屈曲变形速率逐渐增大。

2) 高温高压海底不埋管线的整体屈曲变形包括水平向变形和轴向缩进。水平向变形产生于屈曲管段, 是在初始缺陷的基础上进一步发展形成的, 当管线具有反对称初始缺陷时, 其水平向屈曲形态亦关于中点反对称; 轴向缩进变形产生于滑动管段,变形在中点处改变方向, 变形量随距中点距离的增加而减小。屈曲管段与滑动管段长度之和为两虚拟锚间的距离, 是管线屈曲分析的最小单元。

3) 管线整体屈曲的水平向变形释放了屈曲管段内的轴力, 轴向缩进变形累积了滑动管段内的轴力,轴力的释放量与累积量受管线几何初始缺陷形态v0mL0的影响,v0mL0越大管线起始发生整体屈曲所需的屈曲段轴力Pm越小, 屈曲变形过程中屈曲段释放的轴力越大, 滑动段累积的轴力越小。

[1] 周延东, 刘日柱.我国海底管道的发展现状与前景[J].中国海上油气(工程), 1998, 10(4): 1-5.

[2] 刘润, 郭林坪, 刘文彬, 等.温度应力下海底管线水平屈曲分析[J].科技论文在线, 2011, 6(5): 368-373.

[3] Hobbs R E.In-service Buckling of Heated Pipelines [J].Journal of Transportation Engineering, 1984, 110(2):175-189.

[4] Kerr A D.On the stability of the railroad track in the vertical plane[J].Rail Int., 1978, 5: 131-142.

[5] Taylor N, Gan A B.Submarine Pipeline Buckling Imperfection Studies[J].Thin-Walled Structures, 1986,4: 294-323.

[6] Taylor N, Tran V.Experimental and theoretical studies in subsea pipeline buckling [J].Marine Structures,1996, 9 (2): 211- 257.

[7] James G A C.A simplified model of upheaval thermal buckling of subsea pipelines [J].Thin Walled Structures,1997, 29: 59-78.

[8] Junes A V, Jose F R, Cora M.Buried Pipe Modeling with Initial Imperfections [J].Journal of Pressure Vessel Technology, 2004, 126(2): 250-257.

[9] Edgard P J, Renato S R.Experimental tests and numerical simulation in a reduced model in a pipeline with zig-zag geometry – a parametric study [C]//International Pipeline Conference.International Pipeline Conference.Calgary, 2004.

[10] David B, Malcolm C, Michael C, et al.The safe design of hot on-bottom pipeline with lateral buckling using the design guideline developed by the SAFEBUCK joint industry project[C]//Deep Offshore Technology Conference.Deep Offshore Technology Conference.Vitoria, 2005: 1-26.

[11] Yu S K, Konuk I.Continuum FE modeling of lateral buckling: Study of soil effects[C]// OMAE.Proceedings of the International Conference on Offshore Mechanics and Arctic Engineering – OMAE.2007: 347-354.

[12] Newson T A, Deljoui P.Finite Element Modeling of Upheaval Buckling of Buried Offshore Pipelines in Clayey Soils[C]// ASCE.Soil and Rock Behavior and modeling., 2006: 351-358.

[13] Liu R, Wang W G, Yan S W.Finite element analysis on thermal upheaval buckling of submarine burial pipelines with initial imperfection[J].Journal of Central South University.2013, 20(1): 236-245.

[14] Antunes B R, Solano R F, Vaz M A.Analytical formulation of distributed buoyancy sections to control lateral buckling of subsea pipelines[C]// Proceedings of the ASME 29th International Conference on Ocean, Offshore and Arctic Engineering.Shanghai, 2010, 5: 669-677.

[15] Liu R, Xiong H, Wu X L et al.Numerical studies on global buckling of subsea pipelines[J].Ocean Engineering, 2014, 78: 62-72.

[16] Liu R, Liu W B, Yan S W.Global lateral buckling analysis of idealized subsea pipelines[J].Journal of Central South University, 2014, 21(1): 416-427.

[17] Liu R, Liu W B, Hong Z H, et al.Triggering methods research on pipeline global buckling basing on reasonable length of VAS finite element model[J].Journal of pipeline Systems Engineering and Practice,(under review).