石油与天然气化工  2017, Vol. 46 Issue (6): 1-7
超声速喷管内CO2气体凝结特性研究
曹学文 , 赵西廓 , 孙文娟     
中国石油大学(华东)储运与建筑工程学院
摘要:基于欧拉——欧拉双流体模型,建立气相及液相流动控制方程组,结合凝结成核与液滴生长理论,对喷管内CO2气体的凝结特性进行了数值模拟研究。结果表明,采用的数学模型和数值计算方法可较准确地反映喷管内气体的凝结流动过程。CO2气体凝结潜热较小,凝结冲波现象不明显;气体进入喷管特别是在经过喉部之后,在马赫数增大的同时,压力和温度降低,过冷度增加,最大可至30 K左右,并于凝结发生后快速下降至约5 K;CO2气体成核过程在时间和空间上表现出急剧性。凝结起始位置距喉部约2.21 mm,成核率由0激增至2.04×1021 m-3·s-1,液滴数目达到1015的数量级;凝结核心形成后,气体分子在一定的过冷度下在液滴表面团聚、液化,液滴半径和湿度迅速增加。成核过程结束后,已有凝结核心仍能不断生长,至喷管出口处液滴半径增至1.46×10-7 m,湿度可达0.093 5。
关键词超声速    CO2气体    喷管    凝结    酸性组分    天然气    
Study on condensation properties of CO2 gas in supersonic nozzle
Cao Xuewen , Zhao Xikuo , Sun Wenjuan     
College of Pipeline and Civil Engineering, China University of Petroleum(East China), Qingdao, Shandong, China
Abstract: Based on Euler-Euler two-fluid model, the gas and liquid flow control equations were established. Combined with the theory of condensation nucleation and droplet growth, the condensation properties of CO2 in a Laval nozzle were simulated. The results showed that both the mathematical models and numerical calculation methods used could reflect the process of gas condensation in the Laval nozzle accurately. The condensation wave of CO2 was not obvious due to its small latent heat of condensation. As the gas flow into the nozzle, the Mach number increased while pressure and temperature decreased, the undercooling increased up to 30 K, and then declined quickly to about 5 K after the occurrence of condensation. The process of nucleation showed a sharp both in time and space. Wilson point is about 2.21 mm from the throat, the nucleation rate increased sharply from 0 to 2.04×1021 m-3·s-1, while the number of droplets reached the order of 1015 instantly. After the formation of condensation core, gas molecules under certain degree of undercooling accumulated and liquefied on the droplets surface, causing droplet radius as well as humidity to increase rapidly. After the process of nucleation, the existing coagulation cores kept growing with the droplet radius increasing to 1.46×10-7 m and the humidity up to 0.093 5 at the exit of the Laval nozzle.
Key Words: supersonic    CO2 gas    Laval nozzle    condensation    sour components    natural gas    

天然气是一种清洁、高效、丰富和廉价的主要能源,占全球一次能源消费的三分之一,加速天然气勘探开发,发展天然气产业,对保障我国能源安全、推动经济增长及改善环境具有重大意义。随着天然气在我国能源结构中地位的不断上升,天然气净化工业的地位也越发重要。针对高酸性组分天然气及高CO2含量二氧化碳驱油田伴生气的处理,迫切需要开发新的技术。

近年来,超声速旋流分离技术开始应用于天然气处理加工领域,国内外研究学者先后将其应用于天然气脱水、脱重烃和液化中,利用超声速流动条件下气体的低温凝结效应,结合旋流分离技术,实现天然气中某些组分的冷凝分离[1-15]。该技术结合了气体动力学、工程热力学和流体力学理论,将膨胀降温、旋流式气液分离、再压缩等处理过程集中在一个密闭紧凑的装置中完成,具有密闭无泄漏、无需化学药剂、结构紧凑轻巧、节能环保等优点,特别适用于高酸性组分天然气的净化处理。在这个过程中,酸性组分在超声速喷管内的凝结是实现其脱除的前提和关键。因此,以天然气中的主要酸性组分CO2为例,开展了超声速喷管内CO2气体凝结特性的研究,对CO2气体在喷管内的凝结流动过程进行数值模拟,分析凝结流动参数的分布规律。该研究对于今后进一步开展天然气酸性组分超声速旋流分离机理研究,完善超声速旋流分离技术,丰富和发展天然气净化处理工艺具有重要的指导意义。

