基于群 体平衡理论的管道内水合物颗粒聚集过程模拟

宋光春，李玉星，王武昌，姚淑鹏，魏 丁，房祥鹏

（1.山东省油气储运安全省级重点实验室 中国石油大学（华东），山东 青岛 266580；2.泰能天然气有限公司，山东 青岛 266580；3.青岛新奥燃气有限公司，山东 青岛 266580）

研究管道内水合物颗粒的聚集过程对深水油气管道的流动安全保障具有重要意义［1-2］。目前，对管道内水合物颗粒聚集过程的实验多局限在静态条件下［3-4］，而研究流动条件下水合物颗粒的聚集过程则因为缺乏有效检测手段而大大受限。因此，随着近年来计算流体力学的发展，数 值模拟已成为研究水合物颗粒聚集过程的有效手段。

李莹玉等［5］利用弹性碰撞理 论描述了水合物颗粒间的聚集，进而模拟了流动过程中水合物颗粒的粒径变化，然而该模拟只考虑了对心碰撞且忽略了由碰撞导致的颗粒破碎。引入水合物颗粒聚集粒径计算模型［6-7］进行水合物颗粒聚集模拟是国内外学者普遍采用的一种方法［8-9］，但该方法取水合物颗粒间的聚集力为常数，且并未考虑水合物颗粒间的聚集效率。Balakin等［10-11］建立了水合物颗粒数量密度的群体平衡模型，并通过求解该模型模拟了管道内水合物颗粒的粒径变化。该模拟方法由于较为全面地考虑了水合物颗粒在流动过程中的聚集和破碎，因此模拟结果与工程实际更为接近。但他们在模型建立过程中进行了简化，假设了发生聚集的两个水合物颗粒粒径是相同的，因此该模型还有待进一步完善。

针对以上数值模型模拟的不足和缺陷，本工作引入了基于水合物颗粒聚集动力学的群体平 衡模型，通过碰撞频率和聚集效率描述水合物颗粒间的聚集，通过破碎频率和破碎后子颗粒的粒径分布函数描述水合物颗粒的破碎。再利用FLUENT 14.5 软件，对固液两相流模型和群体平衡模型进行联合求解，实现对管道内水合物颗粒聚集过程的数值模拟。

1 数值模型

1.1 几何模型

以 Balakin等［12］开展的管内 CCl3F（R11）水合物浆流动特性实验为模拟对象，将模拟结果与实验结果进行对比。根据实验所用环道建立了三维水平圆管模型（见图1），管道长度为3 m，管道内径为45.2 mm。对几何模型进行六面体网格划分，进口壁面以1 mm为步长进行划分；在近壁面处，为处理边界层效应，采用8层网格加密，其余网格均以1 mm为步长进行划分。图1几何模型总计267 696个六面体网格，网格质量0.913。对网格数量进行独立性检验（见图2）。从图2可看出，当网格数分别为527 616个和267 696个时，管流（2 m/s）充分发展段近壁面处速度梯度的相对偏差较小（6.3%），可见网格数量满足独立性要求。

1.2 多相流模型及湍流模型

在建模过程中采用以下假设：流动过程等温且忽略相间质量传递，即不考虑水合物的生成及分解；水合物浆由水相和水合物颗粒相构成，不考虑气相和油相；水合物颗粒粒径连续分布；不考虑水合物颗粒在管壁上的黏附过程。在上述假设下，采用的物理模型主要包括多相流模型、湍流模型和群体平衡模型。其中，多相流模型采用欧拉-欧拉双流体模型，由连续性方程（见式（1））、动量方程（见式（2））控制方程和若干用于封闭方程组的本构方程构成。

式中，i为水相或水合物颗粒相；t为时间，s；ρ为密度，kg/m3；α为体积分数，%；∇为拉普拉斯算子；u为速度矢量，m/s；p为压力，Pa；τi为应力张量，Pa；Mi则为相间动量交换项，kg/（m·s）2。

