0 引　言

随着海洋资源的大力发展，海底管道在海底油气集输、能源供应、水资源管理、跨海交通通信等诸多领域应用越来越广泛^[1]。在海底地形变化或海水流动冲刷过程中，由于失去支撑而使管道处于悬空状态，这类管道被称为海底悬跨管道^[2]。与常规管道相比，海底悬跨管道更容易受到由海水中的氧气、温度、压力、盐分和CO₂等因素引起的应力腐蚀影响，从而增加了管道失效风险^[3]。

海底悬跨管道应力腐蚀是一种局部腐蚀，同时也是机械应力和腐蚀反应的综合和协同结果^[4]。管道的应力腐蚀微裂纹一旦形成，其扩展速度很快，且在破坏前没有明显的预兆，因此，应力腐蚀被认为是所有腐蚀类型中破坏性和危害性最大的一种^[5]。迄今为止，部分学者对管道腐蚀行为进行了研究^{[6 − 8]}。Sun等^[9]在考虑管道拉伸性能和应力腐蚀等的情况下，通过多物理场耦合技术将应力分析和电化学计算引入模型并建立了管线钢应力腐蚀三维有限元模型。Santos等^[10]通过对API X70钢试件在循环合成乙醇溶液和空气中进行慢应变速率、恒载、交错恒载和疲劳裂纹扩展试验，评价了应力腐蚀敏感性。上述管道的应力腐蚀研究未考虑管道内外温差产生的热应力对应力腐蚀的影响，无法模拟真实海底环境下的管道应力腐蚀。

海底管道在受到应力腐蚀的同时，还会受到由温度变化产生的热应力^{[11 - 12]}。温度的升高会增强电解质溶液对流和腐蚀性离子扩散的速度，从而加速应力腐蚀。然而管道表面氧化膜生长速率的增加和电解质溶液中溶解氧含量的降低又会抑制腐蚀速率。因此，海底管道的应力腐蚀与温度之间存在着复杂的关系。Raafat等^[13]借助Ansys软件，采用有限元法对钢管内、外部热应力分别进行组合分析和评估。Jeong等^[14]采用有限元分析了双层套管系统在高温条件下的热应力特性，考虑了热传导、对流和辐射条件下的热应力。孙颖等^[15]采用有限元模拟方法建立空间模型，计算了管道中产生的热应力，并对影响管道热应力的因素进行分析。上述研究仅单纯对管道热应力及其影响因素展开分析，未采用多物理场耦合模型模拟多物理场条件下的管道温度场。

综上所述，现有针对海底悬跨管道应力腐蚀研究主要集中于一系列应力腐蚀试验和有限元模拟应力腐蚀管道电-热耦合等方面。由于海底悬跨管道运行过程中存在复杂的热-电-力多物理场耦合效应，且考虑热-电-力多物理场耦合下的海底管道应力腐蚀研究处于空白。因此，为保障海洋油气的安全输送，本文利用COMSOL建立海底悬跨管道应力腐蚀三维有限元模型，并对管道的温度场进行数值模拟，通过多物理场耦合方法对管道温度场、应力场和电化学腐蚀场进行耦合分析，研究热-电-力多物理场耦合效应下海底悬跨管道的应力腐蚀行为。

1 材料与参数选择

当海底管道防腐层因机械损伤或阴极保护剥离时，管道钢直接与海水接触，从而形成一个电化学反应体系，海底管道在海水、溶解氧、CO₂以及H₂S等因素下在管外表面发生电化学腐蚀反应，电化学腐蚀参数如表1所示。

正常运行的海底输油管道，管内原油、管道以及海水共同构成了一个热力系统，因此分析海底管道的温度场时，应考虑管内原油对管道的传热和海水对管道的对流换热。本文热分析做如下假设：1）假设原油为不可压缩流体；2）与对流换热相比，忽略原油中的导热；3）忽略原油和海底土壤沿管线轴向的传热。

X70钢的弹性模量、屈服强度和热膨胀系数随温度变化所对应数值如表2所示。

原油与管道之间的传热系数由下式得出，当原油平均流速为0.5 m/s时，不考虑原油粘度随温度的变化，管内流体换热Sieder-Tate准则关系式为^[16]：

$ Nu = 1.86{\left(Re \cdot Pr\frac{d}{l}\right)^{1/3}}，$

(1)

$ R{e_1} = \frac{{d{u_1}{r_1}}}{{{m_1}}} ，$

(2)

$ P{r_1} = \frac{{{c_1}{m_1}}}{{{l_1}}}。$

(3)

管道外壁与海水的传热系数由下式得出，当海水流速为0.01 m/s时，有

$ R{e_2} = \frac{{D{u_2}{r_2}}}{{{m_2}}}，$

(4)

