改进HB模型与液滴动力学耦合的垂直气井流动特性模拟

Simulation of Flow Characteristics in Vertical Gas Wells Coupling Improved HB Model and Droplet Dynamics

摘要: 针对气井积液诊断中井筒压力及微观液滴行为预测的难点，提出一种基于宏观-微观单向耦合的数值模拟方法。采用考虑加速压降项与算术平均粘度的改进Hagedorn-Brown (HB)模型计算垂直井筒气液两相宏观流动参数；基于Hinze湍流破碎理论与Saffman湍流碰撞聚合理论，建立沿井筒高度分布的液滴破碎与聚并动力学模型。利用文献实验数据验证了压力预测模型的准确性，相对误差在3%以内。参数敏感性分析表明：管径增大导致产能非线性剧增，进而显著提高井筒内流体湍流强度，使液滴平均直径有减小的趋势；井壁粗糙度增加会限制产能，由于流速降低导致的湍流耗散减弱效应超过了粗糙度增加带来的湍流增强效应，使液滴直径有增大的趋势；井口回压升高会产生“背压限制”与“气体压缩”双重效应，致使流速大幅下降，使液滴直径有显著的增大趋势。研究揭示了宏观工况通过改变产能进而影响微观流场剪切强度的物理机制。

Abstract: To address the challenges in predicting wellbore pressure and microscopic droplet behavior for liquid loading diagnosis in gas wells, a numerical simulation method based on macro-micro one-way coupling is proposed. An improved Hagedorn-Brown (HB) model, incorporating the acceleration pressure drop term and arithmetic mean viscosity, is adopted to calculate the macroscopic flow parameters of gas-liquid two-phase flow in vertical wellbores. Furthermore, based on the Hinze turbulence breakup theory and the Saffman turbulence collision-coalescence theory, a dynamics model for droplet breakup and coalescence distributed along the wellbore depth is established. The accuracy of the pressure prediction model is validated using experimental data from the literature, with relative errors controlled within 3%. Parameter sensitivity analysis indicates that increasing the tubing diameter leads to a non-linear surge in productivity, which significantly enhances the turbulence intensity of the fluid and consequently causes a decreasing trend in the average droplet diameter. Conversely, increasing wellbore wall roughness limits productivity; the effect of weakened turbulent dissipation caused by reduced flow velocity outweighs the turbulence enhancement effect brought by increased roughness, resulting in an increasing trend in droplet diameter. Additionally, elevated wellhead back pressure induces a dual effect of “back-pressure constraint” and “gas compression”, causing a substantial drop in flow velocity and a significant increase in droplet diameter. This study elucidates the physical mechanism by which macroscopic operating conditions regulate the shear strength of the microscopic flow field by altering gas production capacity.

1. 引言

在天然气开采过程中，随着地层能量的衰竭，井底积液成为制约气井寿命的关键问题。准确预测井筒内的压力梯度及液相存在形式(液膜或液滴)，是判断气井是否满足临界携液条件的前提。

目前，工业界广泛采用Hagedorn-Brown等[1]经验相关式进行压力预测，但这些模型通常将流体视为均质混合物，忽略了微观分散相(液滴)的动力学演变。然而，液滴的直径分布直接决定了其终端沉降速度，是判定气流能否有效携液的关键微观参数。传统的携液模型及其最新的修正研究虽然在一定程度上考虑了液滴受力，但通常假设液滴为固定形状或单一临界尺寸，忽略了井筒沿程流动条件变化对液滴尺寸的动态重塑[2] [3]。Liao等[4]提出了一种考虑液滴与液膜共存的临界携液计算方法，通过分析微观气液界面形态获取液滴尺寸参数，弥补了单一临界流速的不足。