水合物浆的液固耦合是多相流模型建立时考虑的重点，在模拟过程中，液固耦合通过相间动量交换实现。在计算Mi时，主要考虑相间 拖曳力（Mdi）和湍流扩散力（Mti），并以此建立相间作用力模型以封闭多相流模型，见式（3）。

式中，ur为相间相对流速，m/s；μtm为湍流黏度，kg/（m·s）；σd为Planck扩散系数；αs为水合物体积分数，%；kls为动量传递系数，本工作采用Gidaspow模型［13］计算：当αs≤20%时，由Wen-Yu公式（见式（4））计算；当αs＞20%时，由Ergun公式（见式（5））计算。

式中，L为水合物颗粒粒径，m；μl为水相动力黏度，kg/（m·s）；CD为曳力系数。

除了相间作用力模型，多相流模型的封闭还需确定水合物颗粒相的黏度μs，计算式见式（6）［14］：

式中，αl为水相体积分数，%；μm为混合相黏度，kg/（m·s），可用 Roscoe-Brinkman 方程［15-16］（见式（7））计算得到。

根据式（6）～（7）即可编制用户定义函数（UDF）并以此计算水合物颗粒相的黏度μs。湍流模型则采用标准k - ε模型，近壁面采用标准壁面函数进行处理。

1.3 群体平衡模型

根据本工作假设，不考虑水合物颗粒的生成分解及水合物颗粒的管壁黏附且认为水合物颗粒粒径连续分布，因此可采用式（8）的群体平衡方程［17］：

式中，n（L，t）表示粒径为L的水合物颗粒在t时刻的数量密度，m-3；β（L-L′，L′）表示粒径分别为L-L′和L′两水合物颗粒的碰撞频率，m3/s；a为两水合物颗粒发生碰撞后的聚并效率；S（L′）表示粒径为L′的水合物颗粒的破碎频率，s-1；b（L|L′）则表示粒径为L′的水合物颗粒破碎后产生粒径为L水合物颗粒的概率。

基于水合物颗粒聚集动力学，对式（8）中关键参数的计算公式进行选取。

碰撞频率方面，主要考虑由差速沉降和流动剪切造成的碰撞并取两者碰撞频率之和作为水合物颗粒的实际碰撞频率。其中，差速沉降碰撞频率（βDS）采用式（9）［18］计算：

式中，V为沉降速率，m/s，可由式（10）计算得到：

对于流动剪切碰撞频率，当水合物颗粒小于Kolmogorov尺度时，用式（11）［19］进行计算；当水合物颗粒大于Kolmogorov尺度时，用式（12）［20］进行计算：

式中，G为流场局部的绝对速度梯度，s-1；u为水合物颗粒的平均速度，m/s。

采用曲线模型计算水合物颗粒间的聚并效率，见式（13）［21］：

式中，H为表征范德华力大小的哈梅克常数，J；R为发生碰撞两水合物颗粒的调和半径，m。

计算破碎频率时则主要考虑由流动剪切导致的破碎，计算式见式（14）：

式中，E和m为经验常数。

采用二元分布作为水合物颗粒破碎后的粒径分布函数［22-23］，表达式见式（15）：

根据式（9）～（15）编制UDF并以此计算群体平衡模型中的关键参数。

1.4 模型求解

利用FLUENT 14.5软件，对固液两相流模型和群体平衡模型联合求解。管道入口边界条件设定为速度入口，管道出口边界条件设定为压力出口，以二阶迎风差分格式离散动量方程，压力-速度耦合使用SIMPLE算法，群体平衡方程使用离散法进行求解，瞬态模拟时间步长设置为0.1 s，残差设置为10-5。模型求解时所涉及的部分参数见表1。

表1 模型参数Table 1 Model parameters

2 结果与讨论

2.1 模拟工况

本工作主要模拟不同条件下水合物颗粒的聚集过程，并分析管内流 速和水合物体积分数对水合物颗粒聚集过程的影响。模拟时的流速和水合物体积分数取值均参考文献［12］，具体见表2。

表2 模拟工况Table 2 Simulated working conditions

