引 言

气井积液是气井生产中需要克服的重要课题之一。研究发现，气液流体以螺旋状态向上运动可以获得最小的压力损失^[1]。在这一理论基础上，气井井下涡流工具被设计出来。涡流工具将自发的气液流动状态转变为具有一定稳定性的螺旋环状流，从而增加气井产量，提高携液能力。Ali等首先进行了涡流工具的实验研究^[2]，指出涡流可以降低流动压降和临界携液流速。Hein介绍了2002——2006年美国能源部进行涡流工具现场试验^[3]。Surendra等对下入涡流工具后的两相流动过程进行了数值模拟研究^[4]。随后，国内也逐渐进行了大量涡流工具的现场试验^[5-8]，开展了对于涡流工具的实验及理论研究^[9-12]。吴晓东等根据物理实验数据给出了旋流条件下气体临界携液模型^[13]。石怡帆等使用Fluent模拟了不同流型下，气液经过涡流工具前后的流动过程^[14]。张昭等使用CFD方法模拟了不同涡流工具结构下的流动过程，分析了涡流工具叶片数量和分布方式对携液性能的影响^[15]。

由前期的现场试验也可以得到，涡流排液方法存在适用条件以及最优工作状态。只有通过合理的选井以及设计才能发挥涡流工具的排液效果。由于目前对于涡流工具排液机理的研究和认识均不够全面，没有建立针对涡流工具工作特性以及流动特性的理论模型^[16-19]，限制了涡流工具的应用。因此，在前人的多相流研究中，没有针对于旋流流态的两相流动模型。本研究开展了涡流工具的气液两相流动实验，并在两相流动理论和实验结果的基础上，建立旋转环空流动条件下，气液两相流动参数计算的机理模型。

1 涡流工具气液两相涡流实验

气液涡流实验装置如图 1所示。实验中使用水和空气分别作为液相和气相流动介质。其中，液体通过螺杆泵加压，最大流量为6.3 m³/h。注气系统使用空压机最大排气量为1.2 m³/min，最大的出口压力为1.2 MPa。流动实验环路中安装了垂直测试管道，垂直管道下部可安装和更换不同结构参数的涡流工具。测试管道总长度为7.5 m，由两段3.0 m的不锈钢管和3段0.5 m的可视有机玻璃管组成。实验中加工了3个具有不同螺旋角的涡流工具以模拟不同旋流条件，如图 2a所示。工具1，2，3的螺旋角分别为30°，45°和55°。

如图 1所示，该实验平台可以研究气液两相旋流规律以及涡流工具的工作特征。在涡流工具前后以及管道入口和出口处安置了压力表，分别用来测量气液流经工具产生的节流损失和在垂直管道中由于重力和摩擦产生的压降。在测试段的供气和供液通道上分别安装有气体流量计和液体转子流量计。在管道初始段、中间段和结束段加入了可视段。

如图 2所示，流经涡流工具后，可视段的管壁上可以观测到螺旋状的液环，根据测量液环痕迹的角度可以测出液相旋流强度。通过采集3个不同位置可视段的流动图片，可以观测到由于液膜与管壁摩擦导致的液相旋流速度下降，即液环与水平面的夹角沿着流动方向逐渐增大。

2 流动模型建立 2.1 涡流水力学参数预测模型

气液相流入涡流工具后，液体汇集到螺旋流道的上部，因此，在流出工具后在管壁上形成如图 2b中的液环，一部分液体由于重力作用从液环下落形成液膜。最终液环重力与气体对其的拖拽作用达到平衡。与其他两相流模型一样，流型是研究两相流动机理的基础^[20]。图 2b是对旋流状态下气液环空流流型的简化。由于旋流的存在，旋流模型不仅需要考虑轴向上的水力学机理，对于横截面上周向和径向的运动也需要给出相应的控制方程。

r方向(径向)上，假设气芯和液膜的径向速度为0，由于旋转速度的存在，向心加速度造成管道横截面上，外部的液相压力要高于位于内部的气芯压力。在柱坐标下的N-S径向方程中，令径向速度u_r=0，有

