钻具旋转时同心窄间隙环空流动的数值模拟

摘要

钻井过程中，当循环钻井液在井眼环空中作轴向流动时，由于钻具旋转而导致周向运动的出现，整个环空中因此形成螺旋流动，从而可能会对环空（特别是窄间隙）摩擦压降存在潜在的影响。对此，本文借助数值模拟软件，对同心窄间隙环空中钻具旋转下的层流流动进行直接数值模拟计算并分析。其中，钻具理想化为钻杆或圆柱体；钻井液则选用3种不同流变参数的宾汉流体。当前数值模拟计算结果表明，在不同钻具/杆转速下，PLR（压损比）几乎保持不变；这就表明，钻具/杆旋转过程中，并不会引起环空中摩擦压降的显著变化；然而，对比一些经验公式后发现，除了Ooms（1999）的公式能与当前计算吻合外，其余公式均不能反应该现象。进一步地，通过观察环空压力沿径向的变化过程时，结果表明，在当前180 rad/min范围内，压力沿径向几乎不变；即使当转速增大到600 rad/min，这种径向变化也十分有限，占比大约是2×10^-5量级。因此，在当前工程应用范围（钻具转速在180 rad/min以内）内，钻具旋转既不会导致轴向摩擦压降增大，也不会引起压力沿着径向显著变化。本文研究结果为实际钻井过程中对钻具旋转引起的摩擦压降的估算提供了一些参考，有助于优化钻井工艺和设备设计，提高钻进效率和安全性。

关键词

钻具旋转; 窄间隙; 环空流动; 宾汉流体; 摩擦压降; 数值模拟

0 引言

钻井过程中，钻井液在环空中循环时，除了主要的轴向流动之外，还会由于钻具旋转而引起周向运动；这两种运动耦合起来构成螺旋状流动。另一方面，窄间隙环空中，由于间隙过小，相较常规环空中的流动，相同体积排量下的钻井液流动加速，从而导致流动阻力大大增加。由此，窄间隙环空下的钻具旋转由于环空尺寸小、摩擦面积大以及钻井液粘度高等特点会对环空摩擦压降产生潜在的重要影响^［1-3］，因此开展窄间隙环空中循环钻井液时钻具旋转的影响研究很有必要。

对于钻具/杆类旋转引起的摩擦压降变化的研究，罗敏等^［4］考虑了接触非线性特性和动力学特性，研究了旋转钻柱流固耦合的效果，使钻柱的受力状态得到合理描述，为钻柱优化设计和井眼轨迹控制提供可行的方法。Hemphill等^［5］在压降和当量循环密度（Equivalent Circulating Density，ECD）的研究过程中发现层流状态下，随着钻具/杆旋转速度增大，环空压降几乎呈线性增大，相较于环空间隙逐渐变窄导致ECD增大，钻具/杆旋转导致ECD增大的影响更为明显。通过MODFLOW/MT3DMS数值模拟方法，研究现有模型的可能存在的误差，并提出改进的方法，从而提高SWIW试验的模拟精度。也有部分国外学者^［4-10］认为钻具/杆旋转速度和摩擦压降之间的关系是复杂的，常常受到流体属性、流动区域、钻杆旋转条件、直径比和偏心度的影响，对此进行了计算模型的修正。其中Ooms等^［11］、田野等^［12］、蔡萌^［13］则针对偏心条件下钻具/杆旋转进行研究，并给出了摩擦压降比例（Pressure Loss Ratio-PLR）公式。钻井工程通常认为旋转管柱会对摩阻计算造成影响，产生的原因可能为其他因素导致，例如管柱运动导致的偏心度变化、摆动频率变化等^［14-17］，导致流动状态发生变化，有待进一步研究。然而，除了数值模拟外，现有一些经验公式也经常用于此类问题的预测计算，这些经验公式通常基于实验测量结果或理论计算结果，使用起来较为简单方便，这些工程经验公式众多，计算参数复杂，并且存在特定的适用范围，特别是同心环空下钻具/杆旋转与摩擦压降的数值计算；对于非线性方程，通常需要近似公式进行初步计算^［18-19］，且不难发现这些公式的计算结果并不一致，这就给实际钻井过程中对钻具旋转引起的摩擦压降的估算带来了诸多不便，甚至使计算的精度不足。反观数值模拟，能够更加准确地预测钻具/杆旋转对宾汉流体摩擦压降的变化规律，可以提供更为可靠的预测结果。通过数值模拟的结果，我们可以了解旋转速度、轴间距和流体粘度等因素对摩擦压降变化的影响，并可以得到摩擦压降随这些参数变化的定量关系。

