天然气是一种清洁、高效、丰富和廉价的主要能源,占全球一次能源消费的三分之一,加速天然气勘探开发,发展天然气产业,对保障我国能源安全、推动经济增长及改善环境具有重大意义。随着天然气在我国能源结构中地位的不断上升,天然气净化工业的地位也越发重要。针对高酸性组分天然气及高CO2含量二氧化碳驱油田伴生气的处理,迫切需要开发新的技术。
近年来,超声速旋流分离技术开始应用于天然气处理加工领域,国内外研究学者先后将其应用于天然气脱水、脱重烃和液化中,利用超声速流动条件下气体的低温凝结效应,结合旋流分离技术,实现天然气中某些组分的冷凝分离[1-15]。该技术结合了气体动力学、工程热力学和流体力学理论,将膨胀降温、旋流式气液分离、再压缩等处理过程集中在一个密闭紧凑的装置中完成,具有密闭无泄漏、无需化学药剂、结构紧凑轻巧、节能环保等优点,特别适用于高酸性组分天然气的净化处理。在这个过程中,酸性组分在超声速喷管内的凝结是实现其脱除的前提和关键。因此,以天然气中的主要酸性组分CO2为例,开展了超声速喷管内CO2气体凝结特性的研究,对CO2气体在喷管内的凝结流动过程进行数值模拟,分析凝结流动参数的分布规律。该研究对于今后进一步开展天然气酸性组分超声速旋流分离机理研究,完善超声速旋流分离技术,丰富和发展天然气净化处理工艺具有重要的指导意义。
采用Girshick等[16-17]推导的内部一致经典成核理论(ICCT)进行成核率的计算,见式(1)。
式中:S为气体过饱和度,$S{\rm{ = }}\frac{{{\rho _{\rm{v}}}}}{{{\rho _{\rm{l}}}}}$;rc为临界半径,${r_{\rm{c}}}{\rm{ = }}\frac{{2\sigma }}{{{\rho _{\rm{l}}}{R_{\rm{M}}}T\left( {{\rm{In}}S} \right)}}$,m;θ为无因次表面张力,$\theta = \frac{{\sigma {a_0}}}{{{k_{\rm{B}}}T}}$;a0为分子表面积,a0=(36π)1/3(vl)2/3,m2。
采用Gyarmathy液滴生长模型计算液滴生长速率[18],见式(2)。
式中:Kn为无量纲数,其定义为$Kn = \frac{{\bar l}}{{2{r_{\rm{d}}}}}$;l为蒸汽分子的平均自由程,$\bar l = \frac{{2\eta \sqrt {{R_{\rm{M}}}T} }}{{{p_{\rm{v}}}}}$,m。
基于欧拉——欧拉双流体模型,忽略气、液相间速度滑移,同时认为液滴温度为当地压力下的饱和温度,建立气相及液相流动控制方程组。气相包括质量、动量及能量方程,见式(3)~式(5);液相包括质量、液滴数目及液滴半径方程,见式(6)~式(8)。
式中:$m = J{\rho _{\rm{l}}}\frac{{4{\rm{ \mathsf{ π} }}r_{\rm{c}}^3}}{3} + N{\rho _{\rm{v}}}{\rho _{\rm{l}}}4{\rm{ \mathsf{ π} }}r_{\rm{d}}^2\frac{{{\rm{d}}{r_{\rm{d}}}}}{{{\rm{d}}t}}$,为单位时间单位体积内凝结的液滴质量,kg/(m3·s)。
以上方程中添加了由于液滴凝结而产生的源项:Sm=-m;Su=-mu;Sh=m(hlv-h);SY=m。
由于凝结成核模型中表面张力σ以立方的形式出现,其微小变化可能对成核率的计算带来重大影响,进而影响数值计算结果的准确度。因此,有必要选择更为准确的表面张力计算方法。由于微小液滴的表面张力难以测得,根据平面液体表面张力对比计算值和实验值,一般认为,基于对比态原理和等张比容的表面张力计算方法最为精确[19]。因此,本研究优先选用对比态关联式进行表面张力的计算,该关联式适用于非极性液体,见式(9)。
式中:$Q = 0.119{\rm{ }}6(1 + \frac{{{T_{{\rm{br}}}}{\rm{ln}}({p_{\rm{c}}}/1.013)}}{{1 - {T_{{\rm{br}}}}}}) - 0.279$;${T_{{\rm{br}}}} = \frac{{{T_{\rm{b}}}}}{{{T_{\rm{c}}}}}$,为气体在正常沸点下的对比温度。
UDS为用户自定义标量(User-Defined Scalar)的简称,Fluent软件允许用户自己定义一些标量来解决标准模块中不能解决的问题。本研究采用UDS输运方程添加液相流动控制方程组。
UDF为用户自定义函数(User-Defined Function)的简称,它是Fluent软件提供的一个用户接口,用户可通过UDF,以Fluent软件为平台进行二次开发和个性化设置。本研究利用C语言编写了相应的UDF程序,将源相表达式嵌入Fluent软件中进行计算。考虑真实气体效应,在UDF中对真实气体密度、比热容、黏度、导热系数、凝结潜热及液滴密度等参数进行计算,查阅相关文献选择合适方法。
Laval喷管中的气体流动属于高速可压缩流动,采用密度基方法求解。湍流模型选用k-ω方程。气、液相流动控制方程组、湍流动能方程、湍流耗散率方程均采用二阶迎风格式进行离散。进口设置为压力进口,指定总压、静压、总温及湍流参数。出口设置为压力出口,因超声速流动的所有流动属性均从内部推算得到,故不进行相应设置。固壁设置为无滑移、无渗流、绝热边界。
通过目前文献中已有的喷管内CO2气体自发凝结实验数据,验证前述数学模型及数值计算方法,并对喷管内的CO2气体凝结特性进行分析。喷管结构及实验数据选自文献[20]。由于CO2气体液相区范围较窄,实验时难以保证其完全处于液相区而没有固相出现,故模型中不对液固两相进行区分。
实验所用喷管总长132.08 mm,喉部半径1.35 mm,位于35.56 mm处。其结构及尺寸见图 1。采用非结构化网格对喷管进行网格划分,并进行无关性验证。网格总数确定为7 338,网格划分情况见图 2。
选取文献[20]中nozzle I Run 1/2/5/6共4组实验数据,按照前述模型和方法进行数值计算,将数值计算所得的喷管内压力分布与实验结果进行对比。4组实验中喷管的入口参数列于表 1。
喷管内CO2凝结数值计算压力分布与实验结果对比见图 3。由图 3可知,数值计算所得喷管内压力变化趋势与实验结果基本一致,在可接受的误差范围内,数值计算值略小于实验值,表明所采用的数学模型和数值计算方法可较准确地反映喷管内气体自发凝结流动过程。值得说明的是,数值计算和实验结果均未捕捉到明显的凝结冲波现象,即由于蒸汽凝结导致蒸汽流量减小,释放大量潜热,喷管内气体在凝结的瞬间出现压力和温度的突跃。经分析,喷管内气体的自发凝结与喷管结构、气体性质、入口参数及环境条件等多种因素有关。其中,气体性质、入口参数对凝结冲波有较大影响。在前述实验工况下,CO2气体凝结潜热较小,凝结时释放的热量少,对周围气体的加热作用不明显,因而未对周围气体的压力和温度产生显著影响。从实验数据来看,喷管内压力下降较平稳,没有出现突跃,与数值计算结果基本一致,也从侧面证明了上述分析的合理性。
以Run 1实验工况为例,对喷管内CO2气体自发凝结过程中凝结流动参数的分布进行分析。
图 4(a)~(d)为喷管内气体压力、马赫数、温度和过冷度的分布。CO2气体进入喷管后,随着喷管的收缩及扩张,气体马赫数不断增大,压力和温度逐渐降低,过冷度增加。特别是气体经过喉部之后,压力和温度迅速降低,过冷度急剧增加,最大可至约30 K,产生自发凝结现象。凝结发生后,过冷度快速下降,于喷管出口降至约5 K。
图 4(e)~(h)分别为成核率、液滴数目、液滴半径和湿度的分布。结合图 5气体凝结参数沿喷管轴线的分布可以看出,在距离喉部约2.21 mm处(Wilson点),成核过程急剧发生。在极短的时间和距离内,成核率由0迅速增至2.04×1021 m-3·s-1,气流中出现大量临界尺寸的凝结核心,液滴数目也随之由0激增至1015的数量级。凝结核心的出现和环境中较大的过冷度使得气体分子不断在液滴表面团聚、液化,促使液滴不断生长,液滴半径和湿度迅速增加。此后,由于凝结释放了潜热,过冷度迅速下降,成核率由峰值急剧减小为0,喷管内不再有新的凝结核心生成,成核过程结束。由于数值计算中未考虑液滴的碰撞、聚并和破碎,液滴数目基本保持不变。此后,由于过冷度仍能维持在5 K以上,气体分子在液滴表面的团聚与液化过程得以继续,已有的凝结核心不断增长,液滴半径和湿度继续增加,直至喷管出口。最终,液滴半径可增至1.46×10-7 m,出口湿度可达0.093 5。
本研究基于欧拉——欧拉双流体模型建立气相及液相流动控制方程组,结合凝结成核与液滴生长理论,在对CO2气体物性进行充分研究的基础上,选取合适的数值计算方法,对喷管内CO2气体凝结特性进行研究,分析了喷管内气体凝结流动参数的分布情况,主要结论如下:
(1) 数值计算所得喷管内压力变化趋势与文献中实验结果基本一致,在可接受的误差范围内。数值计算值略小于实验值,说明所采用的数学模型和数值计算方法可较准确地反映喷管内气体自发凝结流动过程。由于实验工况下CO2气体凝结潜热较小,凝结释放的热量对周围气体的加热作用不明显,因而没有表现出明显的凝结冲波现象。
(2) CO2进入喷管,特别是在经过喉部之后,气体马赫数不断增大,压力和温度逐渐降低,过冷度不断增加,最大可至30 K左右。凝结发生后,过冷度迅速下降,在喷管出口处降至5 K左右。
(3) 喷管内CO2气体的成核过程在时间和空间上表现出急剧性。凝结起始位置(Wilson点)位于距离喉部约2.21 mm处,成核率由0突跃至2.04×1021 m-3·s-1,液滴数目达到1015的数量级。此后,凝结释放的潜热使过冷度快速下降,喷管不再具备成核条件,液滴数目基本保持不变。
(4) 凝结核心形成后,气体分子在一定的过冷度下不断在液滴表面团聚、液化,液滴半径和湿度迅速增加。成核过程结束后,由于气体过冷度仍能维持在5 K以上,已有的凝结核心继续生长,最终液滴半径可增至1.46×10-7 m,出口湿度可达0.093 5。
按照文中出现的先后顺序,特对符号及所采用的计算公式说明如下:
J,液滴成核率,m-3·s-1;
S,气体过饱和度;
rc,临界半径,m;
θ,无因次表面张力;
a0,分子表面积,m2;
pv,气体压力,Pa;
pS,气体饱和压力,Pa,CO2饱和蒸汽压拟合公式[21];
ρv,气相密度,kg/m3,SRK方程;
ρl,液相密度,kg/m3,Rackett方程[19];
σ,液滴表面张力,N/m,对比态关联式;
m0,单个气体分子质量,kg;
kB,Boltzmann常数,取为1.380 6×10-23 J/K;
RM,气体常数,J/(kg·K);
T,气体温度,K;
vl,单个液滴分子体积,m3;
$\frac{{{\rm{d}}{r_{\rm{d}}}}}{{{\rm{d}}t}}$,液滴生长速率,m/s;
Kn,无量纲Knudsen数;
l,蒸汽分子的平均自由程,m;
λv,气体导热系数,W/(m·K),改进Eucken关联式、Stiel-Thodos法[19];
hlv,气体凝结潜热,J/kg,Watson公式[19];
TS,饱和温度,K,NIST数据拟合;
rd,液滴半径,m;
Prv,气体Prandtl数;
γ,气体比热比,BWRS方程[22-23];
η,气体动力黏度,Pa·s,Neufeld方程、剩余黏度关联法[19];
t,时间,s;
xi,轴向位置,m;
xj,径向位置,m;
ui,轴向速度,m/s;
uj,径向速度,m/s;
Sm,质量源项,kg/(m3·s);
Su,动量源项,kg/(m2·s2);
Sh,能量源项,J/(m3·s);
SY,湿度源项,kg/(m3·s);
m,单位时间单位体积内凝结的液滴质量,kg/(m3·s);
δij,Kronecker delta数;
ui',轴向速度波动,m/s;
uj',径向速度波动,m/s;
E,总能,J/kg;
keff,有效导热系数,W/(m·K);
τeff,有效应力张量;
h,气体总焓,J/kg;
ρ,混合相密度,kg/m3;
Y,湿度;
N,液滴数目,kg-1;
pc,临界压力,bar;
Tc,临界温度,K;
Tr,对比温度;
Tb,正常沸点,K。