1 数学模型与数值计算方法
1.1 凝结成核模型

采用Girshick等[16-17]推导的内部一致经典成核理论(ICCT)进行成核率的计算,见式(1)。

$J = \frac{1}{S}\frac{{\rho _{\rm{v}}^2}}{{{\rho _{\rm{l}}}}}\sqrt {\frac{{2\sigma }}{{{\rm{ \mathsf{ π} }}m_{\rm{o}}^3}}} {\rm{exp}}\left( { - \frac{{4{\rm{ \mathsf{ π} }}\sigma r_{\rm{c}}^3}}{{3{k_{\rm{B}}}T}}} \right){\rm{exp}}(\theta )$ (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

1.2 液滴生长模型

采用Gyarmathy液滴生长模型计算液滴生长速率[18],见式(2)。

$\frac{{{\rm{d}}{r_{\rm{d}}}}}{{{\rm{d}}t}} = \frac{{{\lambda _{\rm{v}}}}}{{{\rho _{\rm{l}}}{h_{{\rm{lv}}}}}}\frac{{({T_{\rm{s}}} - T)\left( {1 - \frac{{{r_{\rm{c}}}}}{{{r_{\rm{d}}}}}} \right)}}{{{r_{\rm{d}}}\left( {1 + \frac{{2\sqrt 8 {\rm{ \mathsf{ π} }}}}{{1.5P{r_{\rm{v}}}}} \cdot \frac{\gamma }{{1 + \gamma }}Kn} \right)}}$ (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。

1.3 流动控制方程组

基于欧拉——欧拉双流体模型,忽略气、液相间速度滑移,同时认为液滴温度为当地压力下的饱和温度,建立气相及液相流动控制方程组。气相包括质量、动量及能量方程,见式(3)~式(5);液相包括质量、液滴数目及液滴半径方程,见式(6)~式(8)。

$\frac{{\partial {\rho _{\rm{v}}}}}{{\partial t}} + \frac{\partial }{{\partial {x_{\rm{j}}}}}\left( {{\rho _{\rm{v}}}{u_{\rm{j}}}} \right) = {S_{\rm{m}}}$ (3)
$\begin{array}{l} \quad \quad \quad \frac{\partial }{{\partial t}}({\rho _{\rm{v}}}{u_{\rm{i}}}) + \frac{\partial }{{\partial {x_{\rm{j}}}}}({\rho _{\rm{v}}}{u_{\rm{j}}}{u_{\rm{i}}}) = - \frac{{\partial {\rho _{\rm{v}}}}}{{\partial {x_{\rm{i}}}}} + \\ \frac{\partial }{{\partial {x_{\rm{j}}}}}\left[ {\eta \left( {\frac{{\partial {u_{\rm{j}}}}}{{\partial {x_{\rm{i}}}}} + \frac{{\partial {u_{\rm{i}}}}}{{\partial {x_{\rm{j}}}}} - \frac{2}{3}{\delta _{{\rm{ij}}}}\frac{{\partial {u_{\rm{j}}}}}{{\partial {x_{\rm{i}}}}}} \right)} \right] + \frac{\partial }{{\partial {x_{\rm{j}}}}}\left( { - {\rho _{\rm{v}}}\overline {u_{\rm{i}}^\prime u_{\rm{j}}^\prime } } \right) + {S_{\rm{u}}} \end{array}$ (4)
$\begin{array}{l} \frac{\partial }{{\partial t}}\left( {{\rho _{\rm{v}}}E} \right) + \frac{\partial }{{\partial {x_{\rm{j}}}}}\left( {{\rho _{\rm{v}}}{u_{\rm{j}}}E + {u_{\rm{j}}}{p_{\rm{v}}}} \right) = \\ \quad \quad \frac{\partial }{{\partial {x_{\rm{j}}}}}\left( {{k_{{\rm{eff}}}}\frac{{\partial T}}{{\partial {x_{\rm{j}}}}} + {u_{\rm{i}}}{\tau _{{\rm{eff}}}}} \right) + {S_{\rm{h}}} \end{array}$ (5)
$\frac{\partial }{{\partial t}}\left( {\rho Y} \right) + \frac{\partial }{{\partial {x_{\rm{j}}}}}\left( {\rho {u_{\rm{j}}}Y} \right) = {S_{\rm{Y}}}$ (6)
$\frac{{\partial \rho N}}{{\partial t}} + \frac{\partial }{{\partial {x_{\rm{j}}}}}\left( {\rho N{u_{\rm{j}}}} \right) = J$ (7)
${r_{\rm{d}}} = \sqrt[3]{{3Y/\left( {4{\rm{ \mathsf{ π} }}{\rho _{\rm{l}}}N} \right)}}$ (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=-mSu=-muSh=m(hlvh);SY=m

1.4 表面张力模型

由于凝结成核模型中表面张力σ以立方的形式出现,其微小变化可能对成核率的计算带来重大影响,进而影响数值计算结果的准确度。因此,有必要选择更为准确的表面张力计算方法。由于微小液滴的表面张力难以测得,根据平面液体表面张力对比计算值和实验值,一般认为,基于对比态原理和等张比容的表面张力计算方法最为精确[19]。因此,本研究优先选用对比态关联式进行表面张力的计算,该关联式适用于非极性液体,见式(9)。

$\sigma = \frac{{p_{\rm{c}}^{2/3}T_{\rm{c}}^{1/3}Q{{(1 - {T_{\rm{r}}})}^{11/9}}}}{{1{\rm{ }}000}}$ (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}}}}}$,为气体在正常沸点下的对比温度。