综合以上分析，并且考虑到作为非牛顿流体的钻井液，主要具有宾汉塑性流变特征，因此，本文主要选取了3种流变参数不同的宾汉塑性模型（Bingham模型、Herschel-Bulkley模型、Power-Law模型），这些公式基于流体的粘度和剪切应力来推导摩擦压降的计算公式，并对钻具/杆施加不同的旋转速度，针对同心窄间隙情况下，当钻具/杆旋转时对宾汉流体循环过程中的摩擦压降进行直接数值模拟，探究分析旋转速度与摩擦压降的关系。并考察当前经验公式及其应用在此类问题上预测的精确性，判断是否符合数值模拟结果定性或定量上的变化。

1 计算模型

为进一步分析比较钻具/杆旋转对摩擦压降的影响，本文从直接数值模拟角度，对宾汉流体的摩擦压降进行模拟计算。

1.1 控制方程及边界条件

在直角坐标系下研究同心环空流动。其中，假设流态为层流，流体不可压缩，并且忽略重力，而钻具/杆则抽象为环空内侧固体壁面。由此，该流体动力控制方程为质量守恒方程和Navier-Stokes方程，分别表达为：

$$\nabla \cdot u=\mathrm{0}$$

（1a)

$$\rho \frac{\partial u}{\partial t}+\rho \left(u\cdot \nabla \right)u=-\nabla \cdot p+\nabla \cdot \tau$$

(1b)

式中：$\nabla$——梯度算子；$u$——速度矢量；$p$——压力；t——时间；$\tau$——应力张量，其具体表达式和非牛顿流体的流变模型有关。

当前研究的是宾汉流变模型：

$$\tau ={\tau }_{\mathrm{y}}+{\mu }_{\mathrm{p}}\dot{\gamma }$$

式中：τ_y——屈服应力；μ_p——塑性粘度系数；$\dot{\gamma }$——剪切速率。

初始条件中，当t=0时，除了入口边界流速分布外，整个计算域内流速为零。

边界条件中，不考虑钻具/杆旋转时，计算域入口为均匀速度边界，即轴向流速等于入口剖面平均流速，而其它分量均为零；出口为自由出流边界条件，即流速分量沿法向梯度为零；固壁为无滑移边界条件，即流速恒为零。参考压力为入口中心压力为零。当考虑钻具/杆旋转时，环空内壁面周向速度具有钻杆线速度，而环空外侧壁面的周向速度则为零。

1.2 数值计算方法

本文选用宾汉流体模型在同心环空层流流动工况下进行摩擦压降计算，环空统一设计为典型的窄间隙环空，其中环空内壁、外壁直径分别为D_e=0.1143 m与D_w=0.1219 m，环空当量直径为D_h=0.0076 m。

对于网格无关性的验证以及环空轴向尺寸的选取，参考Ahmed等^［8］关于层流模型固井ECD模拟计算，选用了30D_h的轴向尺寸环空作为计算区域，网格尺寸调节参数环空剖面网格尺寸δ与轴向尺寸网格尺寸Δ_y分别为δ=0.001 mm、Δ_y=2 mm，如图1（a）所示。其中环空横截面由两侧边界层向横截面中心以R=1.1膨胀率生成网格，贴近避免处的网格采用0.001量级（如图1b所示），总网格数量为546000个，此情况下可较好保证数值模拟的精度。