$ P{r_2} = \frac{{{c_2}{m_2}}}{{{l_2}}}，$

(5)

$ N{u_2} = 2 \times \frac{{0.3387Re_2^{\frac{1}{2}}Pr_2^{\frac{1}{3}}}}{{{{\left[ {1 + {{\left( {0.0468/P{r_2}} \right)}^{\frac{2}{3}}}} \right]}^{\frac{1}{4}}}}} 。$

(6)

式中：Nu为努塞尔数；Re为雷诺数；Pr为普朗特数；d为管道内径，m；l为管道长度，m；Nu=hd/λ，Re=ud/v，Pr=v/a，a=λ/ρc；h为流体与管道之间的传热系数，${\text{W}}/({{\text{m}}^2} \cdot {\text{K}})$；λ为流体导热系数，W/(m·K)；u为流体特征流速，m/s；ρ为流体密度，kg/m³；c为流体比热容，J/(kg·K)。海水和原油的物理参数见表3。

2 热-电-力多物理场耦合模型建立 2.1 几何模型

本文数值模拟采用的管道壁厚为14 mm，外径为1016 mm，有限元模拟管道长度为2 m^[9]。由于真实的腐蚀缺陷形状复杂，难以精确描述，研究人员通常认为腐蚀缺陷是规则的，因此将初始腐蚀缺陷简化为椭球体^{[17 − 19]}。如图1所示，利用COMSOL Multiphysics软件建立X70海底悬跨管道三维模型，在管道中间建立一个长240 mm、宽58 mm的椭圆形缺陷，腐蚀缺陷深度取管壁厚度的40%，由于管道的对称性，为减少计算量和计算时间，仅对管道的1/2进行仿真分析。本文在管道外部建立了一个2 m×1.5 m×2 m的计算域，用于模拟电解质溶液。为了提高计算效率，在远离腐蚀缺陷和管壁的地方使用了较粗的网格。由于腐蚀缺陷与电解质溶液边界处存在应力集中，因此腐蚀缺陷处采用精细化网格，以提高有限元模拟结果的精确度。

2.2 多物理场耦合模型

建立弹塑性应力模型，对有限元模拟X70钢的应力-应变曲线^[9]进行分段三次插值，采用小应变塑性模型对管道域进行弹塑性应力仿真，选择各向同性硬化模型，硬化函数σ_yhard定义如下：

$ \sigma\mathrm{_{yhard}}=\sigma_{\exp}\left(\varepsilon_p+\frac{\sigma_e}{E}\right)-\sigma_{ys} 。$

(7)

式中：σ_exp为X70钢应力-应变曲线；σ_ys为X70钢的屈服强度；ε_p为塑性应变；σ_e为von Mises应力；E为杨氏模量；σ_e/E为弹性应变。

阳极Tafel表达式用于模拟铁溶解反应，阴极Tafel表达式用于模拟析氢反应，局部阳极和阴极电流密度定义如下：

$ {i_a} = {i_{0,a}}\exp \left(\frac{{{\eta _a}}}{{{b_a}}}\right)，$

(8)

$ {i_c} = {i_{0,c}}\exp \left(\frac{{{\eta _c}}}{{{b_c}}}\right) ，$

(9)

$ \eta = \varphi - {\varphi _{eq}} 。$

(10)

式中：a和c分别为阳极和阴极反应；i_a和i_c为电荷转移电流密度；i_0,a和i_0,c为交换电流密度；η_a和η_c为超电势；b_a和b_c为塔菲尔斜率；φ为电极电势；φ_eq为平衡电极电势。

机械-电化学效应是指固体电极的弹性和塑性应力会影响电极反应的平衡电位。

由弹性应变引起的阳极平衡电位变化：

$ \Delta \varphi _{a,eq}^e = - \frac{{\Delta {P_m}{V_m}}}{{ZF}} 。$

(11)

由塑性应变引起的阳极平衡电位变化：

$ \Delta \varphi _{a,eq}^p = - \frac{{TR}}{{ZF}}\ln \left(\frac{{na}}{{{N_0}}}{\varepsilon _p} + 1\right) 。$

(12)

Gutman^[20]定量地表达了机械-电化学效应，在机械-电化学效应下阳极反应的平衡电势表示为：

$ \begin{split} {\varphi _{eq,a}} = &{\varphi _{eq0,a}} + \Delta \varphi _{e,a}^s = \\ &{\varphi _{eq0,a}} - \frac{{\Delta {P_m}{V_m}}}{{ZF}} - \frac{{TR}}{{ZF}}\ln \left( {\frac{{va}}{{{N_0}}}{\varepsilon _p} + 1} \right) 。\end{split} $

(13)

