油水分离器中液滴流动模拟方法研究

在石油行业普遍存在着油水两相流动，具体到油气储运领域也不例外［1］.我国陆上注水开采油田大多进入开发后期，采出液含水量高，油水两相流并不鲜见于油田集输管道和储存设施中.与此同时，在远洋深水等地域大力拓展的石油开采也由于成本原因很难在井口进行油水分离，也导致了油水两相流广泛存在于长距离集输管道中.

然而，储运设施中油水两相流的存在会带来如下影响：①由于油水密度的差异，会导致水相沉积形成分离层润湿管道底部，造成管道内化学和电化学腐蚀加剧，从而导致管道穿孔，石油泄漏扩散，进而造成土壤和水环境污染，并带来火灾、爆炸等灾害风险［3］.②原油在收集、矿场加工、输送和储存过程中，不时需要加热升温，而水相的存在，由于其比热容较高（通常为油相的两倍），增大了燃料消耗，占用了集油、加热资源，增加了原油生产和储运成本.综上所述，必须对油品在储运过程中进行油水分离处理.

1 实验模拟方法分析

要想“百战不殆”地进行原油脱水，就需要对脱水原理和油水两相流动过程“知己知彼”.虽然对多相流通过实验进行流型分析的研究多年来一直在进行，但不可否认的是以流型为依据对多相流进行分析研究的方法一直存在着流型定义不明确、不统一和不定量的问题.具体来说有以下三点：①受观察条件影响，目前常用的实验用油都是十六烷烃或其他长链烷烃，大多透明且轻，这与我们常见的油品种类有很大差别，尤其是油水黏度的巨大差异会对液滴运动和融合过程带来显著影响；②部分实验采用的光源不是严格的冷光源，例如激光等光源条件更会对尺寸极小的液滴带去足以影响其原本物性的能量影响，光所产生的波动也会起到对表面的扰动作用；③实验室出于摄像对焦需要，所选用的流道尺寸往往只有毫米级，这就会导致此时流道的边界效应作用被急剧放大，影响到了实验的真实性.

此外，工业计算仍主要以基于实验的、按照流型区分的经验和半经验关联式为数学模型，所导致结果是：①对实验结果背后的非平衡动力学、化学等复杂的热物理过程不能准确理解，对相关机理无法进行深入分析研究.而对这些过程和相关机理的完整理解将为具体的工业应用提供更坚实的理论基础和极有价值的参考信息.②经验公式往往并不具有普遍性，解决具体问题时通常需要先做实验确定某些关键系数的值，而直接采用经验值又会导致计算结果偏离实际情况.

填补这两个缺憾的最佳途径为直接数值模拟方法，即不需要建立依据流型的计算模型，而是直接通过求解控制方程来实现对流动的数值模拟.此外，由于在流动中存在搅拌、剪切等作用，水相大多以液滴的形式分散在油相中.因此，我们从直接求解液滴的控制方程入手，分析水滴在油流中的位置、变形和融合过程，弄清油水两相流动的基本规律，从而为储运设备的设计以及安全、高效、可靠运行提供理论参考和决策依据.

2 数学模型构建

用宏观计算流体力学方法模拟多相流动过程，通常的做法是在经典的Stokes方程的基础上，引入描述相分布的控制方程，联立求解得到两相流体的运动参数变化规律.

2.1 Cahn-Hilliard方程简介

Cahn-Hilliard方程的一般形式为：

式中：M为迁移率，有时是定值，有时是ϕ的函数；而ϕ称为序参量，其定义是：

显然，ϕ具有如下性质：

式中：ρ1，ρ2为1和2两种流体的当地密度.显然，ϕ=0的位置应当在界面处，这也是我们检测模拟结果是否合理的一个判据.

f()ϕ表示单位体积内的自由能，可以用下式表示：

系统的总能量为

式中，ε为界面参数，而积分内第一项表示的是界面能.

2.2 Cahn-Hilliard方程离散

化学势作为两相流动中引起相间分离和流体相扩散的主要驱动力，是对各相运动规律研究的基础，其用μ表示，

序参量ϕ对时间的偏导为：

有些文献中定义了J，

那么由式（7）可以得到：

序参量ϕ对空间的偏导为：

即：

那么可以得到：

积分得：

所以，序参量可以写成：

其中，tanh函数的表达式为：

进行数值模拟时，我们首先把单位体积内的自由能f()