图1  网格尺寸

Fig.1  Mesh size diagram

采用FLUENT软件进行直接数值模拟计算，涉及到控制方程的求解和收敛问题。为了准确地模拟实际情况，需要采用符合问题特点的控制方程。针对本问题，采用了非定常的PISO算法，较好地解决了从初始流动到最终稳定流动的过渡阶段问题。同时，整个控制方程组的收敛精度/误差控制在10^-4以内，以保证数值模拟的精度。

在控制方程中，压力-速度的耦合算法采用了非定常的PISO算法，其通过对速度场和压力更为准确地模拟，能够降低模拟误差并实现稳定的控制方程求解。此外，压力控制方程采用二阶离散格式，能够更为准确地捕捉到液流在腔体中的压力分布情况。动量方程则采用二阶迎风格式，能够更好地模拟多种复杂流动现象。

1.3 计算结果检验

按照以上参数设置，开展了计算域和网格无关性检验，相关工作已经完成，详细细节参见文献［10］。通过模拟发现：出口界面的最大轴向速度达到稳定值，且能较好满足实际要求精度，网格的数量也在可控范围以内，而且在各个钻具/杆的转速下都能实现速度稳定且较快达到收敛，可以较好满足仿真模拟的要求。

2 结果及分析

为进一步排除流体流变参数在模拟中的影响因素，选用了3种不同的流变参数的宾汉流体（如表1所示）进行比较。

表1  3种宾汉流体的参数

Table 1  Parameters of three Bingham fluids

流体编号	$\rho$/(g·cm-3)	Recr	ReMR	${\mu }_{\mathrm{p}}$/(Pa·s)	${\tau }_0$/Pa
1	1.74	3003	626	0.1981	15.89
2	2.10	3001	522	0.2882	9.85
3	1.94	3001	622	0.2237	3.72

环空钻井液循环时计算的体积排量均为0.8 m³/min，对应剖面平均流速为9.457 m/s；根据文献［10］从而得到循环雷诺数$R{e}_{\mathrm{M}\mathrm{R}}$。如表1所示，均小于临界雷诺数$R{e}_{\mathrm{c}\mathrm{r}}$，因此该轴向循环流动属于层流流态。为了考察不同钻具/杆转速对循环时的摩擦压降影响，根据文献中报道的测试范围^［4-7］，对3种宾汉流体都分别取转速（T）为0、10、20、30、60、120和180 rad/min（1 r/min=2π rad/min，下同）；其中，当最大转速为180 rad/min时，对应最大线速度为0.17 m/s，小于轴向流速，因此也是层流运动。

环空轴向上的循环摩擦压降计算方法，首先在环空入口下方15D_h处与出口上方10D_h处取2个横截面；由于该2个横剖面之间的流动完全稳定^［8］，因而可以计算该横截面的剖面压力平均值；从而，依据压力梯度定义，计算2个横剖面间的摩擦压降梯度值。

为更好分析比较前人经验公式在此工况下的优劣情况，将各个转速下获取的摩擦压降与零转速时的相比获得各个转速下（含零转速）的PLR（压损比），并将每组宾汉流体数据与以往的经验公式（本文主要选用Hemphill等^［6］、Ahmed等^［8］、Enfits等^［10］和Ooms等^［11］的研究结果作为参考）在同种工况下求得的PLR（压损比）进行比对，比对情况如图2所示。

图2  3种宾汉流体与4种计算模型压损比随转速变化的对比

Fig.2  Comparison of pressure‑loss ratio of three Bingham fluids with four calculation models with rotational speed