在使用群体平衡模型时，需给出水合物颗粒的初始粒径分布。为保证可对比性，模拟的7组工况 的水合物颗粒初始粒径均相同，即水合物初始粒径均为5 μm左右［25-26］（不可破碎），初始粒径分布见表3。

表3 水合物颗粒 初始粒径分布Table 3 Initial particle size distribution of hydrate particles

2.2 实验验证

依据相关实验［12］，取流动压降和水合物颗粒粒径分布对数值模型的可行性和有效性进行验证。表4为单位管长压降实验值与模型模拟值的对比，由表4可知，工况5～7单位管长压降实验值与模拟值的变化趋势相同且两者间的相对误差小于18%。因此，建立的数值模型能较好地模拟水合物颗粒在管道内流动时的压降变化规律。

工况1～5水合物颗粒粒径分布实验值与模拟值的对比见图3。由图3可知，由数值模型模拟得到的颗粒粒径分布与通过实验获得的颗粒粒径分布规律相同，均近似呈对数正态分布。因此，建立的数值模型也可较好地模拟水合物颗粒在管内流动时的粒径分布情况。

表4 实验和模拟条件下单位管长压降对比Table 4 Comparison of pressure gradient between experimental results and simulation results

图3 水合物颗粒粒径分布实验值与模拟值对比Fig.3 Comparison of experimental value and simulation value of hydrate particle size distribution.

2.3 低流速下水合物颗粒聚集过程分析

采用工况4～7分析低流速下水合物颗粒的聚集过程。图4为较低流速（1.5 m/s）下，使用数值模型模拟和以浓度差为推动力的水合物颗粒生长模型分别计算得到的管道内水合物颗粒平均粒径的变化情况。由图4可知，在较低流速下，由数值模型模拟得到的水合物颗粒平均粒径变化过程可大致分以下4个阶段：1）缓慢增长阶段。该阶段水合物颗粒粒径较小，碰撞频率较低，因此平均粒径增长缓慢。在该阶段，随管道内水合物体积分数的升高，水合物颗粒平均粒径的增长速度逐渐增大。这是因为高体积分数下水合物颗粒间更容易发生碰撞和聚集。2）快速增长阶段。该阶段水合物颗粒粒径由小变大，碰撞频率升高，水合物颗粒平均粒径开始快速增长。在该阶段，管道内水合物的体积分数越高，水合物颗粒平均粒径开始快速增长的时间就越早，同时水合物颗粒平均粒径的增长速度和所能达到的最大粒径也越大。3）破碎主导阶段。在阶段1和阶段2，颗粒聚集占据了主导地位，因此水合物颗粒平均粒径呈增长趋势。当水合物颗粒粒径增大到一定程度时，破碎过程将超过聚集过程占据主导地位，此后水合物颗粒进入破碎主导阶段，平均粒径开始逐渐降低。值得注意的是，工况4由于水合物浓度较低、颗粒粒径较小，故在平均粒径变化过程中水合物颗粒的聚集过程和破碎过程能逐渐达到平衡，因此并未出现明显的破碎主导阶段。4）动态平衡阶段。当水合物颗粒粒径降至一定程度时，聚集过程重新占据主导地位，水合物颗粒平均粒径开始逐渐增大，直到管道内水合物颗粒的破碎过程和聚集过程达到动态平衡。此后，管道内水合物颗粒的平均粒径逐渐维持稳定，不再发生明显变化。在该阶段，管道内水合物体积分数越高，水合物颗粒所能维持的最终稳定粒径就越大。

图4 低流速下由数值模型模拟和生长模型分别计算得到的管道内水合物颗粒平均粒径变化情况Fig.4 Variations of average particle sizes of hydrate particles in pipeline calculated by numerical model simulation and growth model at low flow rates.

图4中使用的以浓度差为推动力的水合物颗粒生长模型见式（16）［11］：

式中，r为水合物颗粒粒径的生长速率，m/s；kS为水合物颗粒的生长速率常数，m4/(s·mol)；cb为液相主体中R11的浓度，mol/m3；ceq为水合物颗粒与液相主体交界面处R11的浓度，mol/m3，当水合物颗粒生长过程结束时ceq=cb。假设R11可全部溶解在水中且水合物生长过程中R11呈指数规律消

