天然气中硫化氢超声速凝结特性

曹学文，王春智，孙文娟,3，牟林升

(1.中国石油大学(华东)储运与建筑工程学院，山东 青岛 266580；2.滨州学院 化工与安全学院，山东 滨州 256600；3.中海油研究总院有限责任公司 规划研究院，北京 100028)

能源结构低碳化是世界能源发展的趋势，与石油和煤炭相比，天然气更加清洁低碳，中国《能源发展“十三五”规划》提出，到 2020 年，天然气占一次能源消费量的比重力争达到10%[1]，能源生产和消费革命将进一步激发天然气需求。含硫化氢天然气是天然气资源的重要组成部分[2]，目前我国已探明天然气储量中有相当部分为含硫天然气，且有些气田中硫化氢含量较高。根据天然气气质特点选择高效的脱硫工艺，是高含硫气田开发需要解决的重要问题。目前，胺法仍是应用最广泛的酸性天然气净化方法，但其存在设备系统复杂、成本和运行费用高等缺点，尤其是硫化氢含量较高时，脱硫溶剂循环量大、能耗高[3-5]。为降低含硫天然气开发成本，应大力研究发展能耗低、污染小的新型天然气脱硫技术。

超声速旋流分离是一种新兴的混合气体分离技术，最初主要应用于空调空气中水分的分离，后来荷兰Shell石油公司和俄罗斯ENGO石油公司将其引入天然气加工处理领域[6-7]。超声速旋流分离装置集膨胀降温、旋流式气-液分离、再压缩等过程于一体，具有结构简单可靠、无转动部件、无需化学药剂、支持无人值守等优点[8]，国内外学者针对其在天然气脱水及重烃方面的应用开展了大量理论及实验研究[9-10]，在天然气液化[11-12]、天然气脱碳[13-14]方面也取得了一定的研究进展，若能将其用于天然气中硫化氢的脱除，对于完善和发展天然气净化工艺、降低含硫天然气开发成本意义重大。

超声速旋流分离器主要由超声速喷管、旋流装置、扩压段等组成，酸性天然气进入超声速旋流分离装置后，在Laval喷管中高速膨胀至超声速，温度降低，当达到一定过饱和状态时，天然气中的硫化氢组分将发生凝结，形成液滴，同时在旋流场作用下实现气-液分离。硫化氢气体在喷管内的自发凝结是天然气超声速旋流分离脱硫技术的前提条件。目前关于高速膨胀过程中气体自发凝结特性已开展较多理论、模拟及实验研究[15-18]，对数学模型的建立、求解及实验方法等均有探讨，对于本文研究的开展具有较好的借鉴意义，但超声速气流凝结特性研究多以湿蒸汽、醇类、烃类气体为介质，目前尚缺乏对于硫化氢气体凝结特性的研究。因此，笔者结合凝结相变模型，利用计算流体力学软件FLUENT对天然气中硫化氢气体在超声速喷管内的凝结相变过程进行数值模拟研究，并分析入口压力、温度及背压对凝结流动过程的影响。

1 喷管设计

基于流动特征，Laval喷管可划分为亚声速收缩段、达到声速的喉部和超声速扩张段[19]。计算喉部尺寸是喷管设计的首要任务，考虑真实气体效应，选用BWRS方程进行喉部临界参数的计算，进而确定喉部截面积[20]；采用双三次曲线法设计喷管收缩段;采用圆弧加直线法设计喷管扩张段。同时考虑喷管边界层黏性修正，认为边界层厚度从喉部开始沿轴向线性发展，修正角取为0.5°[19]。

选取天然气主要组分甲烷和硫化氢构成二元物系，硫化氢的摩尔分数设为10%。根据上述设计方法，利用MATLAB软件编制喷管结构设计程序，针对入口压力5 MPa、入口温度273 K、入口流量5000 m3/h的工况，对喷管尺寸进行计算，设计所得喷管结构如图1所示。喉部半径为5.65 mm，位于x=128.27 mm处。