1.5 数值计算方法

UDS为用户自定义标量(User-Defined Scalar)的简称,Fluent软件允许用户自己定义一些标量来解决标准模块中不能解决的问题。本研究采用UDS输运方程添加液相流动控制方程组。

UDF为用户自定义函数(User-Defined Function)的简称,它是Fluent软件提供的一个用户接口,用户可通过UDF,以Fluent软件为平台进行二次开发和个性化设置。本研究利用C语言编写了相应的UDF程序,将源相表达式嵌入Fluent软件中进行计算。考虑真实气体效应,在UDF中对真实气体密度、比热容、黏度、导热系数、凝结潜热及液滴密度等参数进行计算,查阅相关文献选择合适方法。

Laval喷管中的气体流动属于高速可压缩流动,采用密度基方法求解。湍流模型选用k-ω方程。气、液相流动控制方程组、湍流动能方程、湍流耗散率方程均采用二阶迎风格式进行离散。进口设置为压力进口,指定总压、静压、总温及湍流参数。出口设置为压力出口,因超声速流动的所有流动属性均从内部推算得到,故不进行相应设置。固壁设置为无滑移、无渗流、绝热边界。

2 喷管内CO2气体凝结特性

通过目前文献中已有的喷管内CO2气体自发凝结实验数据,验证前述数学模型及数值计算方法,并对喷管内的CO2气体凝结特性进行分析。喷管结构及实验数据选自文献[20]。由于CO2气体液相区范围较窄,实验时难以保证其完全处于液相区而没有固相出现,故模型中不对液固两相进行区分。

2.1 喷管物理模型