$ \frac{{u_\varphi ^2}}{r}{\rm{ = }}\frac{1}{\rho }\frac{{\partial p}}{{\partial r}} $

(1)

定义ω_g和ω_L分别为气芯和液膜的角速度，假设在同一管道横截面上，气液相角速度为定值；p_g和p_L分别为气芯和液膜的平均压力。将式(1)积分，可得

$ {p_{\rm{L}}} - {p_{\rm{g}}} = \frac{1}{2}\left( {{\rho _{\rm{g}}}\omega _{\rm{g}}^2\int_0^{D/2 - 2 - \delta } {r{\rm{d}}r} + {\rho _{\rm{L}}}\omega _{\rm{L}}^2\int_{D/2 - \delta }^{D/2} {r{\rm{d}}r} } \right) = \frac{{{D^2}}}{{16}}\left[ {{\rho _{\rm{g}}}\omega _{\rm{g}}^2{{\left( {1 - 2{\delta _{\rm{D}}}} \right)}^2} + {\rho _{\rm{L}}}\omega _{\rm{L}}^24\left( {{\delta _{\rm{D}}} - \delta _{\rm{D}}^2} \right)} \right] $

(2)

其中，δ_D=δ/D，并假设在单元段内液膜厚度不发生变化。

z方向(轴向)上，可看作气液界面存在波动的环空-雾状流。因此，轴向上气芯和液膜的动量方程分别为^[21-22]

$ - {A_{\rm{g}}}\frac{{{\rm{d}}{p_{\rm{g}}}}}{{{\rm{d}}z}} - {\rho _{\rm{g}}}{\rm{g}}{A_{\rm{g}}} - {\tau _{z{\rm{i}}}}{S_{\rm{f}}} = 0 $

(3)

$ - {A_{\rm{L}}}\frac{{{\rm{d}}{p_{\rm{L}}}}}{{{\rm{d}}z}} - {\rho _{\rm{L}}}{\rm{g}}{A_{\rm{L}}} - {\tau _{z{\rm{w}}}}{S_{\rm{w}}} + {\tau _{z{\rm{i}}}}{S_{\rm{f}}} = 0 $

(4)

φ方向(周向)上，考虑界面上的角动量平衡，即横截面上，气相角动量沿z方向的减少量用以克服气液界面上摩擦力产生的力矩。同样地，液相角动量的轴向变化量为相间摩擦力矩以及与管壁的摩擦力矩共同作用的结果。因此，可以得到

$ \frac{{\rm{d}}}{{{\rm{d}}z}}\int_0^{D/2 - \delta } {{\rho _{\rm{g}}}{u_{z{\rm{g}}}}2\pi r{\rm{d}}r\left( {{r^2}{\omega _{\rm{g}}}} \right)} = - {\tau _{\varphi {\rm{i}}}}{S_{\rm{f}}}\frac{D}{2}\left( {1 - 2{\delta _{\rm{D}}}} \right) $

(5)

$ \frac{{\rm{d}}}{{{\rm{d}}z}}\int_{D/2 - \delta }^{D/2} {{\rho _{\rm{L}}}{u_{z{\rm{L}}}}2\pi r{\rm{d}}r\left( {{r^2}{\omega _{\rm{L}}}} \right)} = - {\tau _{\varphi {\rm{w}}}}{S_{\rm{w}}}\frac{D}{2} + {\tau _{\varphi {\rm{i}}}}{S_{\rm{f}}}\frac{D}{2}\left( {1 - 2{\delta _{\rm{D}}}} \right) $

(6)

定义气相和液相的旋流强度W_g和W_L

$ {W_{\rm{g}}}{\rm{ = }}\frac{{D\left( {1 - 2{\delta _{\rm{D}}}} \right){\omega _{\rm{g}}}}}{{2{u_{z{\rm{g}}}}}} $

(7)

$ {W_{\rm{L}}}{\rm{ = }}\frac{{D{\omega _{\rm{L}}}}}{{2{u_{z{\rm{L}}}}}} $

(8)

定义无因次变量 Γ_g和 Γ_L。W_g和W_L定义与文献[22]中相同

$ {\mathit{\Gamma }_{\rm{g}}}{\rm{ = }}\frac{{{\tau _{\varphi {\rm{i}}}}{S_{\rm{f}}}}}{{{A_{\rm{g}}}}} $