式中：$\Delta \varphi _{e,a}^s$为由于应力作用导致阳极平衡电位的变化；ΔP为超压，当钢的变形为弹性变形时等于单轴拉应力的1/3，当钢的变形为塑性变形时等于屈服应力的1/3；V_m为钢的摩尔体积（7.13×10⁻⁶ m³/mol），Z为电荷数（Z=2）；F为法拉第常数（96485 C/mol）；T为溶液温度；R为理想气体常数（8.314 J/mol·K）；v为方向因子（v=0.45）；α为系数（1.67×10¹¹ cm⁻²）；N₀为塑性阶段初始位错密度（管道钢为1×10⁸ cm⁻²）；ε_p为本文模型中通过应力模拟计算得到的塑性应变。

Xu等^[17]将应力对阴极反应的影响简化为应力对交换电流密度的影响：

$ i_{0,c}^s = {i_{0,c}} \times {10^{\frac{{{V_m}{\sigma _M}}}{{6{b_c}F}}}}。$

(14)

式中：$i_{0,c}^s$为受外加应力影响阴极反应的交换电流密度；i_0,c为无外加应力时在电极表面发生阴极反应的交换电流密度；σ_M为电极表面的von Mises应力，通过应力仿真计算得到。

海底管道传热的边界条件如下：

管内原油与管道内壁的对流换热边界为：

$ {\lambda _{{b}}}\frac{{\partial T}}{{\partial r}}{\mid _{r = \frac{d}{2}}} = {h_{{\text{in}}}}\Delta T 。$

(15)

管道外壁与周围环境的边界为:

$ {\lambda _{{b}}}\frac{{\partial T}}{{\partial r}}{|_{_{r = \frac{d}{2} + \delta }}} = \lambda \frac{{\partial T}}{{\partial r}}{|_{_{r = \frac{d}{2} + \delta }}} 。$

(16)

式中：λ_b为海底管道的导热系数，W/(m·K)；T为求解区域某点温度，K；r为求解区域某点至海底管道中心的距离，m；h_in为原油的对流换热系数，W/(m·K)；ΔT为原油与管道内壁温度之差，K；δ为管道壁厚，m。

进行温度影响下的应力场分析时，需要考虑温度产生的变形，其本构方程^[21]为：

$ {\sigma ^e} = {{\boldsymbol{D}}^e}({\varepsilon ^e} - {\varepsilon _0}) = {{\boldsymbol{D}}^e}\left[ {{\varepsilon ^e} - {{\boldsymbol{\alpha}} ^e}({T^e} - {T_0})} \right]。$

(17)

式中：σ^e为单元节点应力；D^e为弹性矩阵；α^e为材料的热膨胀系数矩阵；ε₀为由温度引起的初始应变；T^e为温度变化后的单元节点温度；T_e为参考温度。

海底管道温差引起的应变为：

$ {\varepsilon _{th}} = \alpha (T)({T_0} - {T_{\rm ref}}) 。$

(18)

式中：α(T)为热膨胀系数；T₀为管道的操作温度；T_ref为管道的安装温度。

间接耦合法是将第一次场分析的结果作为第二次场分析的载荷来实现2种场的耦合。本文求解流程如图2所示，首先对模型进行温度加载并施加相应的边界条件，求解得到海底悬跨管道在内外温差下的温度场分布情况；利用间接耦合法通过热膨胀模块实现热-力耦合，然后在热-力耦合的基础上，以温度影响下应力应变场的分析结果作为初始条件，利用COMSOL软件中解的继承功能，通过二次电流模块进行海底悬跨管道电化学腐蚀分析，实现了热-电-力三物理场的耦合分析。本研究采用2个步骤，利用带参数化扫描的顺序求解器进行求解。

2.3 多物理场边界条件

本文旨在分析海底悬跨管道热-电-力多物理场耦合下的应力腐蚀行为，因此所建立的模型涉及COMSOL软件的固体力模块、固体传热模块和二次电流模块，定义的边界条件如下：

1）固体力学模块。管道内壁设置边界载荷以模拟原油作用在管道内壁的内压，右端施加0.2%的拉伸应变以模拟管道正常运行过程中海土对管道的作用。为模拟海底悬跨管道受力情况，管道两侧为跨肩支撑，海床与海底管道接触部分施加固定约束。

2）固体传热模块。管道与原油和海水之间的传热通过用户自定义的方式设置传热系数和外部温度，即第三类边界条件。

3）二次电流模块。将管道设置为电极域，管外流体设置为电解质，金属材料域和电解质域的交界处为内部电极表面，内部电极表面会同时发生阳极反应和阴极反应。

2.4 模型验证

为验证本文建立的海底悬跨管道应力腐蚀三维模型的准确性，将本文有限元模拟结果与文献[9]中的计算结果进行对比，对比结果如图3所示。