实验所用喷管总长132.08 mm,喉部半径1.35 mm,位于35.56 mm处。其结构及尺寸见图 1。采用非结构化网格对喷管进行网格划分,并进行无关性验证。网格总数确定为7 338,网格划分情况见图 2

图 1     喷管结构及尺寸示意图 Figure 1     Structure and size schematic diagram of nozzle

图 2     喷管网格划分 Figure 2     Mesh methodology of nozzle

2.2 模型验证

选取文献[20]中nozzle I Run 1/2/5/6共4组实验数据,按照前述模型和方法进行数值计算,将数值计算所得的喷管内压力分布与实验结果进行对比。4组实验中喷管的入口参数列于表 1

表 1    实验中喷管入口参数 Table 1    Nozzle inlet parameters of different experiments

喷管内CO2凝结数值计算压力分布与实验结果对比见图 3。由图 3可知,数值计算所得喷管内压力变化趋势与实验结果基本一致,在可接受的误差范围内,数值计算值略小于实验值,表明所采用的数学模型和数值计算方法可较准确地反映喷管内气体自发凝结流动过程。值得说明的是,数值计算和实验结果均未捕捉到明显的凝结冲波现象,即由于蒸汽凝结导致蒸汽流量减小,释放大量潜热,喷管内气体在凝结的瞬间出现压力和温度的突跃。经分析,喷管内气体的自发凝结与喷管结构、气体性质、入口参数及环境条件等多种因素有关。其中,气体性质、入口参数对凝结冲波有较大影响。在前述实验工况下,CO2气体凝结潜热较小,凝结时释放的热量少,对周围气体的加热作用不明显,因而未对周围气体的压力和温度产生显著影响。从实验数据来看,喷管内压力下降较平稳,没有出现突跃,与数值计算结果基本一致,也从侧面证明了上述分析的合理性。

图 3     喷管内CO2凝结数值计算压力分布与实验结果对比 Figure 3     Comparison of pressure distribution between numerical simulation and experimental results of CO2 condensation in nozzle

2.3 计算结果分析

以Run 1实验工况为例,对喷管内CO2气体自发凝结过程中凝结流动参数的分布进行分析。

图 4(a)~(d)为喷管内气体压力、马赫数、温度和过冷度的分布。CO2气体进入喷管后,随着喷管的收缩及扩张,气体马赫数不断增大,压力和温度逐渐降低,过冷度增加。特别是气体经过喉部之后,压力和温度迅速降低,过冷度急剧增加,最大可至约30 K,产生自发凝结现象。凝结发生后,过冷度快速下降,于喷管出口降至约5 K。

图 4     喷管内CO2气体凝结流动参数分布 Figure 4     Flow and condensation parameters distribution of CO2 in nozzle

图 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。

图 5     气体凝结参数沿喷管轴线的分布 Figure 5     Condensation parameters distribution of gas along nozzle axis

3 结论

本研究基于欧拉——欧拉双流体模型建立气相及液相流动控制方程组,结合凝结成核与液滴生长理论,在对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。

