高温高压高产气井油管柱三维非线性流致振动模型

郭晓强，李潇，柳军，方达科，黄亮，魏安超，王国荣; GUO Xiao-qiang; LI Xiao; LIU Jun; FANG Da-ke; HUANG Liang; WEI An-chao; WANG Guo-rong

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

高温高压高产气井油管柱三维非线性流致振动模型 PDF

- ORCID：
郭晓强 ¹
✉
- ORCID：
李潇 ¹
- ORCID：
柳军 ¹
✉
- ORCID：
方达科 ²
- ORCID：
黄亮 ²
- ORCID：
魏安超 ²
- ORCID：
王国荣 ¹

1. 西南石油大学机电工程学院，四川成都 610500； 2. 中海石油（中国）有限公司湛江分公司，广东湛江 524057

中图分类号： O322； TE931

最近更新：2022-08-30

DOI：10.16385/j.cnki.issn.1004-4523.2022.04.016

摘要

针对高温高压高产（三高）气井流体诱发油管柱振动失效问题，采用微元法、能量法结合哈密顿原理，考虑了油管自重、井眼轨迹变化、井筒温度压力变化、油套管接触碰撞等因素，建立了三高气井油管柱三维非线性流致振动模型。采用拉格朗日和三次埃尔米特函数离散控制方程，并使用增量形式的Newmark‑β和Newton‑Raphson联合迭代求解离散后的控制方程。根据中国南海三高气井参数，采用相似原理，设计了油管柱流致振动模拟试验装置，测得油管柱振动响应数据，与理论模型计算结果对比，验证了振动模型的正确性。采用实例井参数，分析了油管柱振动特性，结果表明，上部油管三维空间横向振动明显，而中下部油管长时间紧靠于套管壁上发生振动，导致下部油管易发生纵向和横向的磨损，同时油管纵向振动出现混沌现象，且中下部位置油管混沌现象最明显，因此，现场设计人员需重点关注下部油管柱的安全性。

关键词

非线性振动; 管柱力学; 混沌现象; 高温高压高产气井; 振动模拟试验

引言

随着浅层油气资源的日益减少，中国钻井、完井工艺不断向高压、高温和复杂的深部地层方向发展，以高产的开采方式才能够满足当前的需求。与常规气井油管柱相比，三高气井（地层孔隙压力大于69 MPa，温度高于150 ^oC，产量高于120万方/天^［

1］）油管柱将面临更大的风险，主要表现为管内高速流体诱发油管的非周期性剧烈振动，将增加其轴向载荷，引起屈曲变形，更加容易导致油管发生疲劳、磨损等失效^{［参考文献 2

百度学术}2］。一旦油管柱发生破坏，将被迫井下作业，甚至导致井筒报废，造成重大的经济损失和环境污染。

国内外学者针对内流引起的管道振动问题开展了相应的研究，并取得了一些成果。早期学者针对内流作用下的输流管道振动问题开展了初步研究^［

3］，通过试验验证了内流对输流管道的动力响应有显著影响^{［参考文献 4

百度学术}4］。随后，诸多学者针对管道振动模型开展了研究，建立了流体作用力的计算方法^{［参考文献 5

百度学术}5］、管柱纵向振动^{［参考文献 6

百度学术}6］、横向振动^{［参考文献 7

百度学术}7］以及流固耦合振动模型^{［参考文献 8

百度学术}8］，以上模型主要集中在单个方向的振动，且主要适用于输流管道，无法有效分析油气井开采管柱的振动响应。近年来，桑松等^{［参考文献 9

百度学术}9］发现细长结构的纵横向耦合效应明显，不可忽略；邢誉峰等^{［参考文献 10

百度学术}10］针对结构的纵横向振动问题，建立了结构的纵横向耦合振动方程；Liu等^{［参考文献 11‑12}11‑12］考虑几何非线性因素，建立了海洋立管的纵横向耦合振动模型；针对油套管非线性接触碰撞问题，文献［13‑15］主要集中于钻柱‑套管、油管‑套管、换热管束‑支承板之间的接触载荷研究，提出了管柱屈曲后接触力的静态计算方法和换热管束动态碰撞力计算方法。同时在前期研究中^{［参考文献 16

百度学术}16］，建立了高产气井油管柱双重非线性振动模型，重点考虑了管内流体和油套管接触碰撞作用，且属于二维振动模型，能够有效分析由于井斜角变化引起的流致振动响应，忽略了高温高压对管柱振动的影响，以及方位角变化引起的流致振动响应，具有一定的局限性。对于高温高压气井管柱安全性的研究，文献［17‑18］主要考虑井筒温度和压力变化引起的附加应力，提出了管柱“四大效应”变形理论。但忽略了油管柱内部流致振动引起的附加应力。

随着中国南海高温高压高产曲井（定向井和水平井）的大力开发，面临着作业困难、管柱失效率高等问题，亟需开展三高气井油管柱三维非线性流致振动模型研究，以指导现场管柱的设计及作业参数的优化配置，保障中国南海高温高压高产气井的安全开采。所建模型包含的主要创新为：①模型考虑油管柱的三维振动，即纵向振动和两个横向振动。主要原因是目前中国南海东方、陵水以及乐东等区块主要采用的深井（井深大于4500 m）、曲井（定向井、大斜度井以及水平井），在这些复杂的气井中，油管柱属于大长径比结构，易发生纵横向耦合效应，为此不能忽略管内气体引起的纵横向耦合振动；同时其井眼轨迹复杂，方位角和井斜角变化较大，容易影响管柱两个横向方向的振动，二维振动模型无法满足其精度要求。②模型考虑高温高压对管内气体的影响，中国南海高温高压气田井底温度达151 ^oC，压力达65 MPa，易导致管内气体密度和流速随井深方向发生变化，进而影响管柱的振动，如考虑管内气体为恒定的流速将与现场不相符，导致分析的振动响应也不准确。

1 油管柱三维非线性流致振动模型

1.1 振动控制方程

三高气井以定向井和水平井为主（如图1（a）所示），井斜角和方位角随井深发生变化，当管内高速气田流动时，将产生冲击载荷诱发油管纵横向三维振动。为此，采用能量法结合哈密顿原理建立油管微段的三维纵横向耦合振动控制方程，由于油管微段很短，视为直段，建立以竖直向下为z轴的正方向，水平向右方向为x轴的正方向，y轴满足右手定则的直角坐标系（如图1（b）所示）。在建模之前，作了如下假设：

图1 三高气井油管柱井身结构和微段坐标系

Fig.1 Wellbore structure and micro segment coordinate system of tubing string in 3H gas wells