耗，则cb和ceq可分别通过式（17～18）进行计算：

式中，cin为液相主体中R11的初始浓度，mol/m3；tmax为水合物颗粒生长过程结束所需的时间，s。

在使用上述生长模型时，本工作假设水合物颗粒的生长驱动力始终处于最大状态，即始终有cb=cin。模型计算过程中，kS取1.3×10-9m4/(s·mol)［11］，cin则可根据文献［12］报道的实验数据计算得到。从图4可看出，即使在保持最大生长驱动力的情况下，各工况内由水合物颗粒聚集导致的粒径增长仍较由水合物颗粒生长导致的粒径增长更为快速、剧烈。由水合物颗粒生长带来的粒径增长速率仅相当于聚集过程中缓慢增长阶段的粒径增长速率。这说明，相对于颗粒生长，水合物颗粒间的聚集更容易导致颗粒粒径增大、是流动过程中水合物颗粒粒径增大的主要原因。

2.4 高流速下水合物颗粒聚集过程分析

利用工况1～3分析高流速下水合物颗粒的聚集过程。较高流速（2.4～3.9 m/s）下，通过数值模型模拟得到的管道内水合物颗粒平均粒径的变化情况见图5。由图5可知，不同于低流速（1.5 m/s）下的工况，高流速下管道内水合物颗粒的平均粒径变化过程只包含快速增长（阶段Ⅰ）和动态平衡（阶段Ⅱ）阶段。

图5 高流速下由数值模型模拟计算得到的管道内水合物颗粒平均粒径变化情况Fig.5 Variations of average particle size of hydrate particles in pipeline calculated by numerical model simulation at high flow rates.

根据水合物颗粒聚集动力学可知，流动剪切是造成水合物颗粒发生碰撞聚集的主要原因。管道内流体流速越高，流动剪切就越强，因此高流速下水合物颗粒往往能更快、更剧烈地发生碰撞聚集。故高流速工况下水合物颗粒平均粒径增长较快，粒径变化过程中并未出现明显的缓慢增长阶段。另一方面，同工况4类似，工况1～3与工况5～7相比，水合物浓度较低、颗粒粒径较小，故在粒径变化过程中水合物颗粒的聚集过程和破碎过程能逐渐达到平衡，因此并未出现明显的破碎主导阶段。此外，从图5还可看出，管道内流体流速越高，快速增长阶段水合物颗粒平均粒径的增长速度和其所能达到的最大粒径值就越大，这同样是由高流速下水合物更容易发生碰撞聚集所致的。

3 结论

1）建立的数值模型能较好地模拟水合物颗粒在管道内流动时的压降变化规律以及粒径分布情况。

2）低流速下管道内水合物颗粒的聚集过程可大致分为缓慢增长、快速增长、破碎主导和动态平衡四个阶段。高流速下管道内水合物颗粒的聚集过程则只包括快速增长和动态平衡两个阶段。

3）在相同流速下，管道内水合物体积分数越高，水合物颗粒的平均粒径就越大。在相同水合物体积分数下，管道内水合物颗粒的平均粒径随流速的升高而增大。

4）在管内流动条件下，相对于颗粒生长，水合物颗粒间的聚集更容易导致颗粒粒径增大、是流动过程中水合物颗粒粒径增大的主要原因。

［1］ Sloan E D. Natural Gas Hydrates in Flow Assurance［M］. New York：Elsevier Science Ltd，2010：1-36.

［2］ Song Guangchun，Li Yuxing，Wang Wuchang，et al. Investigation of hydrate plugging in natural gas+diesel oil+water systems using a high-pressure fl ow loop［J］.Chem Eng Sci，2017，158：480-489.

［3］ Dieker L E，Aman Z M，George N C，et al. Micromechanical adhesion force measurements between hydrate particles in hydrocarbon oils and their modifications［J］.Energy Fuels，2009，23（12）：5966-5971.