图1 Laval喷管结构图Fig.1 Diagram of Laval nozzle

2 自发凝结数值模型

2.1 流动控制方程

基于欧拉-欧拉模型，建立气、液相流动控制方程组。气相流动控制方程包括连续性方程、硫化氢气体质量守恒方程、动量守恒方程及能量守恒方程，如公式(1)～(4)所示，诸方程中添加由于硫化氢气体自发凝结而产生的源项。

(1)

(2)

(3)

(4)

由于凝结所形成的液滴粒径较小(亚微米级)，故忽略气、液相间速度滑移，认为气、液两相具有相同的流动速度。液相控制方程包括液滴数目及液相质量守恒方程，如式(5)～(6)所示，通过用户自定义标量(User-defined scalar,UDS)输运方程添加。

(5)

(6)

流动控制方程中的源相表达式如式(7)～(10)所示，公式中m按公式(11)计算。源相利用C语言编写用户自定义函数(User-defined function,UDF)程序嵌入到FLUENT中，计算凝结相变对质量、动量以及能量方程的影响。

Sm=-m

(7)

Su=-mu

(8)

Sh=m(hlv-h)

(9)

SY=m

(10)

(11)

2.2 凝结模型

成核是气体达到一定过饱和度后继而自发产生的结果。成核率定义为单位时间、单位体积过饱和气体中自发产生的临界液滴数目。由于凝结相变机理复杂，众多学者开展了大量研究，提出了许多理论和半经验模型。目前，Girshick等提出的内洽经典成核理论(Internally consistent classical theory，ICCT)[21-22]因其相对较高的准确性得到了较为广泛的应用。Lamanna[23]在实验研究的基础上，提出利用经验系数对ICCT进行修正，修正公式预测所得成核率数值与实验更为相符。因此，笔者利用Lamanna修正后的ICCT模型计算成核率，如公式(12)所示。

(12)

凝结核心形成后，更多气体分子将在液滴表面凝结，促使凝结液滴继续生长。考虑液滴与周围气体的传热、传质，采用Gyarmathy模型[24]计算液滴生长速率，如公式(13)所示。

(13)

2.3 数值计算方法

针对喷管中气体高速可压缩流动特征，选用基于密度法的求解器。雷诺应力项的计算采用k-ω湍流模型，该模型适用于墙壁束缚流动计算，且包含了可压缩性的影响。采用二阶迎风格式离散流动控制方程、湍流动能方程和湍流耗散率方程。

喷管入口采用压力入口，指定总压、静压、总温、气体组成及湍流参数；喷管出口采用压力出口，指定静压、回流总温、湍流参数；壁面采用无滑移、无渗流、绝热边界。

喷管内气流为甲烷和硫化氢混合气体，气体状态方程选用石油与天然气工程领域常用的Peng-Robinson (PR)[25]方程，气体混合物黏度、导热系数等热力学参数的计算[26]采用C语言编写UDF加入到FLUENT软件中。

模拟过程中，同时监测残差变量及流量的变化情况，当残差小于10-3，且入口和出口的质量流量(气、液相质量流量之和)相对误差小于0.05%时，才认为计算收敛。

2.4 网格无关性验证

采用CFD专业软件ANSYS ICEM CFD进行网格划分，利用Laval喷管结构的对称性，计算模型采用二维，考虑到喷管结构的不规则性，采用非结构化网格。为排除网格划分情况对数值模拟结果的影响，逐步加密网格，对网格进行无关性验证。网格数为2408、9216、13938及18541情况下喷管中心轴线处液相质量分数分布如图2所示(入口压力5 MPa，入口温度273 K)。从图2可以看出，当网格数大于13938时，计算结果趋于稳定。因此，综合考虑计算精度要求及计算机求解效率，在对喷管内酸性天然气凝结流动特性的数值计算中，网格数量选定为13938，网格划分情况如图3所示。

图2 不同网格数下液相质量分数(Y)沿喷管中心线分布Fig.2 Distributions of liquid-phase mass fraction (Y) along nozzle axis with different mesh cells Tin=273 K;pin=5 MPa