关于气井积液与多相流压降预测的研究，近年来国内外学者在宏观经验模型修正与微观机理探索方面取得了显著进展。在临界携液模型方面，传统的Turner模型因未考虑液滴形变和流型转变而存在预测偏差。刘楠楠等[5]利用实际生产数据对模型系数进行修正，显著提高了临界携液流量的预测准确率，误差降低了47.5%。He等[6]对不同类型的携液模型进行了分类，并分析了不同模型的适用性，指出了大多数模型并未考虑流态对携液能力的影响，在建立新模型时应考虑不同流态下的两相流特性。柯文奇[7]采用多相管流模拟试验，实现在线液膜、液滴测量，明确气井生产过程中液滴、液膜的变化情况。通过系列定量，明确了在气液比233~1463 m³/m³条件下液滴尺寸分布特征。在井筒多相流压降计算方面，Hagedorn-Brown(HB)模型因其鲁棒性仍被广泛应用，但学者们致力于对其进行修正以适应更复杂的工况。Ughulu等[8]提出了一种针对HB相关式的新修正方法，通过引入Beggs-Brill滑脱持液率公式，使模型在水平井和不同流型下的压降预测R²值提升至0.988。此外，机器学习方法也被引入该领域，Salisu等[9]利用非线性回归算法预测凝析气井的积液时机，展示了数据驱动方法的潜力。在微观液滴动力学与多相流耦合方面，研究尚处于起步阶段。董钊等[10]在研究深水气井水合物堵塞时，建立了包含液滴破碎、聚并及沉降行为的动力学模型，证明了考虑微观颗粒演变对宏观流场预测的重要性。国内专利也公开了一种考虑液滴形变和多参数影响的计算方法，试图从受力平衡角度解决积液判断问题[11]。然而，目前大多数研究仍将宏观流场与微观液滴视为静态耦合，缺乏能够沿井筒深度方向实时追踪液滴尺寸演变及其反作用于宏观压降的系统性研究。

文章旨在建立一种宏观流场与微观液滴动力学相结合的数值模拟方法。利用MATLAB编程实现改进的HB模型求解井筒压力场，并以此为边界条件，单向耦合液滴沿井筒上升过程中的破碎与聚并模型。首先利用文献中的实验数据验证模型的可靠性，随后系统研究管径、粗糙度及井口压力变化如何通过影响宏观产能，进而反作用于微观液滴的动力学形态。

2. 数学模型与计算方法

构建的数值计算模型主要包含两个耦合模块：宏观多相流模块与微观液滴动力学模块。研究假设微观液滴尺寸的变化对宏观混合密度及粘度的反作用较小，主要关注宏观流场对微观形态的影响作用。

2.1. 宏观流动模型

为准确描述气液两相在垂直井筒中的流动特性，采用广泛应用于工业界的Hagedorn-Brown经验模型进行求解。基于动量守恒原理，沿垂直井筒深度方向的总压力梯度由重力项、摩擦项和加速项三部分组成：

$\frac{\text{d}P}{\text{d}z}={\rho }_{mix}g+\frac{f{\rho }_{mix}{u}_{mix}^2}{2D}+{\rho }_{mix}{u}_{mix}\frac{\text{d}{u}_{mix}}{\text{d}z}$ (1)

式中，$\text{d}P/\text{d}z$ 为压力梯度，${\rho }_{mix}$ 为气液混合密度，$u_{mix}$ 为气液混合流速，$g$ 为重力加速度，$D$ 为油管内径，$f$ 为摩擦系数。模型显式保留了通常被忽略的加速项，以提高深井高流速工况下的计算精度。

混合密度${\rho }_{mix}$ 取决于滑脱持液率$H_L$ ：

${\rho }_{mix}=H_L{\rho }_l+\left(1-H_L\right){\rho }_g$ (2)

式中${\rho }_l$ 为液相密度，${\rho }_g$ 为气相密度。

滑脱持液率采用Hagedorn-Brown的原始经验相关式计算。该方法引入了四个无量纲数：液体速度数$N_{vL}$ 、气体速度数$N_{vG}$ 、管径数$N_D$ 、和液体粘度数$N_L$ 。通过计算这四个无量纲数，利用Hagedorn-Brown的经验图表求解得到持液率$H_L$ 。

摩擦压降项的计算依赖于两相雷诺数的准确估算：

$R{e}_{mix}=\frac{{\rho }_{mix}{u}_{mix}D}{{\mu }_{mix}}$ (3)