（1）油管挂将海水段与地层段分开，重点建立地层段油管柱流致振动模型；

（2）将管内高速气井视为单一气体，重点考虑气体因井筒温度压力变化而引起其密度、流速的变化，忽略气体在井筒流动过程中的相态变化。

基于以上假设，把油管简化为均匀的Rayleigh梁，考虑纵横耦合，其格林应变为^［

11］

\{\begin{array}{l} ε_{x x} = \frac{\partial u_{1}}{\partial x}, ε_{y y} = \frac{\partial u_{2}}{\partial y}, ε_{z z} = \frac{\partial u_{3}}{\partial z} + \frac{1}{2} [{(\frac{\partial u_{1}}{\partial z})}^{2} + {(\frac{\partial u_{2}}{\partial z})}^{2}] \\ ε_{x z} = \frac{1}{2} (\frac{\partial u_{1}}{\partial z} + \frac{\partial u_{3}}{\partial x}), ε_{y z} = \frac{1}{2} (\frac{\partial u_{2}}{\partial z} + \frac{\partial u_{3}}{\partial y}) \\ ε_{x y} = \frac{1}{2} (\frac{\partial u_{1}}{\partial y} + \frac{\partial u_{2}}{\partial x}) \end{array}

（1）

式中 $ε_{i j} (i, j = x, y, z)$ 为6个应变分量； $u_{1}, u_{2}, u_{3}$ 为与坐标系（ $x, y, z$ ）对应的位移场函数，其表达式为：

\{\begin{array}{l} u_{1} (z, t) = υ_{x} (z, t), u_{2} (z, t) = υ_{y} (z, t) \\ u_{3} (z, x, y, t) = υ_{z} (z, t) - x \frac{\partial υ_{x} (z, t)}{\partial z} - y \frac{\partial υ_{y} (z, t)}{\partial z} \end{array}

（2）

式中 $υ_{x}, υ_{y}, υ_{z}$ 分别表示油管在x，y，z方向上的位移（m）；x，y为横坐标（m）； $t$ 为时间（s）。

管内气体速度 $V$ 的水平分量和垂直分量分别为 $v_{x} = {\dot{υ}}_{x} + V υ_{x}^{'}, v_{y} = {\dot{υ}}_{y} + V υ_{y}^{'}, v_{z} = V$ 。得油管柱总动能 $T$ 、势能 $U$ 和外力做功 $W$ 的表达式为：

\begin{array}{l} T = \int_{0}^{L} \int_{A} \frac{ρ_{υ}}{2} ({\dot{u}}_{1}^{2} + {\dot{u}}_{2}^{2} + {\dot{u}}_{3}^{2}) d A + \frac{1}{2} m_{0} (v_{x}^{2} + v_{y}^{2} + v_{z}^{2}) d z = \\ \frac{1}{2} \int_{0}^{L} \{m_{υ} ({\dot{υ}}_{x}^{2} + {\dot{υ}}_{y}^{2} + {\dot{υ}}_{z}^{2}) + ρ_{υ} I ({\dot{υ}}_{x}^{'}^{2} + {\dot{υ}}_{y}^{'}^{2}) + \\ m_{0} [({\dot{υ}}_{x} + V υ_{x}^{'})^{2} + ({\dot{υ}}_{y} + V υ_{y}^{'})^{2} + V^{2}]\} d z \end{array}

（3）

\begin{array}{l} U = \frac{E}{2} \int_{0}^{L} \int_{A} [{υ_{z}^{'}}^{2} + x^{2} {υ_{x}^{″}}^{2} + y^{2} {υ_{y}^{″}}^{2} + \frac{1}{4} ({υ_{x}^{'}}^{2} + {υ_{y}^{'}}^{2})^{2} - \\ 2 x υ_{z}^{'} υ_{x}^{″} - 2 y υ_{z}^{'} υ_{y}^{″} + 2 x y υ_{x}^{″} υ_{y}^{″} + υ_{z}^{'} ({υ_{x}^{'}}^{2} + {υ_{y}^{'}}^{2}) - \\ x υ_{x}^{″} ({υ_{x}^{'}}^{2} + {υ_{y}^{'}}^{2}) - y υ_{y}^{″} ({υ_{x}^{'}}^{2} + {υ_{y}^{'}}^{2})] d A d z \end{array}

（4）

W = \int_{0}^{L} [p (z, t) d υ_{x} + q (z, t) d υ_{y} + f (z, t) d υ_{z}] d z

（5）

式中 $ρ_{υ}$ 为油管的密度（kg/m³）； $m_{υ}, m_{0}$ 分别为单元长度油管和气体的质量（kg）； ${\dot{υ}}_{x}, {\dot{υ}}_{y}, {\dot{υ}}_{z}$ 分别为油管三个方向位移随时间的一阶导数； $υ_{x}^{'}, υ_{y}^{'}$ 和 $υ_{x}^{″}, υ_{y}^{″}$ 分别为油管横向位移随z的一阶导数和二阶导数； $A$ 为油管横截面积（m²）； $V$ 为管内气体速度（m/s）； $I$ 为油管的惯性矩（m⁴）； $E$ 为油管弹性模量（Pa）； $f (z, t)$ ， $p (z, t)$ 和 $q (z, t)$ 分别为油管的x，y，z方向上的外力（N）。

由于油管是标准的圆柱体，积分满足以下公式

\{\begin{array}{l} \int_{A} x d A = \int_{A} y d A = \int_{A} x y d A = 0 \\ \int_{A} x^{2} d A = \int_{A} y^{2} d A = I \end{array}

（6）

根据哈密顿（Hamilton）变分原理 $δ \int_{t_{1}}^{t_{2}} (T - U + W) = 0$ ，结合公式（3）~（6）化简得到单元油管的x，y，z方向的振动控制方程：

\{\begin{array}{l} (m_{υ} + m_{0}) {\ddot{υ}}_{x} - ρ_{υ} I {\ddot{υ}}_{x}^{″} + 2 m_{0} V {\dot{υ}}_{x}^{'} + m_{0} V^{2} υ_{x}^{″} + \\ E I {υ_{x}}^{″}^{″} - E A {\{[υ_{z}^{'} + \frac{1}{2} ({υ_{x}}^{'}^{2} + {υ_{y}}^{'}^{2})] υ_{x}^{'}\}}^{'} = p (z, t) \\ (m_{υ} + m_{0}) {\ddot{υ}}_{y} - ρ_{υ} I {\ddot{υ}}_{y}^{″} + 2 m_{0} V {\dot{υ}}_{y}^{'} + m_{0} V^{2} υ_{y}^{″} + \\ E I {υ_{y}}^{″}^{″} - E A {\{[υ_{z}^{'} + \frac{1}{2} ({υ_{x}^{'}}^{2} + {υ_{y}^{'}}^{2})] υ_{y}^{'}\}}^{'} = q (z, t) \\ m_{υ} {\ddot{υ}}_{z} - E A [υ_{z}^{″} + \frac{1}{2} ({υ_{x}^{'}}^{2} + {υ_{y}^{'}}^{2})^{'}] = f (z, t) \end{array}

（7）

实际油管上端为油管挂，下端为封隔器，把上下端视为固定端，初始时刻油管不发生振动，可得模型的边界条件和初始条件为：

\{\begin{array}{l} υ_{x} (0, t) = 0, υ_{y} (0, t) = 0, υ_{z} (0, t) = 0, \\ υ_{x} (L, t) = 0, υ_{y} (L, t) = 0, υ_{z} (L, t) = 0 \\ \frac{\partial υ_{x} (0, t)}{\partial z} = 0, \frac{\partial υ_{y} (0, t)}{\partial z} = 0, \\ \frac{\partial υ_{x} (L, t)}{\partial z} = 0, \frac{\partial υ_{y} (L, t)}{\partial z} = 0 \\ υ_{x} (z, 0) = 0, υ_{y} (z, 0) = 0, υ_{z} (z, 0) = 0, \\ \frac{\partial υ_{x} (z, 0)}{\partial z} = 0, \frac{\partial υ_{y} (z, 0)}{\partial z} = 0, \\ \frac{\partial^{2} υ_{x} (z, 0)}{\partial^{2} z} = 0, \frac{\partial^{2} υ_{y} (z, 0)}{\partial^{2} z} = 0 \\ {\dot{υ}}_{x} (z, 0) = 0, {\dot{υ}}_{y} (z, 0) = 0, {\dot{υ}}_{z} (z, 0) = 0 \end{array}

（8）

1.2 载荷分析

三高气井油管在生产作业过程中，诱发其振动的外载荷主要包括管内高速气体的冲击载荷和油套管的接触碰撞载荷，而影响管内气体冲击载荷的因素包括井眼轨迹和井筒温度压力，为此，需建立三高气井油套管接触碰撞载荷计算方法、气体冲击载荷计算方法、井眼轨迹插值方法和井筒温度压力场计算方法，为油管模型的求解奠定载荷基础。

1.2.1 油套管非线性接触载荷

采用弹塑性力学理论^［

19］，建立油套管接触载荷计算方法。其变形结构如图2所示，油管发生碰撞后，油管上的

A_{2}

点变形到套管上的

A_{1}

点。由几何关系得：

R_{i}^{2} = {(R_{i}^{} - z_{i})}^{2} + r^{2} (i = 1,2)

（9）

式中 $R_{1}$ ， $R_{2}$ 分别为套管和油管的半径（m）； $z_{1}$ ， $z_{2}$ 分别为套管和油管接触点到套管内壁的径向距离（m）； $r$ 为接触点到油套管轴心的水平距离（m）。

图2 油套管接触变形示意图

Fig.2 Contact deformation between the tubing and casing

在力的作用下，套管和油管发生变形，它们之间的形变为 $δ$ ，并产生了宽为 $2 b$ 的接触带。套管和油管在接触后沿径向产生的位移分别为 $ω_{1}$ 和 $ω_{2}$ ，由几何关系可得：

\begin{array}{l} δ = (ω_{2} + ω_{1}) + (z_{2} - z_{1}) = \\ (ω_{2} + ω_{1}) + \frac{R_{1} - R_{2}}{2 R_{1} R_{2}} r^{2} \end{array}

（10）

式中 $δ$ 为油套管接触后产生的形变（m）； $z_{2} - z_{1}$ 为油套管接触边缘点移动的径向位移（m）； $ω_{1}$ ， $ω_{2}$ 分别为套管和油管接触后径向产生的位移（m）。

若接触带的宽度比油管的半径小很多时，则每个油管都可以近似地当作弹性半平面来考虑，由文献［

19］可得到

ω_{1}

和

ω_{2}

的计算公式：

\{\begin{array}{l} ω_{1} = - \frac{2 (1 - μ_{1}^{2})}{π E_{1}} \{\int_{- b}^{b} q (x) l n |r - x| d x - [\frac{1}{2 (1 - μ_{1})} - l n R_{1}] F\} \\ ω_{2} = - \frac{2 (1 - μ_{2}^{2})}{π E_{2}} \{\int_{- b}^{b} q (x) l n |r - x| d x - [\frac{1}{2 (1 - μ_{2}^{})} - l n R_{2}] F\} \end{array}

（11）

式中 $F$ 为油套管接触力（N）； $2 b$ 为油套管接触带宽度（m）； $E_{1}$ ， $E_{2}$ 分别为套管和油管材料弹性模量（Pa）； $μ_{1}$ ， $μ_{2}$ 分别为套管和油管材料泊松比； $q (x)$ 为套管和油管均布载荷（N）。

假定 $q (x)$ 是与以直径 $2 b$ 所作的半圆弧的纵坐标成比例，得^［

19］：

q (x) = \frac{2 F}{π b^{2}} \sqrt[]{b^{2} - x^{2}}

（12）

将 $E_{1} = E_{2} = E, μ_{1} = μ_{2} = 0.3$ 以及公式（11），（12）代入式（10），化简得油套管接触载荷和摩擦力计算方法；

\{\begin{array}{l} δ = 1.82 \frac{F}{E} [1 - l n (1.522 \sqrt[]{\frac{F}{π} \frac{R_{1} R_{2}}{R_{2} - R_{1}}})] \\ f = μ F \end{array}

（13）

式中 $f$ 为油套管摩擦力（N）； $μ$ 为油套管的摩擦系数，可由摩擦试验或材料手册确定^［

20］。

详细推导见前期研究工作^［

16］

1.2.2 管内高速流体冲击载荷

当管内高速气体经过井斜角变化区域时，对油管柱产生冲击载荷（如图3所示），引起油管振动。

图3 管内高速气体冲击载荷示意图

Fig.3 Schematic of the impact load by high-speed gas

根据流体力学^［

21］可知管内气体冲击载荷计算方法：

\{\begin{array}{l} F_{x} = - ρ_{0} A_{0} V^{2} s i n (α_{2} - α_{1}) c o s (φ_{2} - φ_{1}) \\ F_{y} = - ρ_{0} A_{0} V^{2} s i n (α_{2} - α_{1}) s i n (φ_{2} - φ_{1}) \\ F_{z} = - ρ_{0} A_{0} V^{2} c o s (α_{2} - α_{1}) \end{array}

（14）

式中 $ρ_{0}$ 为气体密度（kg/m³）； $A_{0}$ 为井筒横截面积（m²）； $α_{1}$ ， $α_{2}$ 分别为上下微段的井斜角（rad）； $φ_{1}$ ， $φ_{2}$ 分别为上下微段的方位角（rad）； $F_{x}$ ， $F_{y}$ ， $F_{z}$ 为 $x$ ， $y$ ， $z$ 方向的气体冲击力（N）。

由于气体冲击载荷与井眼轨迹有关，为此，采用三次样条插值方法^［

22］可求出任意井深处的井斜角和方位角，区间

[s_{k - 1}, s_{k}]

（k=1，2，…，N）上井斜角函数

α (s)

具体表达式为：

\{\begin{array}{l} α (s) = \frac{M_{k - 1} {(s_{k} - s)}^{3}}{6 L_{k}} + \frac{M_{k} {(s - s_{k - 1})}^{3}}{6 L_{k}} + \\ C_{k} (s - s_{k - 1}) + C_{k - 1} (s_{k} - s) \\ φ (s) = \frac{m_{k - 1} {(s_{k} - s)}^{3}}{6 L_{k}} + \frac{m_{k} {(s - s_{k - 1})}^{3}}{6 L_{k}} + \\ c_{k} (s - s_{k - 1}) + c_{k - 1} (s_{k} - s) \end{array}

（15）

式中 $C_{i} = \frac{α_{i}}{L_{k}} - \frac{M_{i} L_{k}}{6}, i = k, k - 1$ ； $c_{i} = \frac{φ_{i}}{L_{k}} - \frac{m_{i} L_{k}}{6}, i = k, k - 1$ ； $M_{i} = α_{i}^{″}, i = k, k - 1$ ； $m_{i} = φ_{i}^{″}, i = k,$ $k - 1$ ； $k$ 为测点序号； $L_{k}$ 为测段长度， $L_{k} = s_{k} - s_{k - 1}$ （m）； $s$ 为插值点处的井深（m）； $N$ 为测点个数。

1.2.3 井筒温度压力场

三高气井井筒温度压力随井深发生变化，引起管内气体流速、密度发生变化，从而导致油管柱不同位置处的冲击载荷发生变化，需建立井筒温度压力场计算方法，确定管内高速气体在不同井深处的流速和密度。根据前期的研究工作^［

23］得到井筒温度、压力、流速和密度耦合计算方法：

\{\begin{array}{l} \frac{d V}{d z} = - \frac{V d ρ_{0}}{ρ_{0} d z} \\ \frac{d p}{d z} = ρ_{0} g c o s α + {(\frac{d p}{d z})}_{f r} - \frac{ρ_{0} V d v}{d z} \\ \frac{d T}{d z} = \frac{g c o s α - \frac{V d V}{d z} + \frac{2 π r_{c o} k_{e} U_{a} (T - T_{e i})}{w_{i} (k_{e} + r_{c o} U_{a} f_{t})}}{c_{p}} \\ \frac{d ρ_{0}}{d z} = \frac{T ρ_{0} g c o s α + T {(\frac{d p}{d z})}_{f r} - \frac{ρ_{0} g c o s α}{c_{p}} - \frac{p}{c_{p}} \cdot \frac{2 π r_{c o} k_{e} U_{a} (T - T_{e i})}{w_{i} (k_{e} + r_{c o} U_{a} f_{t})}}{\frac{Z_{g} R T^{2}}{M} - T V^{2} + \frac{p V^{2}}{c_{p} ρ_{0}}} \end{array}

（16）

式中 $T$ ， $p$ ， $V$ ， $ρ_{0}$ 分别为井筒温度（^oC）、压力（Pa）和气体流速（m/s）、密度（kg/m³）； $k_{e}$ ， $U_{a}$ 分别为地层导热系数和井筒总导热系数（ $W / m \cdot ℃$ ）； $T_{e i}$ 为地层温度（^oC）； $r_{c o}$ 为套管外径（m）； $f_{t}$ 为瞬态传热时间函数； $Z_{g}$ 为天然气压缩因子，采用Gopal方法^［

24］；

R

为气体常量（J/（mol·k））；

M

为气体摩尔质量（kg/mol）；

c_{p}

为定压比热（

J / (k g \cdot K)

）；

w_{i}

为流体质量流量（kg/s）；

α

为井斜角（rad）；

{(\frac{d p}{d z})}_{f r}

为气体摩阻压力梯度。基于现场提供的井底或井口温度压力数据，采用四阶龙格‑库塔法数值求解出井筒不同位置处的温度、压力、流速和气体密度。

1.3 模型求解

1.3.1 有限元离散

采用拉格朗日函数和三次埃尔米特函数表示油管x，y，z方向位移，其有限元离散形式为：

υ_{x} = φ_{x}^{T} d, υ_{y} = φ_{y}^{T} d, υ_{z} = φ_{z}^{T} d

（17）

{[\begin{matrix} d \\ φ_{x} \\ φ_{y} \\ φ_{z} \end{matrix}]}^{T} = ([\begin{array}{l} {υ_{z}}_{1} \\ υ_{x 1} \\ \frac{d υ_{x 1}}{d z} \\ υ_{y 1} \\ \frac{d υ_{y 1}}{d z} \\ {υ_{z}}_{2} \\ υ_{x 2} \\ \frac{d υ_{x 2}}{d z} \\ υ_{y 2} \\ \frac{d υ_{y 2}}{d z} \end{array}], [\begin{array}{l} 0 \\ 1 - \frac{3 z^{2}}{l^{2}} + \frac{2 z^{3}}{l^{3}} \\ z - \frac{2 z^{2}}{l} + \frac{z^{3}}{l^{2}} \\ 0 \\ 0 \\ 0 \\ \frac{3 z^{2}}{l^{2}} - \frac{2 z^{3}}{l^{3}} \\ - \frac{z^{2}}{l} + \frac{z^{3}}{l^{2}} \\ 0 \\ 0 \end{array}], [\begin{array}{l} 0 \\ 0 \\ 0 \\ 1 - \frac{3 z^{2}}{l^{2}} + \frac{2 z^{3}}{l^{3}} \\ z - \frac{2 z^{2}}{l} + \frac{z^{3}}{l^{2}} \\ 0 \\ 0 \\ 0 \\ \frac{3 z^{2}}{l^{2}} - \frac{2 z^{3}}{l^{3}} \\ - \frac{z^{2}}{l} + \frac{z^{3}}{l^{2}} \end{array}], [\begin{array}{l} 1 - \frac{z}{l} \\ 0 \\ 0 \\ 0 \\ 0 \\ \frac{z}{l} \\ 0 \\ 0 \\ 0 \\ 0 \end{array}])

（18）

把位移函数式（18）代入单元油管振动控制方程（7）中，得系统的离散形式动力学方程：

M \ddot{D} + C \dot{D} + K D = F

（19）

式中 $D$ ， $M$ ， $C$ ， $K$ 和 $F$ 分别为系统的位移矩阵、质量矩阵、阻尼矩阵、刚度矩阵和载荷列向量。其表达式为：

\{\begin{array}{l} M = M_{1} + M_{2} + M_{3} + M_{4} + M_{5} \\ C = C_{1} + C_{2} + C_{3} + C_{4} + C_{5} \\ K = K_{1} + K_{2} + K_{3} + K_{4} + K_{5} + K_{6} + K_{7} + \\ K_{8} + K_{9} + K_{10} + K_{11} + K_{12} + K_{13} \end{array}

（20）

\{\begin{array}{l} M_{1} = (m_{υ} + m_{0}) \int_{0}^{l} φ_{x} φ_{x}^{T} d z \\ M_{2} = (m_{υ} + m_{0}) \int_{0}^{l} φ_{y} φ_{y}^{T} d z \\ M_{3} = ρ_{υ} I \int_{0}^{l} φ_{x}^{'} {φ_{x}^{'}}^{T} d z \\ M_{4} = ρ_{υ} I \int_{0}^{l} φ_{y}^{'} {φ_{y}^{'}}^{T} d z \\ M_{5} = m_{υ} \int_{0}^{l} φ_{z} φ_{z}^{T} d z \end{array}

（21）

\{\begin{array}{l} C_{1} = c \int_{0}^{l} φ_{x} φ_{x}^{T} d z \\ C_{2} = 2 m_{0} V \int_{0}^{l} φ_{x} {φ_{x}^{'}}^{T} d z \\ C_{3} = c \int_{0}^{l} φ_{y} φ_{y}^{T} d z \\ C_{4} = 2 m_{0} V \int_{0}^{l} φ_{y} {φ_{y}^{'}}^{T} d z \\ C_{5} = c \int_{0}^{l} φ_{z} φ_{z}^{T} d z \end{array}

（22）

\{\begin{array}{l} K_{1} = E I \int_{0}^{l} φ_{x}^{″} {φ_{x}^{″}}^{T} d z \\ K_{2} = E I \int_{0}^{l} φ_{y}^{″} {φ_{y}^{″}}^{T} d z \\ K_{3} = - m_{0} V^{2} \int_{0}^{l} φ_{x}^{'} {φ_{x}^{'}}^{T} d z \\ K_{4} = - m_{0} V^{2} \int_{0}^{l} φ_{y}^{'} {φ_{y}^{'}}^{T} d z \\ K_{5} = E A \int_{0}^{l} φ_{x}^{'} d^{T} φ_{z}^{'} φ_{x}^{T} d z \\ K_{6} = E A \int_{0}^{l} φ_{y}^{'} d^{T} φ_{z}^{'} {φ_{y}^{'}}^{T} d z \\ K_{7} = \frac{1}{2} E A \int_{0}^{l} φ_{x}^{'} {φ_{x}^{'}}^{T} d d^{T} φ_{x} φ_{x}^{T} d z \\ K_{8} = \frac{1}{2} E A \int_{0}^{l} φ_{y}^{'} {φ_{y}^{'}}^{T} d d^{T} φ_{y}^{'} {φ_{y}^{'}}^{T} d z \\ K_{9} = \frac{1}{2} E A \int_{0}^{l} φ_{x}^{'} {φ_{x}^{'}}^{T} d d^{T} φ_{y}^{'} {φ_{y}^{'}}^{T} d z \\ K_{10} = \frac{1}{2} E A \int_{0}^{l} φ_{y}^{'} {φ_{y}^{'}}^{T} d d^{T} φ_{x}^{'} {φ_{x}^{'}}^{T} d z \\ K_{11} = E A \int_{0}^{l} φ_{z}^{'} {φ_{z}^{'}}^{T} d z \\ K_{12} = \frac{1}{2} E A \int_{0}^{l} φ_{z}^{'} {φ_{x}^{'}}^{T} d {φ_{x}^{'}}^{T} d z \\ K_{13} = \frac{1}{2} E A \int_{0}^{l} φ_{z}^{'} {φ_{y}^{'}}^{T} d {φ_{y}^{'}}^{T} d z \end{array}

（23）

1.3.2 迭代求解

由于模型考虑的非线性因素过多，只用Newmark⁃β法进行逐步积分，会导致求解精度下降，不易收敛。为此，采用增量形式的Newmark⁃β法和Newton⁃Raphson法联合求解。首先给出 $t_{i}$ 和 $t_{i + 1}$ 时刻的控制方程：

\{\begin{array}{l} M {\ddot{u}}_{i} + C {\dot{u}}_{i} + K u_{i} = p_{i} \\ M {\ddot{u}}_{i + 1} + C {\dot{u}}_{i + 1} + K u_{i + 1} = p_{i + 1} \end{array}

（24）

由 $t_{i + 1}$ 减去 $t_{i}$ 时刻的运动方程可得其运动的增量形式的平衡方程：

M Δ {\ddot{u}}_{i} + C Δ {\dot{u}}_{i} + {(Δ f_{s})}_{i} = Δ p_{i}

（25）

式中

\{\begin{array}{l} Δ u_{i} = u_{i + 1} - u_{i} \\ Δ {\dot{u}}_{i} = {\dot{u}}_{i + 1} - {\dot{u}}_{i} \\ Δ {\ddot{u}}_{i} = {\ddot{u}}_{i + 1} - {\ddot{u}}_{i} \\ {(Δ f_{s})}_{i} = {(f_{s})}_{i + 1} - {(f_{s})}_{i} \\ Δ p_{i} = p_{i + 1} - p_{i} \end{array}

（26）

当时间步长 $Δ t$ 取得足够小时，认为在 $t_{i} - t_{i + 1}$ 区间内系统为线性的，如图4所示，即：

(Δ f_{s})_{i} = k_{i}^{s} Δ u_{i}

（27）

式中 $k_{i}^{s}$ 为 $t_{i}$ ， $t_{i + 1}$ 点之间的割线刚度。

图4 Newton-Raphson法迭代原理示意图

Fig. 4 Schematic of Newton-Raphson method iteration principle

则增量形式的平衡方程为：

m Δ {\ddot{u}}_{i} + c Δ {\dot{u}}_{i} + k_{i}^{s} Δ u_{i} = Δ p_{i}

（28）

采用增量形式的Newmark‑β法逐步积分求解：

{\overset{⌢}{k}}_{i} Δ u_{i} = Δ {\overset{⌢}{p}}_{i}

（29）

式中

\begin{array}{l} {\hat{k}}_{i} = k_{i} + \frac{1}{β Δ t^{2}} m + \frac{γ}{β Δ t} c, \\ {\hat{p}}_{i} = p_{i} + [\frac{1}{β Δ t^{2}} u_{i} + \frac{1}{β Δ t} {\dot{u}}_{i} + (\frac{1}{2 β} - 1) {\ddot{u}}_{i}] m + \\ [\frac{γ}{β Δ t} u_{i} + (\frac{γ}{β} - 1) {\dot{u}}_{i} + \frac{Δ t}{2} (\frac{γ}{β} - 2) {\ddot{u}}_{i}] c \end{array}

（30）

求得 $Δ u_{i}$ 后，计算 $t_{i + 1}$ 时刻的总位移：

u_{i + 1} = u_{i} + Δ u_{i}

（31）

利用 $t_{i + 1}$ 时刻的位移，计算 $t_{i + 1}$ 时刻的加速度和速度：

\{\begin{array}{l} {\ddot{u}}_{i + 1} = \frac{1}{β Δ t^{2}} Δ u_{i} - \frac{1}{β Δ t} {\dot{u}}_{i} - (\frac{1}{2 β} - 1) {\ddot{u}}_{i} \\ {\dot{u}}_{i + 1} = \frac{γ}{β Δ t} Δ u_{i} + (1 - \frac{γ}{β}) \frac{1}{β Δ t} {\dot{u}}_{i} + (1 - \frac{γ}{2 β}) {\ddot{u}}_{i} Δ t \end{array}

（32）

为提高Newmark⁃β法的求解精度，采用Newton⁃Raphson法对其位移进行进一步矫正，提高求解精度，图4表示Newton⁃Raphson法迭代原理图，在Newmark⁃β法逐步迭代的过程中，对等式（28）用其可求得 $Δ {u_{i}}^{(1)}$ ， $Δ {u_{i}}^{(2)}$ ，···， $Δ {u_{i}}^{(l)}$ 以后再进行了 $l$ 次迭代计算，位移叠加得：

Δ u = \sum_{j = 1}^{l} Δ {u_{i}}^{(j)}

（33）

若 $\frac{Δ {u_{i}}^{(l)}}{Δ u} < ε$ ，即满足收敛条件，退出Newton⁃ Raphson迭代， $ε$ 为收敛精度，取值为10^-6。采用Fortran语言编写了三高油管柱三维非线性流致振动响应计算代码，具体流程如图5所示。

图5 油管柱三维非线性流致振动模型求解流程图

Fig.5 Flow chart of solving 3D nonlinear flow induced vibration model of tubing string

1.4 模型验证

由于无法从现场获得有效的振动数据验证所建模型的正确性，故借助模拟试验方法验证模型的正确性。根据中国南海某气田的实际情况（产量200万方/天，地层压力54.93 MPa，地层温度151.14 ℃，常温25 ℃，标准大气压0.1 MPa），通过气体状态方程，确定管内模拟气体速度，得到模拟试验参数（如表1所示），搭建了油管柱模拟试验台架，如图6所示。

表1 模拟试验参数

Tab. 1 Simulation test parameters

参数名称

油管

外径/mm

油管

内径/mm

油管

长度/m

套管

外径/mm

套管

内径/mm

气体流速/

(m·s^-1)

油管弹性模量/GPa

油管密度/(kg·m^-3)

气体

密度/(kg·m^-3)

实际值

114.3

100.3

3500

177.8

152.5

7.58

210.0

7850

275

试验值

23.0

20.0

7.3

35.0

30.0

7.58

6.0

1200

275

图6 试验系统结构图

Fig.6 Structure of experimental system

采用建立的非线性振动模型，设置与模拟试验参数一样（如表1所示），采用模拟试验台架结构，其中油管分为300个单元，模拟总时间为70 s，步长为0.001 s，提取了4个位置处（与试验中安装传感器的位置相同）的振动响应。

由图7和8可知，通过非线性振动模型计算和试验测量所得到油管横向振动幅值和振动频率基本一致，其中两种方法都表明位移响应存在50 s的瞬态响应，试验结果的高频分量相对较多，主要由于实验环境因素的干扰。由幅频响应曲线可知油管在0~2.0 Hz范围内的振动能量较大，主频率集中在1.5 Hz左右。如表2所示，理论模型计算结果和实验测量结果精度高于92%，相对于二维非线性振动模型（对比精度高于88%）^［

25］的对比精度提高了4%。验证了所建立的三维非线性流致振动模型的正确性及有效性。

图7 油管不同测点的横向位移时程曲线

Fig.7 Time history curve of lateral displacement of tubing string corresponding to different positions

图8 油管不同测点的横向振动幅频响应曲线

Fig.8 Amplitude-frequency curve of the transverse vibration of tubing string corresponding to different positions

表2 理论计算结果与模拟实验结果对比表

Tab. 2 Comparison of theoretical calculation results and simulation experiment results

类型	振动幅值/mm				主频率 /Hz
类型	测点 2	测点3	测点6	测点7	测点2	测点3	测点6	测点7
理论计算	0.03408	0.07176	0.04039	0.09287	1.475	1.500	1.550	1.550
实验测量	0.03696	0.07412	0.04127	0.0946	1.576	1.615	1.615	1.615
精度	92.2%	96.8%	97.9%	98.2%	93.6%	92.9%	96.0%	96.0%

由图8可知，实验测量结果与理论计算结果在低频段和高频段存在误差，其中低频段存在误差（靠近台架固定位置，测点2和测点7）的主要原因是外界振动的干扰，由于实验场地存在其他实验仪器，在测量时，其他仪器运行产生的振动引起实验结果在低频段的幅值较大，理论模型无法考虑外界因素，其幅值相对较小。而在试验测点3位置处恰好存在一个接头，其外径略大于管柱的外径，阻碍了管柱的振动，导致此位置处的幅值反而更小。在管柱高频段振动区域，实验测量结果小于理论模型计算结果，其主要原因是实验采用了小波分析法处理了实验数据，将能量小的高频振动进行了过滤处理，使得其幅值较小。以上外界原因对主频率的振动幅值影响较小，主要原因是管柱的主频振动包含了整个管柱振动的能量，反映了管柱主要振动响应特性及振动趋势，为此，对比主频率的幅值能够有效验证模型的正确性。

2 案例分析

根据中国南海M三高气井的井眼轨迹（如图9（a）所示）、井身结构（如图9（b）所示）及参数（如表3所示），计算得到井筒温度压力分布曲线（如图9（c）和（d）所示），预测了三高气井油管柱振动情况，揭示了三高气井油管柱非线性振动特性，具体分析如下。

图9 中国南海M三高气井井眼轨迹、井身结构及井筒温压曲线

Fig.9 Well trajectory， wellbore structure and wellbore temperature pressure curve of M 3H gas well in South China Sea

表3 中国南海M三高气井计算参数

Tab. 3 Calculation parameters of M 3H gas well in South China Sea

管长/m

油管

内径/m

油管外径/m

套管

内径/m

套管

外径/m

产量/（10⁴m³·d^-1)

生产封隔器位置/m

计算

时间/s

时间

步长/s

油管密度/

(kg·m^-3)

3900

0.1003

0.1143

0.1658

0.1778

160

3627

0.001

7850

摩擦系数

屈服强度/

MPa

井斜角/（^o）

地层导热系数

井筒总导热系数

气体常量/

[J·（mol·K）^-1]

天然气摩尔质量/(g·mol^-1)

井底

温度/℃

井底

压力/MPa

流体密度/

(kg·m^-3)

0.243

665

0-43.44

2.06

97.14

8.314

151.14

54.93

275

根据管柱的空间振动轨迹（如图10（a）所示）和平面振动轨迹（如图10（b）所示），再结合图9（b）管柱的井深结构，可知直井段管柱的横向振动最明显，其次造斜段，再次下部稳斜段，最后是中部管柱，且中部管柱长时间紧靠套管壁做x⁃y平面的振动，易导致油管柱在此位置发生摩擦磨损失效；同时表明，当现场油气井的方位角发生显著变化时，仅考虑管柱纵横两个方向的振动，对管柱的磨损分析将不太精确，需考虑管柱的三维振动引起的纵向和横向的振动磨损。同时表明中下部管柱的磨损比上部磨损将更加严重，现场需重点关注中下部管柱的磨损安全性。

图10 油管柱不同时刻振动轨迹图

Fig.10 Vibration trajectory of tubing string at different time

图11为不同位置处管柱三个方向的振动响应曲线，图11（a）和（b）同样表明，上部管柱的x方向振动明显，且在套管间来回振动，振动频率显著高于其他位置；存在两主振动频率，中下部管柱的振动情况复杂，但振动频率很小，主要集中于0~0.3 Hz之间，由于双封隔器的作用，底部油管紧靠于套管壁上。

图11 不同位置油管柱振动响应图

Fig.11 Vibration response of tubing string at different positions

由图11（c）和（d）可知，上部油管y轴方向的振动同样最明显，振动频率比x轴方向的振动更加复杂，其主要原因是M气井的方位角变化比井斜角变化明显，使得气体在y轴方向的外激励也变得更加复杂。同时发现油管横向振动的频率都集中于0~1.2 Hz，属于复杂的低频振动。

图11（e）和（f）为三高气井油管柱纵向振动位移和幅频响应曲线，由于上端为油管挂，下端为封隔器，都为固定端，使得管柱的纵向位移呈现上下两端形变小，中间形变大的情形，其中形变量包含了油管自重的静态形变，但振动幅值在中下部最大，导致管柱的轴力及交变应力在中下部最大，将引起下部油管的磨损和疲劳损伤最大。再次表明，中下部油管的安全性最小，现场设计人员需重点关注。

图12为不同位置油管纵向方向相轨迹图，由图可知，不同位置油管都出现混沌振动现象（其相轨迹图呈现多层的封闭环，且出现局部密集现象，即局部吸引子），其中上端和下端的混沌振动更加明显，其主要原因是所建立的模型是纵横向耦合模型，即纵向振动和横向振动相互影响，相互牵制，由前面分析可知，上下段油管的横向振动更加剧烈，中部及中下部管柱横向振动很小，导致上下两端的混沌现象更加明显，由此表明，下一步研究将重点采用混沌理论分析管柱的非线性行为，找出外界参数对油管Hopf分岔性质的影响，有效控制油管柱的混沌振动，为现场油管柱优化设计奠定理论基础。

图12 不同位置油管柱纵向方向相轨迹图

Fig. 12 Phase trajectories in longitudinal direction of tubing string at different positions

3 结论与建议

针对三高气井流体诱发油管柱振动失效问题，采用理论与实验相结合的方法，开展了油管柱三维非线性流致振动模型研究，探究了油管柱振动响应特性，形成了以下结论及建议：

（1）采用能量法、微元法结合哈密顿原理建立了三高气井油管柱三维非线性流致振动模型，模型考虑了井眼轨迹变化、井筒温度压力变化和油套管非线性接触碰撞等因素，由于模型考虑的非线性因素复杂，其求解异常复杂，故采用拉格朗日函数和三次埃尔米特函数离散振动控制方程，并借助增量形式的Newmark⁃β法和Newton⁃Raphson法联合迭代求解离散后的控制方程，实现了三高气井油管柱三维非线性振动模型的数值求解；搭建了模拟实验台架，测得油管柱横向振动位移和幅频曲线，验证了油管三维非线性流致振动模型的正确性。

（2）采用所建立的油管柱三维非线性振动模型，借助中国南海M三高气井实例井参数，分析了油管柱的振动特性，发现上部油管三维空间横向振动明显，而中下部油管长时间紧靠于套管壁上发生振动，导致下部油管易发生纵向和横向的磨损，需重点关注此位置的磨损。由于套管的约束作用，油管柱横向振动幅值远小于纵向振动幅值，但其横向振动频率比纵向振动频率复杂，故不能忽略油管横向振动。同时发现油管纵向振动出现混沌现象，且中下部位置油管混沌现象最明显，因此，现场设计人员需重点关注下部油管柱的安全性。

（3）目前研究集中建立了油管柱三维非线性振动模型，对油管柱振动失效机理还未开展，建议后期可重点开展三高气井油管柱疲劳、稳定性及摩擦磨损分析方法研究，探究三高气井油管柱振动失效机理，提出更加全面的安全控制方法。

参考文献

范青，练章华，邓玮，等.三高气井管柱损伤研究现状［J］.重庆科技学院学报（自然科学版）， 2010， 12（2）：63-67. [百度学术]

Fan Qing， Lian Zhanghua， Deng Wei， et al. Three high gas wells column injury status research and development［J］. Journal of Chongqing University of Science and Technology （Natural Sciences Edition）， 2010， 12（2）：63-67. [百度学术]

刘巨保，王秀文，岳欠杯，等.受压油管屈曲变形与内外层杆管接触分析［J］.力学与实践， 2011， 33（3）：25-29. [百度学术]

Liu Jubao， Wang Xiuwen， Yue Qianbei， et al. The buckling deformation and contact with inner rod and outer casing for compressional tubing［J］. Mechanics in Engineering， 2011， 33（3）：25-29. [百度学术]

Aitken J. An account of some experiments on rigidity produced by centrifugal force［J］. Philosophical Magazine， 1878（5）：81-105. [百度学术]

Shilling R， Lou Y K. An experimental study on the dynamic response of a vertical cantilever pipe conveying fluid［J］. Journal of Energy Resources Technology， 1980， 102（3）：129-135. [百度学术]

邓元洲. 高产气井油管柱振动机理分析及疲劳寿命预测［D］. 成都：西南石油大学， 2006. [百度学术]

Deng Yuanzhou. Vibration mechanism analysis and fatigue life prediction of tubing string in high production gas well［D］. Chengdu： Southwest Petroleum University， 2006. [百度学术]

巨全利，仝少凯.高产气井完井管柱纵向振动特性分析［J］.钻采工艺，2014， 37（2）： 79-82. [百度学术]

Ju Quanli， Tong Shaokai. Analysis of vertical vibration characteristics of completion string of high production gas well［J］. Drilling and Production Technology， 2014， 37 （2）： 79-82. [百度学术]

Bagdatli S. Non-linear transverse vibrations of tensioned nanobeams using nonlocal beam theory［J］. Structural Engineering and Mechanics， 2015， 55（2）：281-298. [百度学术]

于凯强.高产气井完井管柱流固耦合振动分析及安全性评价［D］. 西安：西安石油大学，2017. [百度学术]

Yu Kaiqiang. Fluid solid coupling vibration analysis and safety evaluation of completion string of high production gas well［D］. Xi'an： Xi'an University of Petroleum， 2017. [百度学术]

桑松，褚兆飞，曹爱霞，等. 深海TTR立管涡激-参激耦合振动特性［J］. 哈尔滨工程大学学报， 2019，40（11）： 1823-1832. [百度学术]

Sang Song， Chu Zhaofei， Cao Aixia， et al. Vortex-parametric coupled vibration characteristics of deep-sea top-tension riser［J］. Journal of Harbin University of Engineering， 2019， 40（11）： 1823-1832. [百度学术]

邢誉峰，梁昆.梁纵向与横向耦合非线性振动分析［J］.北京航空航天大学学报，2015，41（8）：1359-1366. [百度学术]

Xing Yufeng， Liang Kun. Nonlinear vibration analysis of longitudinal-transverse coupled beam［J］. Journal of Beijing University of Aeronautics and Astronautics， 2015， 41（8）：1359-1366. [百度学术]

Liu J， Zhao H， Liu Q， et al. Dynamic behavior of a deepwater hard suspension riser under emergency evacuation conditions［J］. Ocean Engineering， 2018， 150：138-151. [百度学术]

Liu J， Zhao H， Yang S， et al. Nonlinear dynamic characteristic analysis of a landing string in deepwater riserless drilling［J］. Shock and Vibration， 2018： 8191526. [百度学术]

谈梅兰.三维曲井内钻柱的双重非线性静力有限元法［D］. 南京：南京航空航天大学， 2005. [百度学术]

Tan Meilan. Dual nonlinear static finite element method for drill string in three-dimensional curved well［D］. Nanjing： Nanjing University of Aeronautics and Astronautics， 2005. [百度学术]

王文明，熊明皓，陈钱荣，等. 深海垂直管中管载荷传递模拟分析［J］. 石油矿场机械， 2015， 44（3）：1-5. [百度学术]

Wang Wenming， Xiong Minghao， Chen Qianrong， et al. Simulation analysis of load transfer in deep sea pipe-in-pipe vertical system［J］. Oil Field Equipment， 2015， 44（3）：1-5. [百度学术]

Li C N， Samuel R. Buckling of concentric string pipe-in-pipe［C］. The SPE Annual Technical Conference and Exhibition， San Antonio， Texas， USA， October 2017. [百度学术]

Liu J， Guo X， Wang G， et al. Bi-nonlinear vibration model of tubing string in oil & gas well and its experimental verification［J］. Applied Mathematical Modelling， 2020， 81： 50-69. [百度学术]

Gao D L， Gao B K. A method for calculating tubing behavior in HPHT wells［J］. Journal of Petroleum Science & Engineering， 2004， 41（1）， 183-188. [百度学术]

Zhang Z， Wang H. Sealed annulus thermal expansion pressure mechanical calculation method and application among multiple packers in HPHT gas wells［J］. Journal of Natural Gas Science and Engineering， 2016， 31：692-702. [百度学术]

徐秉业，黄炎. 弹塑性力学及其应用［M］. 北京：机械工业出版社， 1984. [百度学术]

Xu Bingye， Huang Yan. Elastoplastic Mechanics and Its Application［M］. Beijing： Machinery Industry Press， 1984. [百度学术]

温诗铸，黄平，田煜，等. 摩擦学原理［M］. 5版. 北京：清华大学出版社， 2018. [百度学术]

Wen Shizhu， Huang Ping， Tian Yu， et al. Principles of Tribology ［M］. 5th ed. Beijing： Tsinghua University Press， 2018. [百度学术]

罗惕干. 流体力学［M］. 北京：机械工业出版社， 2017. [百度学术]

Luo Tigan. Fluid Mechanics［M］. Beijing： Machinery Industry Press， 2017. [百度学术]

刘清友，何玉发. 深井注入管柱力学行为及应用［M］. 北京：科学出版社， 2013. [百度学术]

Liu Qingyou， He Yufa. Mechanical Behavior of Injection String in Deep Well and Its Application［M］. Beijing： Science Press， 2013. [百度学术]

Liu J， He X， Huang C， et al. Calculation of shaft temperature field in deep-water testing system［J］. Mathematical and Computer Modelling， 2013， 17：102-111. [百度学术]

Gopal V. Gas z-factor equations developed for computer［J］. Oil and Gas Journal， 1977， 75（32）： 58-60. [百度学术]

Guo X， Liu J， Dai L， et al. Friction-wear failure mechanism of tubing strings used in high-pressure， high-temperature and high-yield gas wells［J］. Wear， 2021， 468-469： 203576. [百度学术]