参考文献
[1]
LAMANNA G. On nucleation and droplet growth in condensing nozzle flows[D]. The Netherlands: Eindhoven University of Technology, 2000. https://www.narcis.nl/publication/RecordID/oai%3Alibrary.tue.nl%3A539104
[2]
JASSIM E, ABDI M A, MUZYCHKA Y. Computational fluid dynamics study for flow of natural gas through high-pressure supersonic nozzles: part 1. real gas effects and shockwave[J]. Petroleum Science and Technology, 2008, 26(15): 1757-1772. DOI:10.1080/10916460701287847
[3]
JASSIM E, ABDI M A, MUZYCHKA Y. Computational fluid dynamics study for flow of natural gas through high-pressure supersonic nozzles: Part 2. nozzle geometry and vorticity[J]. Petroleum Science and Technology, 2008, 26(15): 1773-1785. DOI:10.1080/10916460701304410
[4]
KARIMI A, ABDI M A. Selective dehydration of high-pressure natural gas using supersonic nozzles[J]. Chemical Engineering and Processing: Process Intensification, 2009, 48(1): 560-568. DOI:10.1016/j.cep.2008.09.002
[5]
MALYSHKINA M M. The procedure for investigation of the efficiency of purification of natural gases in a supersonic separator[J]. High Temperature, 2010, 48(2): 244-250. DOI:10.1134/S0018151X10020161
[6]
曹学文, 陈丽, 林宗虎, 等. 用于超声速旋流分离器中的超声速喷管研究[J]. 天然气工业, 2007, 27(7): 112-114.
[7]
曹学文, 肖荣鸽, 杜永军, 等. 喷管内高速流动天然气相变特性数值研究[J]. 西安石油大学学报(自然科学版), 2008, 23(1): 56-60.
[8]
马庆芬. 旋转超音速凝结流动及应用技术研究[D]. 大连: 大连理工大学, 2009. http://cdmd.cnki.com.cn/Article/CDMD-10141-2009116110.htm
[9]
蒋文明. 多组分凝结性超音速流传热传质理论及实验研究[D]. 北京: 北京工业大学, 2010. http://cdmd.cnki.com.cn/Article/CDMD-10005-2010108019.htm
[10]
文闯. 湿天然气超声速旋流分离机理研究[D]. 青岛: 中国石油大学(华东), 2014. http://industry.wanfangdata.com.cn/yj/Magazine?magazineId=sydxxb&yearIssue=2010_4
[11]
杨文. 天然气高速膨胀凝结相变特性研究[D]. 青岛: 中国石油大学(华东), 2016.
[12]
程霖, 额日其太, 计维安, 等. 天然气超音速分离器中漩涡发生器及喷管的数值模拟研究[J]. 石油与天然气化工, 2011, 40(3): 232-235.
[13]
刘清, 王超. 天然气超音速旋流脱水装置设计及凝结特性分析[J]. 石油与天然气化工, 2016, 45(3): 1-5.
[14]
张明, 王春升. 超音速分离管技术在海上平台的应用分析[J]. 石油与天然气化工, 2012, 41(1): 39-42.
[15]
高晓根, 计维安, 刘蔷, 等. 超音速分离技术及在气田地面工程中的应用[J]. 石油与天然气化工, 2011, 40(1): 42-46.
[16]
GIRSHICK S L, CHIU C P. Kinetic nucleation theory: A new expression for the rate of homogeneous nucleation from an ideal supersaturated vapor[J]. The Journal of Chemical Physics, 1990, 93(2): 1273-1277. DOI:10.1063/1.459191
[17]
GIRSHICK S L. Comment on: self-consistency correction to homogeneous nucleation theory[J]. Journal of Chemical Physics, 1991, 94(1): 826-827. DOI:10.1063/1.460309
[18]
GYARMATHY G. The spherical droplet in gaseous carrier streams: review and synthesis[J]. Multiphase Science and Technology, 1982, 1(1-4): 99-279. DOI:10.1615/MultScienTechn.v1.i1-4
[19]
童景山. 流体热物性学[M]. 北京: 中国石化出版社, 2008.
[20]
D K M. Non-equilibrium condensation of carbon dioxide in supersonic nozzles[J]. Revista Chilena De Pediatría, 1966, 61(2): 74-77.
[21]
丁国良, 黄东平. 二氧化碳制冷技术[M]. 北京: 化学工业出版社, 2007.
[22]
吴玉国, 陈保东. BWRS方程在天然气物性计算中的应用[J]. 油气储运, 2003, 22(10): 16-21. DOI:10.3969/j.issn.1000-8241-D.2003.10.004
[23]
李玉星, 姚光镇. 输气管道设计与管理[M]. 2版. 东营: 中国石油大学出版社, 2009.