(9)

$ {\mathit{\Gamma }_{\rm{L}}}{\rm{ = }}\frac{{{\tau _{\varphi {\rm{w}}}}{S_{\rm{w}}} - {\tau _{\varphi {\rm{i}}}}{S_{\rm{f}}}\left( {1 - 2{\delta _{\rm{D}}}} \right)}}{{{A_{\rm{L}}}}} $

(10)

由式(5)，式(6)得到

$ \frac{{{\rm{d}}{W_{\rm{g}}}}}{{{\rm{d}}z}}{\rm{ = }}\frac{{ - 2{\mathit{\Gamma }_{\rm{g}}}}}{{{\rho _{\rm{g}}}u_{z{\rm{g}}}^2}} $

(11)

$ \frac{{{\rm{d}}{W_{\rm{L}}}}}{{{\rm{d}}z}}{\rm{ = }}\frac{{ - {\mathit{\Gamma }_{\rm{L}}}}}{{{\rho _{\rm{L}}}u_{z{\rm{L}}}^2}}\frac{1}{{\left( {1 - 2{\delta _{\rm{D}}}{\rm{ + }}2\delta _{\rm{D}}^2} \right)}} $

(12)

2.2 求解液膜平均厚度 δ

由式(2)，式(3)和式(4)可以消去 p_g和 p_L，得到

$ - \frac{{{\tau _{z{\rm{w}}}}{S_{\rm{w}}}}}{{{A_{\rm{L}}}}} + {\tau _{z{\rm{i}}}}{S_{\rm{f}}}\left( {\frac{1}{{{A_{\rm{g}}}}} + \frac{1}{{{A_{\rm{L}}}}}} \right) - \left( {{\rho _{\rm{L}}} - {\rho _{\rm{g}}}} \right){\rm{g}} - W = 0 $

(13)

根据式(13)可以求解液膜平均厚度 δ。其中，变量W反映了旋流产生的离心力对 δ的影响，若W=0，则式(13)便与环空-雾状流模型相同^[23-24]。螺旋流随着距离会慢慢消散，即dω/dz<0。由式(14)可知，W<0，这样式(13)得到的液膜厚度要小于W=0的情况，因此，旋流可以使得重力压力损失下降，这是涡流工具的工作机理之一。W的计算如下

$ W{\rm{ = }} - {\mathit{\Gamma }_{\rm{g}}}{W_{\rm{g}}} - \frac{{2\left( {{\delta _{\rm{D}}} - \delta _{\rm{D}}^2} \right)}}{{1 - 2{\delta _{\rm{D}}} + 2\delta _{\rm{D}}^2}}{\mathit{\Gamma }_{\rm{L}}}{W_{\rm{L}}} $

(14)

2.3 气液相间作用

式(3)~式(14)中轴向剪切应力与周向剪切应力分别为

$ {\tau _{z{\rm{w}}}}{\rm{ = }}0.5{f_{z{\rm{w}}}}{\rho _{\rm{L}}}u_{z{\rm{L}}}^2 $

(15)

$ {\tau _{\varphi {\rm{w}}}}{\rm{ = }}0.5{f_{\varphi {\rm{w}}}}{\rho _{\rm{L}}}{\left( {{W_{\rm{L}}}{u_{z{\rm{L}}}}} \right)^2} $

(16)

$ {\tau _{z{\rm{i}}}}{\rm{ = }}{\tau _{{\rm{gf}}}}\left[ {1{\rm{ + }}\frac{{2\left( {\eta - 1} \right){h_{{\rm{aD}}}}}}{\pi }} \right] $

(17)

$ {\tau _{\varphi {\rm{i}}}}{\rm{ = 0}}.{\rm{5}}{f_{\varphi {\rm{i}}}}{\rho _{\rm{g}}}{\left( {{W_{\rm{g}}}{u_{z{\rm{g}}}}} \right)^2} $

(18)