ϕ写成如下形式：

式中：κ,r,u可以由下式得到：

那么，化学势可以由下式得到

2.3 Cahn-Hilliar-Stokes联立方程组的建立

在一个完整的数值模拟中，相分布的控制方程需要和动量方程以及质量方程相联立，以Cahn-Hilliard-Stokes联立方程为例，控制方程组应该写为：

式中：fe为系统所受外力，而σ=η(∇v+(∇v)T)，表示牛顿流体所受的黏性力.

边界条件的定义方面，除了对Stokes方程常见的边界条件限制外，还需要加上针对C-H方程的边界条件.具体可分为：

无通量的边界条件为：

即

自然边界条件为：

初始条件为：

2.4 Shift-Matrix算法

Shift-Matrix方法是孙树瑜等［3］于2012年提出的一种针对Stokes方程、基于编程语言简洁化的目的、在周期性边界条件限制下，尤其适用于Matlab语言编程使用的离散方法.章涛等［4］在2015年验证了其对Navier-Stokes方程及非周期性边界条件的适用性，并且用数值实验证明了其对偏微分方程的求解具有较快的速度.由于上一节所建立的联立方程组中含有多个偏微分项，因此将该算法引入模型的求解具有显著作用.

该方法定义了如下的偏移算子：

那么，x方向动量方程就可以写为：

为了计算简便，上式等号左边第一项乘

下一步，我们定义如下算子：

那么，x方向上的动量方程离散形式就可以写为：

同理，y方向上的动量方程为：

连续性方程可以写为：

那么，整个数学模型可以写成如下形式：

3 模拟结果与分析

3.1 水滴的融合过程

在一个静止的充满了油的管道中，有两个圆形水滴相切，且圆心连线过切点.两水滴大小相同，半径均为12.红色区域在此处代表水；蓝色区域在此处代表油.此外，控制相界面厚度的β值取为0.8.油水密度比为0.8.所得的初始相分布可用图1表示.为了证实图中颜色代表流体种类的性质，点击深蓝色的点，可以读出此时的序参量为-1，代表此处为纯油，符合初始条件设定.

从图1中可以看出，由于水滴相接，在表面张力作用下，两水滴开始融合.连接面积不断扩大，从交接点处向东西两侧扩展，直到两水滴整体呈“两头圆、中间方”的形状.此后，继续由表面张力作用，水滴团向东西方向持续扩展，南北方向开始明显压缩，整体向圆形发展.针对同一初始条件、边界条件和时间步长，用LBM方法中的颜色阶梯模型作为数值模型也能得到同样的结果，在此不做赘述，只比较其计算耗时不同（表1）.

从表中可以看出，收敛标准越严苛，颜色阶梯模型计算耗时的劣势越明显.从误差收敛图中可以发现，宏观尺度的S-C-H模型的数值模拟收敛较快也较稳定，而LBM方法中的颜色阶梯模型误差收敛较慢，且存在一定的波动现象.这也说明了引入Shift-Matrix算法的S-C-H模型更适合于在Matlab语言编程条件下对油水两相流进行模拟.

图1 模拟开始时刻相分布图及t=100，500，3000时相分布图
Fig.1 Phase distribution at the beginning and change during the merge process when t of 100，500，3000

前文已经对两个相同的水滴相接时的融合过程模拟，根据结果可以推测：如果三个水滴呈直线状相接，那么在无重力作用环境下的纯融合过程也会呈现相似的结果，即从相接处开始向东西两侧扩展，南北方向压缩直至水滴团整体呈圆形，如图2～图4所示.

表1 两种方法计算耗时比较
Tab.1 Comparison between simulation time of the two methods

图2 两种方法的误差收敛过程比较
Fig.2 Comparison between error reduction of the two methods

图3 三水滴融合的初始时刻相分布
Fig.3 Phase distribution of three water droplets at the beginning

3.2 水滴的沉降过程

以长为1、宽为1的静止油流作计算区域，网格划分成100×100，其中有一圆形水滴，其圆形初始时刻位于计算区域中心处.水滴半径为0.18，油水密度比为0.8，如图5所示.

由图6可知，模拟开始后，由于油水密度差的存在，受重力作用水滴开始下降，下降过程中伴随着形状的变化，最终水滴呈扁平状摊在区域底端.