图3 网格划分Fig.3 Mesh display

3 结果与讨论

3.1 模型验证

为验证前述模型的准确性，设计矩形截面可视化喷管及相应的测试系统，开展Laval喷管内气体凝结相变实验研究。室内实验过程中，天然气存在较大危险，硫化氢是剧毒气体，实验室条件下难以采用真实的天然气组分进行实验研究，因此只能采用其他实验介质替代。结合实验室实际状况，采用湿空气为介质进行超声速流动凝结特性的研究。采用压力分布测试系统测量喷管内压力分布，多波长消光法颗粒测试系统测量喷管内凝结液滴数目和粒径，设计喷管总长为177.3 mm，喉部尺寸为8 mm×8 mm，入口尺寸为40 mm×8 mm，出口尺寸为12.51 mm×8 mm。实验工况为：入口压力0.574 MPa，入口温度292.95 K，入口湿度90.3%。

实验测得喷管沿程压力、液滴数量密度、液滴半径数据与模拟所得结果对比如图4所示。从图4(a)中可以看出，实验测得喷管沿程压力变化趋势与数值模拟结果基本一致，随着气体在喷管内的流动，压力逐渐降低，由于凝结液量非常小，凝结时释放的潜热对周围气体的加热作用被高速膨胀引起的温降效应所掩盖，并未捕捉到明显的凝结冲波现象。从图4(b)中可以看出，实验测得凝结液滴数目分布与模拟结果呈现出相同规律，液滴数目达到一定数量级后基本不变。从图4(c)中可以看出，实验测得凝结液滴粒径与数值模拟结果均在亚微米级。但由于喷管扩张段内气体流速达到超声速，测试装置对于液滴的捕捉存在一定困难。另外，虽然实验系统中采用了两级过滤器进行固体、液体颗粒的过滤，但仍无法完全消除外界颗粒对气体凝结的影响，使得实验测得数据与模拟结果之间存在一定偏差。以lg(Nexp/Nsim)衡量液滴数目之间的偏差在-1.78～1.01个数量级，与目前国内外常用气体凝结参数测试结果的偏差相当[17,27]，且模拟所得凝结起始位置与实验结果较为吻合，进一步说明了所建气体自发凝结数值模型能较为准确地描述喷管内气体的超声速凝结流动过程。

图4 喷管内凝结流动参数实验数据与模拟结果对比Fig.4 Comparison of flow and condensation parameters in nozzle between experimental data and simulation resultsTin=292.95 K;pin=0.574 MPa(a)p；(b)N；(c)r

3.2 凝结参数分布

基于上述凝结流动数学模型，模拟计算甲烷-硫化氢混合气体在喷管中的凝结流动过程，入口压力5 MPa，入口温度273 K工况条件下，模拟所得凝结参数沿喷管轴线分布如图5所示。

图5 凝结参数沿喷管中心线分布Fig.5 Distributions of condensation parameters along nozzle axisTin=273 K;pin=5 MPa(a)ΔT，J；(b)N，r；(c)Y，ys

由图5(a)看到，随着气流在喷管内的膨胀，硫化氢气体由不饱和状态逐渐进入过饱和状态，过冷度不断增加，最大可达40 K左右，气体处于极度热力学不平衡状态。从图5(a)可以看出，气体达到饱和状态后并不会立即发生凝结，而是达到一定过冷度后才出现大规模成核现象，凝结核在较短距离内急剧产生，成核率峰值可达1021数量级，并释放凝结潜热，过冷度迅速降低，成核环境遭到破坏，成核率由峰值剧降为0，成核过程结束。由于成核过程的急剧性，液滴数量密度也急剧升高至1015数量级，成核过程结束后，液滴数目保持不变，如图5(b)所示，但由于气体仍处于过冷状态，硫化氢气体在凝结核心表面继续凝结，使得液滴继续生长，液滴半径和液相质量分数持续增加，混合气体中硫化氢的摩尔分数相应降低，如图5(c)所示。