式中${\mu }_{mix}$ 为混合粘度，针对原始Hagedorn-Brown模型中混合粘度采用几何平均导致雷诺数估算偏高的问题，采用了更符合两相流物理机制的算术平均粘度模型，从而提高了摩擦压降预测的可靠性采用基于持液率的加权平均值计算：

${\mu }_{mix}=H_L{\mu }_l+\left(1-H_L\right){\mu }_g$ (4)

摩擦系数$f$ 采用Colebrook-White [12]提出的隐式方程计算。该方程通过将光滑管定律与粗糙管定律在对数项中叠加，从而覆盖了从水力光滑区到完全湍流粗糙区的全过渡区域：

$\frac{1}{\sqrt{f}}=-2{\log }_{10}\left(\frac{k/D}{3.7}+\frac{2.51}{R{e}_{mix}\sqrt{f}}\right)$ (5)

式中$k$ 为绝对粗糙度。

流体物性方面，模拟流体设定为甲烷与水的混合体系。为准确模拟井筒内的压力与温度分布，充分考虑流体物性随压力和温度的动态变化，气相密度采用真实气体状态方程计算：

${\rho }_g=\frac{P{M}_g}{ZRT}$ (6)

式中，$M_g$ 为甲烷摩尔质量，R为通用气体常数，Z为气体偏差因子，P为实际压力，T为实际温度。

气体动力粘度采用Lucas提出的对应态原理模型计算。该模型综合考虑了对比温度和对比压力的影响，能够精确描述高压环境下气体的粘度变化：

${\mu }_g=\frac{Z_1\left(T_r\right)\cdot Z_2\left(T_r,P_r\right)}{\xi }$ (7)

式中，$\xi$ 为Lucas参数，$Z_1\left(T_r\right)$ 和$Z_2\left(T_r,P_r\right)$ 分别为低压项与高压修正因子，具体系数取自Poling等[13]的物性手册。

液相的动力粘度随温度的变化采用Andrade [14]经验方程描述：

${\mu }_l=A\cdot \exp \left(\frac{B}{T}\right)$ (8)

式中，A，B为公式系数基于水的标准粘度数据回归得到。

相较于气体物性的显著变化，文章忽略了液体密度和表面张力随压力和温度的波动，采用了水–气体系的典型值，分别取值为1000 kg/m³和0.07 N/m。

本研究采用宏观流场向微观动力学的单向耦合求解策略。虽然液滴尺寸的变化理论上会通过改变相间曳力进而影响滑脱持液率，但在低液气比的垂直气井中，气体携带液滴的体积分数较低，液滴尺寸微小的变化对混合流体宏观物性的非线性反作用较弱。因此，本文主要聚焦于宏观工况参数对微观流场结构的影响机制，暂时忽略微观参数对宏观压降的迭代反馈，这在保证计算效率的同时能够满足工程对趋势预测的需求。

2.2. 微观液滴动力学模型

假设液滴在湍流气核中运动，其尺寸演变受湍流剪切破碎和碰撞聚并影响。液滴在高速气流中的破碎行为主要受控于外部气动压力与内部表面张力的相互作用。文章基于Hinze [15]提出的经典湍流破碎理论，采用韦伯数作为判断液滴是否发生破碎的关键无量纲参数。韦伯数定义为惯性力与表面张力之比：

$We=\frac{{\rho }_g{u}_{rel}^{2}d}{\sigma }$ (9)

式中，${\rho }_g$ 为气体密度，$u_{rel}$ 为气液两相相对速度，d为液滴直径，$\sigma$ 为气液界面张力。

根据Pilch和Erdman [16]对气液两相流破碎机理的研究，液滴的破碎模式随韦伯数变化。当We超过临界值时，液滴失稳并发生破碎。参考Pilch等[16]的研究，对于低粘度液体在气流中的剪切破碎，取破碎时的临界韦伯数$W{e}_{break}=12$ 。由于井筒尺度巨大，直接耦合复杂的湍流DNS模型计算成本过高，并且考虑到井筒内湍流流动的随机性，文章并未假设当时液滴必然立即破碎，而是引入了基于过余韦伯数的简化线性随机概率模型。破碎概率$P_{break}$ 定义为：