其中，τ_zw，τ_φw，τ_zi和τ_φi分别为液膜与管壁在轴向和周向上的剪切应力，液相与气相之间在轴向和周向上的剪切应力。由于涡流工具的作用，液流在工具内部汇集进入流动管道，并沿着螺旋上升的液环流道流动。凸出的液环进入并占据了原本气芯的过流截面，气液界面的黏性剪切作用变为气芯流过液环产生的拖拽力。为了与环空流动模型产生联系，将拖拽力等效为剪切力即可得到式(17)。

计算结果显示，η >1和τ_zi ≥τ_gf，因此，在涡流工具的作用下，气相对液相向上的作用力增加，最终根据式(13)计算得到液膜厚度下降，这也使得存在旋流速度时的重力压降下降。

式(15)~式(18)中摩擦系数计算如下。其中，气液界面的轴向摩擦系数f_zi使用文献[15]中的关联式，如式(20)所示，其中，气相速度取u_g=u_zg(1+W_g)^0.5。

$ {f_{z{\rm{w}}}} = 0.0460Re_{\rm{L}}^{ - 0.2000} $

(19)

$ {f_{z{\rm{i}}}} = 0.2400{f_{\rm{g}}}{\rm{ }}Re_{\rm{L}}^{0.3050}{\left[ {\frac{\sigma }{{D{\rho _{\rm{g}}}u_{\rm{g}}^2\left( {1 - 2{\delta _{\rm{D}}}} \right)}}} \right]^{0.085}} $

(20)

$ {f_{\varphi {\rm{w}}}} = 0.0175{f_{\rm{w}}}W_{\rm{L}}^{ - 0.2940} $

(21)

$ {f_{\varphi {\rm{i}}}}{\rm{ = }}0.0175{f_{\rm{g}}}W_{\rm{g}}^{ - 0.2940} $

(22)

2.4 模型的求解

模型的求解需要其他辅助方程。由物质平衡可知，管道沿程的气液相质量流量不会发生变化。假设气液相的表观速度均已知，这样，气液相的平均轴向速度u_zg和u_zL为

$ {u_{z{\rm{g}}}} = \frac{{{A_{\rm{t}}}}}{{{A_{\rm{g}}}}}{u_{{\rm{sg}}}} $

(23)

$ {u_{z{\rm{L}}}} = \frac{{{A_{\rm{t}}}}}{{{A_{\rm{L}}}}}{u_{{\rm{sL}}}} $

(24)

几何参数为

$ {A_{\rm{t}}} = \frac{1}{4}\pi {D^2} $

(25)

$ {A_{\rm{L}}}{\rm{ = }}{A_{\rm{f}}} + {A_{\rm{r}}}{\rm{ = }}\pi {D^2}\left( {{\delta _{\rm{D}}} - \delta _{\rm{D}}^2} \right){\rm{ + }}\pi h_{\rm{a}}^2 $

(26)

$ {A_{\rm{g}}}{\rm{ = }}{A_{\rm{t}}} - {A_{\rm{L}}}{\rm{ = }}\frac{1}{4}\pi {D^2}{\left( {1 - 2{\delta _{\rm{D}}}} \right)^2} - {\rm{ }}\pi h_{\rm{a}}^2 $

(27)

$ {S_{\rm{f}}}{\rm{ = }}\pi D\left( {1 - 2{\delta _{\rm{D}}}} \right) $

(28)

$ {S_{\rm{w}}}{\rm{ = }}\pi D $

(29)

并有管道截面持液率

$ {H_{\rm{L}}}{\rm{ = 4}}\left( {{\delta _{\rm{D}}} - \delta _{\rm{D}}^2} \right){\rm{ + }}h_{{\rm{aD}}}^2{\left( {1 - 2{\delta _{\rm{D}}}} \right)^2} $

(30)

目前计算长管道沿程压力分布的方法均需对长管道进行分段，再选取合理的压力计算模型来计算单元管段中的压力梯度，从而根据入口压力和单元段长度得到出口压力^[25-26]。

3 结果与讨论 3.1 模型结果与已有实验对比