计算结果表明，随着钻具转速增大，直接数值模拟得到的摩擦压降梯度并没有发生任何显著变化，PLR（压损比）几乎保持常数为1，即非零转速下的摩擦压降与零转速下的摩擦压降相同。而对于各种经验公式，却表现出不完全相同的变化趋势。比如，随着旋转速度的增大，Ahmed算法的压损比在施加旋转速度的初始时便出现急剧上升的现象，并随着转速的增加逐渐呈线性增长。而与Ahmed算法相反，Enfits提供的计算模型在施加转速初始时压损比出现大幅度陡降，但随着转速的增加也逐渐呈线性上升。相反，Hemphill算法全程变化较为平缓，但由图可看出，随着旋转速度增大，环空压降全程几乎呈线性增大。而本文数据与Ooms算法得到的数据更为相似，Ooms等^［11］通过环柱空间流体运动方程解析求解得到轴向速度分布$v_{\mathrm{z}}^{\left(\mathrm{0}\right)}\left(\mathrm{无纲}\mathrm{量形}\mathrm{式}\right):$

$$\frac{v_{\mathrm{z}}^{\left(\mathrm{0}\right)}}{-\frac{d^{\mathrm{2}}}{\mu }\frac{\partial p^{\left(\mathrm{0}\right)}}{\partial z}}=\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{h}{d}\right)}^{\mathrm{2}}\left(\frac{y}{h}\right)\left(\mathrm{1}-\frac{y}{h}\right)$$

（2）

式中：$\mu$——流体动态粘度；$p^{\left(\mathrm{0}\right)}$——压力；$z$——轴向坐标；$y$——钻杆间距；$h$——环形间隙；$d$——井壁半径与钻杆半径差值。

对式（2）进行积分得到零阶环形截面的体积通量$Q^{(\mathrm{0})}$，横截面上的零阶轴向速度分布（δ=3.0、ε=2.0、Ta=6）：

$$\frac{Q^{(\mathrm{0})}}{-\frac{d^{\mathrm{3}}{R}_{\mathrm{0}}}{\mu }\frac{\partial p^{(\mathrm{0})}}{\partial z}}=\frac{\mathrm{\pi }}{\mathrm{12}}\left(\mathrm{2}+\mathrm{3}{\epsilon }^{\mathrm{2}}\right)$$

（3）

式中：ε——无纲量偏心率；$R_{\mathrm{0}}$——井壁半径。

轴流一阶校正：

$$Q^{(\mathrm{1})}=\mathrm{0}$$

（4）

轴向速度二阶校正后积分，二阶轴流校正$Q^{(\mathrm{2})}$（δ=0.2、ε=0.6、Ta=1）：

$$\begin{array}{l}\frac{Q^{(\mathrm{2})}}{-\frac{d^{\mathrm{3}}{R}_o}{\mu }\frac{\partial p^{(\mathrm{0})}}{\partial z}}=\frac{\mathrm{\pi }}{\mathrm{4435200}}{\delta }^{\mathrm{2}}{(Ta)}^{\mathrm{2}}\frac{{\epsilon }^{\mathrm{2}}}{(\mathrm{2}+{\epsilon }^{\mathrm{2}}{)}^{\mathrm{2}}}\times \\ (-\mathrm{13120}+\mathrm{38112}{\epsilon }^{\mathrm{2}}-\mathrm{10608}{\epsilon }^{\mathrm{4}}-\mathrm{33062}{\epsilon }^{\mathrm{6}}-\mathrm{1221}{\epsilon }^{\mathrm{8}})\end{array}$$

（5）

式中：δ——无纲量间隙宽度；Ta——泰勒数。

结合式（3）、（4）、（5）可以得到钻杆旋转的轴向压降和不旋转的轴向压降之比：