$P_{break}=\min \left(\lambda \frac{We-W{e}_{break}}{W{e}_{break}},P_{\max }\right)$ (10)

式中，$\lambda$ 为破碎概率系数，$P_{\max }$ 为最大破碎概率。

该公式反映了驱动力越强，液滴发生破碎的可能性越大这一物理事实。当判定破碎发生时，为满足流动稳定性要求，生成的子液滴直径$d_{new}$ 依据能量守恒与最大稳定直径准则进行更新，假设破碎过程并非一次完成，而是通过级联过程逐渐减小，单次破碎后的直径按韦伯数比率的立方根缩放，以避免数值计算中的直径突变导致的不稳定性。该处理方法确保了破碎后的子液滴在该流场条件下的动力学稳定性。

$d_{new}=d_{current}{\left(\frac{W{e}_{break}}{We}\right)}^{1/3}$ (11)

式中$d_{current}$ 和$d_{new}$ 分别表示破碎前后的液滴直径。

在垂直气井的高雷诺数流动中，湍流脉动是导致液滴相互碰撞并发生聚并的主要驱动力。文章采用基于群体平衡方程(Population Balance Equation)的简化模型来描述液滴数密度的演变。假设液滴在控制体内空间分布均匀，对Smoluchowski [17]碰撞频率方程关于全尺寸范围进行积分，可得液滴总数密度n随时间的演化遵循二阶动力学衰减规律：

$\frac{\text{d}n}{\text{d}t}=-K_{eff}{n}^2$ (12)

式中，n为液滴数密度，t为时间，$K_{eff}$ 为有效聚并速率常数。

根据Saffman和Turner的湍流碰撞理论[18]，碰撞频率，即有效聚并速率常数主要取决于流体的湍流能量耗散率和液滴尺寸，具体形式如下：

$K_{eff}=\alpha {\epsilon }^{-1/3}{d}^{7/3}$ (13)

式中，$\alpha$ 为聚并效率系数，d为液滴直径，$\epsilon$ 为单位质量的湍流能量耗散率。

井筒内流体的平均湍流能量耗散率基于宏观能量平衡计算。根据Levich [19]在湍流管流理论中阐述的摩擦阻力定律，单位质量流体的能量耗散率可由宏观管流摩阻系数$f$ 和气体流速$u_g$ 计算得出：

$\epsilon =\frac{f{u}_g^3}{2D}$ (14)

2.3. 求解策略

采用自上而下的有限差分法。基于已知井口压力和温度等边界条件，采用迭代法沿井深方向离散化求解压力分布。随后，基于计算得到的流场数据，从井底向井口追踪液滴群的动力学演变。

3. 模型验证

为确保耦合数值模型的可靠性，本文采取先宏观流场、后微观形态的分步验证策略。准确的井筒压力及流速分布预测是液滴动力学计算的前提，因此，首先利用文献中的垂直井实测数据对改进HB宏观流动模型进行精度检验；在此基础上，结合经典的Turner [20]微观理论，对液滴动力学模型的计算结果进行量级与趋势的合理性分析，从而全面评估耦合模型在工程应用中的有效性。

3.1. 宏观参数合理性验证

为验证代码中压降计算模块的准确性，文章选取Reinicke [21]文献中的垂直井实测数据进行验证。表中及后续文中的体积流量均是指在标准状况下的体积流量。计算结果如表1所示，改进的HB模型在典型工况下的相对误差在3%以内，证明了宏观压力求解器的可靠性。

Table 1. Comparison between measured data and simulation results

表1. 实测数据与模拟结果对比

工况序号	气液比/10^4 m^3/m^3	井口压力/MPa	实测井底压力	模拟井底压力	相对误差
1	20	16.9	22.03	21.56	−2.1
2	25	19.2	25.21	24.54	−2.66
3	1.85	12.2	17.9	17.46	−2.46
4	2.38	12.1	17.7	17.24	−2.6