Ali等研究了空气在涡流工具下入前后的临界携液能力，并且给出了临界状态下气液两相流产生的压力降^[2]，本文选取了环空流阶段的实验结果，如图 3a~图 3c中散点所示，可以看出，在相同气液流速下，下入工具后压降要比未下工具前低，入口压力为170.27，239.22和308.17~kPa时，压降平均下降幅度分别为17.0%，8.5%和5.0%。图中还给出了使用Beggs-Brill(B & B)模型，Hagedorn-Brown(H & B)模型和Ansari模型计算正常气液垂直管流的压降，以及本文中旋流模型的计算压降。由于文献[2]中未给出流体的旋流强度，因此，回归得到W_L的取值结果，这里计算旋流压降时W_L取0.88。在入口压力为170.27~kPa时，各个模型计算结果与文献[2]中的实验结果误差较大，而且压降变化趋势均有不同。在入口压力为239.22~kPa时，各个模型计算的压降变化趋势与实验结果趋势吻合度升高。而入口压力为308.17~kPa时，B & B和H & B模型计算结果偏小，Ansari模型与旋流模型的计算结果与实验结果误差很小。整体而言，根据文献[2]中的实验数据，Ansari模型计算的环空气液垂直管流压降准确度较高。在3种入口压力条件下，本文建立的旋流模型结果与Ansari模型的变化趋势保持一致，压降大小比Ansari结果平均分别低22.87%，10.42%和7.83%，这与实验规律基本符合，说明涡流工具降低流动压降的机理在旋流模型中得到了较好的体现。

3种入口压力条件下，旋流模型计算结果与实验结果的最大相对误差分别为30.06%，8.10%和1.76%，平均相对误差分别为11.54%，5.41%和1.08%。入口压力为170.27~kPa时，高流速的计算结果误差较大，但是平均误差在可接受范围内，并且入口压力为239.22~kPa和308.17~kPa的两种情况下，旋流模型计算误差非常小。

3.2 旋流压降实验和模型计算结果

在3.1节中，给出的是Ali等的实验结果，但并没有给出与旋流相关参数的测量，而旋流强度变化可以反映涡流工具的作用距离。并且，Ali等重点研究了涡流工具对于临界携液能力的提升，在文献[2]中未给出超临界携液状态的流动实验结果。为了进一步验证旋流模型，设计了更加详细的实验。本文主要研究对象为旋流状态下的气液环空流，设置的气相速度要高于临界携液流速，具体实验条件对比如表 1所示。

如图 2b所示，从实验可视段可以测量液环的螺旋角度θ_s，根据θ_s可以得到液相的旋流强度，结果如下

$ \frac{1}{{\tan {\theta _{\rm{s}}}}} = \frac{{{u_{\varphi {\rm{L}}}}}}{{{u_{z{\rm{L}}}}}} = \frac{{D{\omega _{\rm{L}}}}}{{2{u_{z{\rm{L}}}}}}{\rm{ = }}{W_{\rm{L}}} $

(31)

在文献[2]中，缺乏相应液相旋流强度参数的数据。而在本文实验中，涡流工具的螺旋角(30°，45°和55°)对应的旋流强度分别为1.73，1.00和0.70。图 4给出了在安装涡流工具1时，不同气液速度条件下的垂直管道压降，其中散点为实验结果，曲线为模型计算结果。随着气流速度增加，一方面会降低截面持液率，降低流动压降；另一方面更大的气流速度会增加流动摩擦阻力，如图所示，液相表观速度为0.05 m/s时，压降从18.21~kPa下降到7.51~kPa后，由于气相速度的增大加剧了摩擦阻力，因此随着气相速度增大而上升到10.51~kPa。在高液流速度下，气相增加降低持液率的作用更为明显，因此垂直管流压差一直呈下降趋势。图中给出了实验点压差结果上下20%的误差线，可以看出，模型计算结果相对于实验结果偏大，但各测试点上的误差基本都在可控范围内。

图 5给出了液相速度为0.05 m/s时，不同涡流工具下气液流动压降随气相速度的变化关系。