$$\begin{array}{cc}& \begin{array}{l}\frac{\left(-\frac{\partial p}{\partial z}\right)}{{\left(-\frac{\partial p}{\partial z}\right)}_{Ta=\mathrm{0}}}=\left[\left(\mathrm{1}+\frac{\mathrm{3}}{\mathrm{2}}{\epsilon }^{\mathrm{2}}\right)\right]\times \\ {\left[\left(\mathrm{1}+\frac{\mathrm{3}}{\mathrm{2}}{\epsilon }^{\mathrm{2}}\right)+\frac{\mathrm{1}}{\mathrm{739200}}{\delta }^{\mathrm{2}}{\left(Ta\right)}^{\mathrm{2}}\frac{{\epsilon }^{\mathrm{2}}}{{\left(\mathrm{2}+{\epsilon }^{\mathrm{2}}\right)}^{\mathrm{2}}}f\left(\epsilon \right)\right]}^{-\mathrm{1}}\end{array}\\ & \begin{array}{l}f(\epsilon )=-\mathrm{13120}+\mathrm{38112}{\epsilon }^{\mathrm{2}}-\mathrm{10608}{\epsilon }^{\mathrm{4}}-\\ \mathrm{33062}{\epsilon }^{\mathrm{6}}-\mathrm{1221}{\epsilon }^{\mathrm{8}}\end{array}\end{array}$$

（6）

当处于同心状态（ε=0）时，由式（6）可以看出压降并未增加，得到的压损比持续为1，并不会随着转速的增大而线性增大，似与本文模拟结果更加吻合。

针对摩擦压降梯度随着钻具/杆旋转速度的增大并未发生变化的现象，进一步分析同心环空中流动压力沿着径向分布特征和钻具/杆转速的定量关系。选取环空L=15D_h处，沿径向的一条直线分别作T为0、120、180和600 rad/min的压力分布曲线图。

以1号宾汉流体为例（如图3所示），在相关文献报道的计算范围（转速<120 rad/min）内^［5-9］，转速对于井壁压力几乎没有影响；而在当前（工程）测试范围（转速<180 rad/min），井壁压力有轻微变化但变化的值域仅在1 Pa的范围内，因此也可以忽略不计。为此，本文补充计算了更高的钻具转速，达到600 rad/min时。如图3（d）所示，当钻具转速达到600 rad/min时，井壁压力变化梯度较为明显，钻具旋转导致环空内侧壁面附近的流体拥有较大的旋转速度，受到的离心力增大，这是由于旋转钻具对钻井液体产生向外的离心力的作用。流体在向井眼深处流动时，会受到内侧壁面的限制而产生强烈的摩擦作用，这导致液体相对静止，流动速度变慢。而在环空外侧或井壁旁边，液体由于受到离心力作用，具有向外的运动趋势，形成一定的流动速度，这使得环空外侧或井壁外的压力增大。然而，即使如此高的转速下，这种压力增长也只是在10 Pa量级，占当地压力比值约为10^-5量级，因此也可以将其忽略。

图3  1号宾汉流体在不同转速下的压力沿径向分布

Fig.3  The radial distribution of pressure of Bingham fluid No.1, at different speeds

3 结论与建议

（1）数值计算可以发现，当钻具/杆旋转速度持续增大到180 rad/min时，同心窄间隙环空的压损比PLR基本保持不变。结合当量循环密度ECD计算可知，钻具/杆旋转过程中，并不会引起环空中摩擦压降的显著变化，钻具旋转导致井壁压力的增大，等同于均匀增大了整体环空的ECD值，但对增大摩擦压降或摩擦压降梯度无关。因此，理论上该工况下压损比保持为1较为合理。同时，通过分析环空压力沿径向的分布曲线，由转速与井壁压力的关系可知，南海油气井钻进过程中当钻具/杆转速<180 rad/min时，环空沿径向压力并不会发生明显变化，因此更不会引起ECD发生较明显的增加。

（2）对于钻具转速达到600 rad/min时，沿径向的环空压力有明显梯度变化，并且计算中注意到轴向方向压力呈现小幅波动，这一现象分析认为可能与Taylor涡胞现象有关。此外，钻井工程中业界通常认为旋转管柱会对摩阻计算造成影响，产生的原因可能为其他因素导致，例如管柱运动导致的偏心度变化、摆动频率变化等，进而导致流动状态发生变化，但本研究中此现象并不非常明显，后续将进一步开展相关研究。