在喷管出口处，气相中硫化氢摩尔分数可由10%降至3.33%，显示了较好的硫化氢液化脱出能力，但尚未达到商品天然气质量指标，需考虑将超声速旋流分离技术与其他脱硫工艺相结合。当天然气中硫化氢摩尔分数较高时，可将超声速旋流分离技术作为第一步脱硫化氢措施，将硫化氢摩尔分数降至较低水平，然后再进入胺法等常规天然气脱硫装置进一步处理，使净化气达到规定指标，在满足产品处理指标的同时降低装置能耗，使传统技术与新技术的优势得以充分发挥。

3.3 入口压力对硫化氢凝结的影响

保持入口温度273 K不变，不同入口压力情况下喷管内过冷度、成核率、液相质量分数分布如图6 所示。从图6可以看出，更高入口压力情况下入口过冷度较大，硫化氢气体能更早达到大规模成核所需要的过冷状态，凝结成核过程发生在更靠近喷管喉部处，更多的硫化氢气体从气相中凝结出来，生成的液相质量分数增加。

凝结起始位置随入口压力变化情况如图7所示。由图7看到，随着入口压力的降低，凝结起始位置向后移动。在入口温度273 K条件下，当入口压力低至2 MPa时，在靠近喷管出口处才开始有凝结核心生成，硫化氢气体不能在喷管内完成自发成核过程，喷管冷凝效果几乎为0，此时超声速旋流分离装置将无法实现对天然气中硫化氢气体的脱除。

3.4 入口温度对硫化氢凝结的影响

保持入口压力5 MPa不变，不同入口温度条件下喷管内过冷度、成核率及液相质量分数分布如图8 所示。从图8可以看出，随着入口温度的升高，入口过冷度降低，更晚达到凝结成核所需要的极限过冷状态，成核起始点后移，且凝结生成的液相质量分数随入口温度的升高而减少。

图6 不同入口压力喷管内过冷度(ΔT)、成核率(J)及液相质量分数(Y)沿喷管中心线的分布Fig.6 Distributions of supercooled degree (ΔT),nucleation rate (J)and liquid-phase mass fraction (Y) in nozzle at different inlet pressures Tin=273 Kpin/MPa:(1)3；(2)4;(3)5(a)ΔT；(b)J；(c)Y

凝结起始点随入口温度变化情况如图9所示。由图9看到，随着入口温度的升高，凝结起始点后移，在入口压力5 MPa条件下，当入口温度高于313 K时，在靠近喷管出口处才开始有凝结核心生成，硫化氢气体无法在喷管内完成自发成核过程，此时喷管冷凝效果几乎为0，超声速旋流分离装置将无法实现对天然气中硫化氢气体的脱除。

图7 凝结起始点位置(X)随入口压力(pin)变化Fig.7 Nucleation onset points (X)vs inlet pressures (pin)Tin=273 K

图8 不同入口温度(Tin)喷管内过冷度(ΔT)、成核率(J)及液相质量分数(Y)分布Fig.8 Distributions of supercooled degree (ΔT),nucleation rate (J)and liquid-phase mass fraction (Y)in nozzle at different inlet temperatures (Tin)pin=5 MPaTin/K：(1)293；(2)283;(3)273(a)ΔT；(b)J；(c)Y

图9 凝结起始点位置(X)随入口温度(Tin)变化Fig.9 Nucleation onset points (X)vs inlet temperatures (Tin)pin=5 MPa

3.5 背压对硫化氢凝结的影响

保持入口温度(273 K)和压力(5 MPa)不变，不同背压情况下喷管内压力、过冷度、成核率和液相质量分数如图10所示。从图10可以看出，随着背压的增大，喷管内产生激波，并逐步向喷管喉部方向移动，激波引起压力突跃，破坏气体的凝结环境。当背压不是很大时，激波未到达成核区间，对成核过程不会产生影响，如图10(c)中(2)所示，但会造成凝结液滴再次蒸发，液相质量分数变为0，如图10(d)中(2)所示。当激波前移至成核区间时，成核过程受到影响，如图10各图中(1)所示，激波的产生导致压力、温度的突升，破坏喷管内的低温环境，过冷度迅速降低，成核环境遭到破坏，使得成核条件不再具备，凝结核心无法继续生成，喷管冷凝效果几乎为0。因此，在超声速旋流分离装置运行过程中，应合理控制背压，使激波不进入喷管，避免激波对凝结过程的破坏。