3.2. 微观液滴尺寸合理性验证

由于高压气井井下液滴尺寸实测数据极难获取，本文采用经典的临界韦伯数理论对模拟结果的量级合理性进行验证。根据Turner等提出的液滴携带理论，在给定工况下，流场中能够稳定存在的临界液滴直径满足以下关系式：

$d_{crit}=\frac{W{e}_{crit}\cdot \sigma }{{\rho }_g{u}_g^2}$ (15)

式中，$W{e}_{crit}$ 为Turner理论下对应的临界韦伯数，取值为30。

将基准工况参数代入计算，得到模拟液滴直径与理论临界液滴直径的关系如图1所示，理论最大的稳定直径约为4.7 mm。本文模拟得到的液滴直径$d_{sim}$ 分布在1.9~2.4 mm之间，与理论估算值处于同一

Figure 1. Relationship curves of simulated droplet diameter and theoretical critical droplet diameter versus well depth

图1. 模拟液滴直径与理论临界液滴直径与井深的关系曲线

数量级，且随井深z的变化趋势与理论值和柯文奇[7]在实验中观测到的规律一致，说明微观模型具有物理合理性。

4. 结果与讨论

基于验证后的模型，文章采用控制变量法的思路，研究井筒内径D、管壁绝对粗糙度k和井口压力P_t对气井产能及微观液滴特性的影响。根据Reinicke [21]文献选取模拟基准工况为：井深3053 m，井筒内径0.1 m，气体流量264,000 m³/d，液体流量1.32 m³/d，井口压力16.9 MPa，井口温度54℃，井底压力22.03 MPa，井底温度106℃，由于文献中没有给出具体的管壁粗糙度，文章先假设管壁绝对粗糙度为0.1 mm。

4.1. 井筒内径对产能和液滴直径的影响

将其他边界条件设置为基准工况，改变井筒内径D，分析井筒内径对井筒产能和整个井筒中的平均液滴直径的影响。

Figure 2. Relationship curves between wellbore inner diameter and macroscopic productivity

图2. 井筒内径与宏观产能关系曲线

宏观产能分析：如图2所示，随着井筒内径的增大，气井的标准状态下产气量和产液量均呈现剧烈的非线性增长。当管径从0.04 m增至0.12 m时，产气量增加了近10倍。由管流流体力学理论可知，摩擦压降与管径的五次方成反比，管径的微幅增加即可显著降低流动阻力，从而大幅提升气井携液能力。在驱动压差一定的情况下，系统需通过大幅增加流量来达到新的压力平衡。

微观液滴分析：模拟结果如图3所示，整个井筒内液滴的最大和最小平均直径均随井筒内径的增大而显著减小。值得注意的是，该结果与“大管径导致低流速”的常规认知存在差异，但在产能受限的节点分析中，流量的剧增起到了主导作用。尽管流通面积随井筒内径增加，但流量的增长速率远超面积增长，导致实际气体流速不降反升。根据Levich[19]理论，湍流耗散率与气体流速的三次方成正比。流速的增加导致湍流剪切力急剧增强，巨大的气动韦伯数促使液滴发生剧烈破碎，从而形成了更小的微观粒径。

进一步结合临界韦伯数对液滴的悬浮稳定性进行判定：如图4所示，当井筒内径为0.1 m时，全井段内的气动韦伯数始终维持在13~18之间，均高于临界破碎韦伯数。从流体动力学角度分析，这意味着液滴处于剪切破碎主导的动力学不稳定区，大尺寸液滴无法在流场中稳定存在，会被高速气流持续剪切破碎成细小的雾滴。然而，正是这种强烈的破碎机制，使得液滴直径始终小于其在当前流速下的临界悬浮直径。因此，虽然液滴尺寸在微观上表现为动态破碎，但在宏观运动状态上，细小的液滴能够满足临界携液条件，处于稳定的气相悬浮输运状态，说明此时井底积液的风险较低。

Figure 3. Relationship curves between wellbore inner diameter and droplet diameter

图3. 井筒内径与液滴直径关系曲线

Figure 4. Relationship curve between Weber number and well depth