［4］ Morrissy S A，Lim V W，May E F，et al. Micromechanical cohesive force measurements between precipitated asphaltene solids and cyclopentane hydrates［J］.Energy Fuels，2015，29（10）：6277-6285.

［5］ 李莹玉，江国业，陈世一. 基于CFD的集输管道内水合物聚集行为仿真分析［J］.计算机应用与软件，2013，30（7）：101-103.

［6］ Muhle K. Flock Stability in Laminar and Turbulent Flow. In：Coagulation and Flocculation：Theory and Applications［M］.New York：Marcel Dekker，1996.

［7］ Camargo R，Palermo T. Rheological properties of hydrate suspensions in an asphaltenic crude oil［C］// International Conference on Gas Hydrates，Yokohama：ICGH4，2002.

［8］ 陈鹏，刘福旺，李玉星，等. 水合物浆液流动特性数值模拟［J］.油气储运，2014，33（2）：160-164.

［9］ Fatnes E D. Numerical simulations of the flow and plugging behavior of hydrate particles［D］.Bergen：University of Bergen，2010.

［10］ Balakin B V，Hoffmann A C，Kosinski P. Experimental study and computational fluid dynamics modeling of deposition of hydrate particles in a pipeline with turbulent water flow［J］.Chem Eng Sci，2011，66（4）：755-765.

［11］ Balakin B V，Hoffmann A C，Kosinski P. Population balance model for nucleation，growth，aggregation，and breakage of hydrate particles in turbulent fl ow［J］.AIChE J，2010，56（8）：2052-2062.

［12］ Balakin B V，Pedersen H，Kilinc Z，et al. Turbulent fl ow of freon R11 hydrate slurry［J］.J Pet Sci Eng，2010，70（3/4）：177-182.

［13］ Ding Jianmin，Gidaspow D. A bubbling fluidization model using kinetic theory of granular fl ow［J］.AIChE J，1990，36（4）：523-538.

［14］ 王继红. 冰浆的管道输送热流动特性［D］.大连：大连理工大学，2013.

［15］ Pabst W. Fundamental considerations on suspension rheology［J］.P R Soc A，2004，48（1）：6-13.

［16］ 赵鹏飞，王武昌，李玉星，等. 管道内水合物浆流动的数值模型［J］.油气储运，2016，35（3）：272-277.

［17］ Hulburt H M，Katz S. Some problems in particle technology：A statistical mechanical formulation［J］.Chem Eng Sci，1964，19（8）：555-574.

［18］ Camp T R，Stein P C. Velocity gradients and internal work in fl uid motion［J］.J Bsn Soc Civ Eng，1943，30（4）：219-237.

［19］ Saffman P G，Turner J S. On the collision of drops in turbulent clouds［J］.J Fluid Mech，1956，1：16-30.

［20］ Abrahamson J. Collision rates of small particles in a vigorously turbulent fluid［J］.Chem Eng Sci，1975，30（11）：1371-1379.

［21］ Ven T G M V D，Mason S G. The microrheology of colloidal dispersions Ⅶ . Orthokinetic doublet formation of spheres［J］.Colloid Polym Sci，1977，255（5）：468-479.

［22］ Zhang Jianjun，Li Xiaoyin. Modeling particle-size distribution dynamics in a flocculation system［J］.AIChE J，2003，49（7）：1870-1882.

［23］ 李振亮. 基于群体平衡的活性污泥絮凝动力学［D］.重庆：重庆大学，2014.

［24］ Li Xiaoyin，Logan B E. Collision frequencies between fractal aggregates and small particles in a turbulently sheared fl uid［J］.Environ Sci Technol，1997，31（4）：1237-1242.

［25］ 吕晓方，吴海浩，史博会，等. 流动体系中二氧化碳水合物堵管时间实验研究［J］.实验室研究与探索，2013，32（11）：197-202.

［26］ Turner D J，Miller K T，Sloan E D. Direct conversion of water droplets to methane hydrate in crude oil［J］.Chem Eng Sci，2009，64（23）：5066-5072.