图10 不同背压条件(pb)下喷管内压力(p)、过冷度(ΔT)、成核率(J)及液相质量分数(Y)分布Fig.10 Distributions of pressure (p),supercooled degree (ΔT),nucleation rate (J)and liquid phase mass fraction (Y)in nozzle at different backpressures (pb)Tin=273 K;pin=5 MPapb/MPa:(1)3.2；(2)2.0;(3)0.5(a)p；(b)ΔT；(c)J；(d)Y

4 结 论

(1)气流进入喷管后高速膨胀，达到一定过冷度后硫化氢气体将发生凝结成核现象，成核过程在较窄区域内完成，此后液滴继续生长，液相质量分数不断增大，气相中硫化氢含量随之降低，在入口压力5 MPa、入口温度273 K条件下，喷管出口气相中硫化氢摩尔分数可由10%降至3.33%，但尚无法满足商品天然气要求，需与其他脱硫方法结合才能达到规定的处理指标。

(2)升高入口压力或降低入口温度将使硫化氢气体更早达到凝结起始条件，从而使成核起始点前移，在喷管出口能达到更大的液相质量分数，因此可通过调节来流入口压力和温度来促进硫化氢气体的凝结过程；过低的入口压力或过高的入口温度将使硫化氢气体无法在喷管内完成自发凝结过程，超声速旋流分离装置将无法实现天然气中硫化氢的脱除。

(3)随着背压的升高，激波在喷管内产生并逐渐前移，激波的产生会破坏凝结所需冷凝环境，造成凝结液滴的再蒸发，为避免激波对硫化氢气体凝结过程的影响，应合理控制背压，防止激波进入喷管。

符号说明：

dr/dt——液滴生长速率，m/s；

E——总能，J/kg；

h——气体总焓，J/kg；

hlv——凝结潜热，J/kg；

J——成核率，m-3·s-1；

kB——Boltzmann常数，1.3806505×10-23J/K；

keff——有效导热系数，W/(m·K)；

Kn——Knudsen数；

m——单位时间、单位体积内凝结的液滴质量，kg/(m3·s)；

mv——单个硫化氢分子质量，kg；

N——液滴数目，kg-1；

p——压力，Pa；

pb——背压，Pa；

pin——入口压力，Pa；

Prv——气体Prandtl数；

r——液滴半径，m；

rc——临界液滴半径，m；

S——过饱和度；

Sh——能量源相，J/(m3·s)；

Sm——质量源相，kg/(m3·s)；

Su——动量源相，kg/(m2·s2)；

SY——湿度源相，kg/(m3·s)；

t——时间，s；

T——气体温度，K；

Tin——入口温度，K；

Ts——饱和温度，K；

u——速度，m/s；

ui、uj——速度分量，m/s；

u′i、u′j——速度波动，m/s；

x——轴向坐标，m;

xi、xj——轴向与径向位置坐标，m;

X——凝结起始点位置，m;

y——径向坐标，m;

ys——硫化氢气相摩尔分数，%;

Y——液相质量分数；

γ——气体比热容比；

δij——Kronecker delta数；

ΔT——过冷度，ΔT=Ts-T，K；

ε——Lamanna给出的修正系数，0.01；

θ——无因次表面张力；

λv——导热系数，W/(m·K)；

μ——黏度，kg/(m·s)；

ρ——混合相密度，kg/m3；

ρs——硫化氢气体密度，kg/m3；

ρl——液相密度，kg/m3；

ρv——气相密度，kg/m3；

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

τeff——有效应力张量。