图4. 韦伯数与井深的关系曲线

4.2. 井壁粗糙度对产能和液滴直径的影响

将其他边界条件设置为基准工况，改变绝对粗糙度k，分析绝对粗糙度对井筒产能和整个井筒中的平均液滴直径的影响。

宏观产能分析：如图5所示，随着粗糙度的增加，产液产气量逐渐减小，且在低粗糙度区间下降较快，随后趋于平缓。粗糙度的增加直接提高了Colebrook摩擦系数，导致沿程阻力损失增加。为维持井底到井口的压力平衡，气井只能通过降低产量来减少摩阻损失，达到新的平衡。

微观液滴分析：如图6所示，整个井筒内液滴的平均直径随粗糙度的增加呈现明显的增大趋势。这里存在两种竞争机制：一方面，粗糙度增加会增大摩擦系数，理论上应增强湍流，促进液滴破碎；但另一方面，产能的下降导致气体流速降低，不利于液滴破碎。由于湍流耗散率对流速是三次方敏感，流速降低带来的剪切减弱效应远超摩擦系数增加带来的扰动增强效应。最终结果是井筒内湍流强度下降，液滴难以破碎，导致平均直径增大。

Figure 5. Relationship curves between absolute roughness and macroscopic productivity

图5. 绝对粗糙度与宏观产能关系曲线

Figure 6. Relationship curves between absolute roughness and droplet diameter

图6. 绝对粗糙度与液滴直径关系曲线

4.3. 井口压力对产能和液滴直径的影响

将其他边界条件设置为基准工况，改变井口压力P_t，分析井口压力对井筒产能和整个井筒中的平均液滴直径的影响。

宏观产能分析：如图7所示，随着井口压力的微小增加，气井产能呈现急速下降趋势。提高井口压力直接增加了系统的背压，在井底流压受地层供液能力限制的情况下，有效举升压差被显著压缩，导致流动能力大幅衰减。

Figure 7. Relationship curves between wellhead pressure and macroscopic productivity

图7. 井口压力与宏观产能关系曲线

Figure 8. Relationship curves between wellhead pressure and droplet diameter

图8. 井口压力与液滴直径关系曲线

微观液滴分析：模拟结果如图8所示，整个井筒内液滴直径随井口压力的升高而急剧增大。首先，产能下降导致气相质量流量降低；其次，高压导致气体密度线性增加。在质量流量受限条件下，显著降低了气相的表观流速。流速的降低直接导致井筒内湍流能量耗散率呈指数级衰减。由于湍流涡旋强度不足，无法提供足够的剪切能量克服液滴表面张力，导致液滴破碎过程受抑，同时低湍流强度加剧了液滴聚并，最终导致宏观上液滴直径显著增大。

5. 结论

文章通过构建宏观Hagedorn-Brown流动与微观Hinze液滴动力学的耦合模型，对垂直气井的流动规律进行了深入研究。主要结论如下：

(1) 宏观工况参数主要通过改变流动阻力和有效压差来影响气井产能。在其他参数不变的条件下，管径增大能大幅释放产能，而粗糙度增加和井口回压升高则显著抑制产能。

(2) 井筒内的液滴直径并非定值，而是由宏观流速决定的湍流耗散率影响。研究发现，产气能力决定了微观流场的剪切强度：大管径工况带来的产能剧增效应显著提高了流体流速与湍流强度，促使液滴发生剧烈破碎进而减小；反之，高井壁粗糙度或高井口回压工况则显著限制了气井产能，导致井筒内流速降低及湍流耗散减弱，气相剪切力不足以克服液滴表面张力，从而使得液滴维持较大的平均直径。

(3) 在进行深层气井积液诊断时，必须动态考虑液滴尺寸的变化。对于高压、低产或高粗糙度的气井，由于液滴直径较大，其所需的临界携液流速可能与传统Turner公式的预测值存在偏差，建议在工艺设计中予以修正。

基金项目

国家自然科学基金“基于流–声耦合作用物理机制的螺旋桨远场噪声数值预报方法研究”(基金号52471351)。