图 5中模型计算结果显示，低速和高速条件下，涡流工具对气液流动的影响存在相反的规律：在低速下(气相表观速度为12 m/s和14 m/s)时，压降随着工具螺旋角的增大而增大；而高气相速度(16 m/s和18 m/s)时，压降随工具螺旋角增大而减小。工具1的螺旋角最小，由式(31)可知，其产生的旋流强度更大，气液界面的作用力更大，在低速时有利于降低持液率，但在高速时会增大摩擦阻力。另外，不同工具下气液流动压降的最小值可以比最大值低15%~35%，因此，对于涡流工具的使用，需要根据气井实际情况来对工具结构参数进行优化。

3.3 旋流强度衰减速度实验和模型计算结果

气液相的角速度会因为与管壁以及气液之间的摩擦而逐渐衰减。气相旋流强度无法测量，液相旋流可以通过管壁液环的螺旋角测得，然后根据出口与入口处旋流强度的差值算出实验管段内旋流强度梯度dW_L/dz的平均值。图 6给出了不同条件下，液相旋流强度的理论计算值和实验值。实验结果与计算结果均显示旋流强度沿管道的衰减速度随液相速度增大而减小，这与式(12)相符。模型中气相速度与旋流强度衰减速度关系并不明显，但从图 6中可以看出，-dW_L/dz随气相速度增加呈上升趋势，这是由于气体速度增加使得液膜厚度δ_D减小，由式(12)计算得到的旋流强度的衰减速度会增大。%如图 6所示，低速时，理论计算结果与实验结果相对误差较大。高速条件下计算结果与实验结果符合度较高，说明此时旋流流态充分发育，与模型假设条件相似度更高。

根据图 6中的绿色实线和虚线可以对比工具1和工具3的旋流强度变化，可知初始旋流速度越低，液相与管壁的摩擦越小，这样也就导致旋流强度梯度越小。可以推断，在实际气井下入井下涡流工具后，流体旋流强度沿井筒向上逐渐降低，但是下降的速度会逐渐变缓。%因此，如图 6中，液相速度为0.05 m/s，气体速度为16 m/s时，旋流强度梯度为0.10 m^-1，工具1初始旋流强度为1.73，其旋流维持的距离要大于1.73/0.10=17.3 m。当然在旋流强度降为0之前，旋流降低压降的作用可能就不明显了。因此，需要根据气井实际生产状况设计涡流工具的级数以及每级之间的间隔。

4 结论

(1) 旋流条件下，气液环空流动由于周向运动产生的离心力造成内部气芯与管壁液膜之间存在压力差，并且气液分布特点增加了气芯对液流的拖拽力。两方面的作用可以降低截面持液率，从而在旋流状态下垂直气液两相流动压降减小。结果表明，压降减小幅度在5%~20%。

(2) 根据旋流条件下气液环空流的水力学原理，分别从轴向、周向和径向上给出了气液相的动量方程以及周向压力梯度计算方法。通过与前人实验和模型的对比，表明所建立的模型预测的旋流条件下气液环空流压降具有较高的精度，并且解释了涡流工具改善气液环空流流动状态的机理。

(3) 不同涡流工具，使得气液两相流动具有不同的旋流强度。在高气流速度条件下，更高的旋流强度会增加流动的摩擦阻力；而低速时，更高的旋流速度更有利于降低重力压降。

(4) 实验与模型计算结果表明：在研究的气液流速范围内，液相旋流强度的衰减速度随液流速度的增大而减小，随气相速度增大而增大。

  符 号 说 明

u——速度，m/s；

ρ——密度，kg/m³；

p——压力，Pa；

ω——管道横截面上的平均角速度，rad/s；

D——直径，m；

δ——液膜厚度，m；

z——管道长度，m；

A——面积，m²；

g——重力加速度，g=9.8 m/s²；

s——周长，m；

τ——剪切应力，Pa；

W——旋流强度，无因次；

f——范宁摩擦系数，无因次；

η——环空面积系数，无因次；

h——液环高度，m；

Re——液相雷诺数，无因次；

σ——表面张力，N/m；

H——截面持液率，无因次；

θ_s——液环与水平面所成夹角，（°）。

下标：

a——液环最高幅度；

D——无因次变量；

f——液膜；

i——中心气芯与液膜之间的界面；

r——液环；

s——表观速度；

t——实验管道；

w——液膜与管道的接触面；

g——气相；

L——液相；

r，z，φ——径向、轴向和周向上的变量。