图4 融合稳定后的相分布
Fig.4 Phase distribution after fusion in a stable state

图5 相分布初场
Fig.5 Phase distribution at the beginning

图6 t=300，3000时的相分布图
Fig.6 Change of phase distribution during the merge process when t of 300，3000

在油气储运行业，油水分离器的设计中常常需要计算水滴的下落速度.按照《油气集输与矿场加工》［14］中提供的方法，水滴的下落速度应由斯托克斯沉降公式确定：

式中：dl代表液滴的直径；ρw代表水的密度；ρo代表油的密度；μ代表动力黏度.按该式计算，本算例中的水滴下降速度应为

然而，对下降过程进行数值模拟后，水滴的降落速度变化如图7所示.

图7 水滴下降速度的变化
Fig.7 Change in dropping velocity of the water droplet

从图7中可以发现，水滴下降的速度并不是恒定于某一值，而是呈现出先升后降的起伏现象，并且囊括了公式（34）计算得到的结果.究其原因，一是由于在下降过程中水滴不断发生变形，受力面积改变，受力大小和方向也随之改变，公式（34）所必需的力平衡条件并不存在；二是由于水滴的形状不规则，其各部分的下降速度并不一样，图7所示的下降速度取的是水滴上各个计算节点的速度的平均值.此外值得一提的是，水滴底部降落到计算区域底部以后，在工程中可能认为此时水滴已经静止了（已经被当作脱离出分离器了）；然而实际情况下，水滴上部仍然受重力作用继续向下运动，只是因为底部被限制住因而只能发生变形.在模拟时，这部分水滴的变形速度会被记录并且平均到底部不动的那部分水滴上，导致计算得到的速度虽然在缓慢减小但始终存在.所以，图7并没有一直显示到水滴静止时的结果，因为此后水滴速度会保持相当长的缓慢减小时间.

3.3 水滴的融合沉降过程

为了进一步证明本文设计的数学模型和计算代码的正确性，下面将对三个油滴的运移和融合过程进行模拟.初始时刻三个水滴并列于一条直线上，大小相同，半径为0.15，油水密度比为0.8.为了更好地显示效果，计算区域大小改为1×2，网格数相应地改为100×200.从图8中可以发现，当有三个油滴同时存在于计算区域时，上面的油滴因为先到达区域顶端从而先开始变形，中间的油滴上升过程中形态也逐渐与单个油滴时的形态产生不同，下部的油滴因为上升时间最长因而最接近单个油滴存在时的情况.所以，多个液滴的运移融合过程，即受到液滴本身的初始位置、大小等因素影响，也会因相互间作用关系的不同而发生变化.

图8 三个油滴在平衡态水流中的上升过程
Fig.8 The rising process of three oil droplets in the equilibrum state of water flow

4 结论

本文针对油气储运领域广泛存在且必须处理的油水两相流动问题设计了数学模型，并进行了直接数值模拟研究.对宏观计算流体力学方法中的Stokes方程与Cahn-Hilliard方程联立模型进行了详细推导和离散，并且引入了偏微分方程的快速求解方法——Shift-Matrix方法，将其应用到了求解S-C-H方程组中.

在油流静止且重力影响忽略的前提下，验证了本文所给出的C-H-S模型的正确性.通过比较其与目前常用的LBM方法计算耗时的不同，得出了引入Shift-Matrix算法的C-H-S模型在Matlab编程中具有较快的计算速度；随后经过对误差收敛过程的比较，并结合Matlab语言的特点，分析了造成计算速度差异的原因.最后利用更适合于Matlab编程的C-H-S模型模拟了三个水滴融合的过程，证明了本方法的普适性.

在引入重力的影响后，本文对液滴的运移及融合问题开展了研究.研究发现，目前油气储运行业在设计重力分离器时常用的斯托克斯沉降公式并不能准确描述水滴的降落速度，而且水滴在油流中下降时还伴随有变形过程，这也是以前很少考虑到的.液滴的排布方式不同，其运移方向、速度及变形规律也不一样.这些结果在实际生产中都会对油水分离时的方法选择和剂量控制产生一定影响.

本文通过对不同条件下水滴和油滴的运移及融合过程的模拟，验证了本文设计的数值模型和编写的计算代码的可靠性和正确性，从而能够为日后油气储运领域实际油水两相流过程（更多的液滴，更多的条件组合）的研究提供理论依据和技术基础.