表4为腐蚀缺陷深度是管壁厚度20%和40%时的von Mises应力结果对比。结果表明，由于数值分析方法及网格划分粗细不同，模拟结果存在着一定的误差，但误差较小，误差均小于5%，证明了有限元模型在预测海底管道应力腐蚀方面的有效性和可靠性。

3 有限元仿真结果 3.1 有无热应力对应力腐蚀的影响

图4和图5为有无热应力状态时缺陷处的von Mises应力和腐蚀电位云图。可知，在有热应力的情况下，缺陷中心应力保持在473 MPa，腐蚀电位越负，缺陷中心处腐蚀电位达到−0.728 V。根据腐蚀热力学原理，相比于不存在热应力的情况下，管道发生腐蚀的倾向性呈现增大的趋势。

在考虑热应力的情况下，腐蚀电位和阴极电流密度比不考虑热应力时更负，有无热应力时阳极电流密度差距更大，缺陷中心处有热应力比无热应力时大529%。因此，热应力对管道阳极电流密度的影响更加显著。具体如图6～图9所示。

3.2 管道内压对应力腐蚀的影响

图10～图13为缺陷深度为40%t时，不同管道内压下沿缺陷长度方向的von Mises应力和腐蚀电位阴阳极电流密度图。结果表明，随着管道内部压力增加，应力更加集中，使得8 MPa时阳极电流密度增加到0.093 A/m²，阴极电流密度减小到0.068 A/m²，腐蚀缺陷处的应力和阳极电流密度逐渐增加，腐蚀电位和阴极电流密度逐渐减小，且缺陷中心处变化相对于两端更加显著。管道与电解质交界处的腐蚀电位分布明显比电解质溶液中的腐蚀电位更负，腐蚀程度增强。管道内压为2 MPa和4 MPa时的von Mises应力未超过材料的屈服应力，材料属于弹性变形阶段，阳极电流密度变化很小；当内压6 MPa时，管道在缺陷处应力超过材料屈服应力，发生塑性变形。缺陷处的腐蚀电位负移，且在缺陷中心处相对于两端更加显著。缺陷中心处阳极电流密度最大，表明腐蚀底部中心区域腐蚀速度快于边缘区域。

3.3 拉伸应变对应力腐蚀的影响

图14～图17为不同拉伸应变下管道的应力、腐蚀电位及电流密度曲线图。随着拉伸应变的增大，缺陷中心位置应力变化小，而管道两端变化显著，呈现出增大的趋势。拉伸应变为0%、0.05%、0.10%、0.15%时，腐蚀电位相差不大，0.20%拉伸应变时，腐蚀电位在缺陷中心处显著减小，最小达到−0.728 V。管内温度80 ℃条件下，随拉伸应变的增大，缺陷中心处阳极电流密度变化为0.0834、0.0841、0.0842、0.0831、0.0832 A/m²，呈现先增大后减小的趋势，阴极电流密度变化为−0.0502、−0.0501、−0.0502、−0.0508、−0.0541 A/m²，呈现逐渐降低的趋势，缺陷中心处的阳极电流密度明显大于缺陷两端，中心阴极电流密度明显小于两端。

3.4 管内温度对应力腐蚀的影响

图18～图19为不同温度下的应力、腐蚀电位曲线图。由图18可知，缺陷处应力随温度的升高逐渐增大，且相同温度间隔时。40℃～60℃下应力对温度更加敏感，缺陷两端处应力相差很小。由图19可知，不同温度时，缺陷中心处腐蚀电位有一定差距，但整体上差距微乎其微。

由图20、图21可知，缺陷中心和两端的阳极电流密度随温度升高呈现稳步增长的趋势，缺陷中心变化为0.0128、0.013、0.013、0.0132 A/m²，缺陷两端变化为0.005、0.0051、0.0051、0.0052 A/m²，阴极电流密度基本保持不变。电流密度增加，物质传输和扩散速度加快，溶液中离子扩散到电极表面的阻力减小，阳极溶解离子向电极表面溶液扩散速度增加，腐蚀进程加快。

4 结　语

本文建立海底悬跨管道热-电-力多物理场耦合下的应力腐蚀有限元模型，对考虑管道内外传热影响的海底悬跨管道应力腐蚀行为进行研究，得出以下结论：

1）有热应力影响情况下，缺陷处的应力更为集中；腐蚀电位越负，管道发生腐蚀的倾向性呈现增大的趋势。

2）随着内压的增大，腐蚀电位负移，缺陷中心处相对于两端更加显著，腐蚀底部中心区域腐蚀速度快于边缘区域。

3）随着拉伸应变的增大，缺陷中心位置应力变化小，而管道两端变化显著；缺陷中心处的阳极电流密度明显大于缺陷两端，中心阴极电流密度明显小于两端。

4）随着温度的增加，管道发生腐蚀的倾向性呈现增大的趋势，电流密度增大，腐蚀进程加快。