黏度是计算流体流动性质的重要物性参数，是天然气物性计算的重点之一，尤其是纯组分黏度的预测精度，对天然气混合物黏度的预测精度具有很大的影响[1]。黏度计算的理论模型主要有两类。第一类为基于对比态原理的关联式，Carr等[2]和Dean等[3]将天然气黏度表示成对比压力、对比温度、对比密度和相对分子量的函数；Lohrenz等[4]将黏度表示为密度的函数，其黏度计算结果对密度的取值非常敏感，尤其在高压下，单质处于液相或者超临界状态，密度计算偏差容易导致黏度计算出现严重误差；TRAPP法[5]和Pedersen关联式[6]需迭代求解，不易编程，且计算耗时很长。第二类为拟立方状态方程的黏度模型，Patel-Teja模型[7]在临界点不能还原为原PT状态方程形式，对高压超临界状态流体预测误差比较大[8]。其他方法例如采用牛顿内摩擦定律公式[9]计算黏度，方法简单，但误差也较大[10]。美国国家标准与技术研究院(NIST)研发的权威工质物性计算软件reference fluid properties(REFPROP)，被很多研究项目用作物性数据源或计算结果准确性的参考数据源[11]，可以准确预测流体的黏度，然而它虽然是开源的，但黏度程序代码极其复杂，涉及的介质数量庞大。
《Matheson气体数据手册》提出单质气体黏度的一元非线性回归计算式[12]，即基于常压下黏度的实验数据对温度唯一因变量进行的回归分析，对常压下单质黏度的预测比较准确，但在加压条件下预测误差却很大。一般天然气长途输送的管道压力，国外大多为8~12 MPa，我国西气东输管道压力为10 MPa[13]，压力对黏度的影响不容忽略，如不考虑高压对于天然气黏度的影响，将使高压下黏度预测出现较大偏差，影响工艺和设备的正确设计。在工程实践中，没有严格的理论可以将黏度表示成温度和压力的函数[14]，同时也没有一种二元回归模型可以准确描述黏度与温度、压力的关系。
为了快速、有针对性地获得天然气混合物的黏度，本文基于多元函数微积分计算原理，提出一种新的二元回归分析法，建立黏度与温度、压力的函数式。该二元回归分析法使用REFPROP软件中相关天然气纯组分在不同温度、压力下的黏度数据作为数据源，通过计算机编程来实现纯组分黏度的准确预测。
1 黏度的回归分析计算 1.1 天然气纯组分分类本文方法所研究的天然气纯组分温度和压力的适用范围为：230 K≤T≤450 K，p≤20 MPa。在适用范围内，根据纯组分黏度曲线的变化规律，以临界温度Tc为基准，将天然气混合物中各组分分为3类。
第一类，Tc≤230 K，例如H2、N2、CO、CH4、AR。如图 1所示为CH4的p-T图，图中方框内为本文所研究的纯组分温度压力的适用范围(下同)。在此范围内，纯组分呈气相或超临界状态，黏度曲线变化有明显规律性。
|
(Tc, pc)—临界点；l—液相区；g—气相区。 图 1 CH4的p-T图 Fig.1 Variation of pressure with temperature for CH4 |
第二类，Tc≥450 K，例如iC5H12、nC5H12、nC6H14、H2O，图 2为iC5H12的p-T图。在适用范围内，纯组分呈气相或者液相(忽略固相)，故黏度曲线有阶跃点存在，需将气相区和液相区分别进行研究。
|
(Tc, pc)—临界点；l—液相区；g—气相区。 图 2 iC5H12的p-T图 Fig.2 Variation of pressure with temperature for iC5H12 |
第三类，230 K＜Tc＜450 K，例如CO2、H2S、C2H4、C2H6、C3H8、iC4H10、nC4H10、C5H10，图 3为CO2的p-T图。在适用范围内，纯组分呈气相、液相或者超临界状态(忽略固相)，黏度曲线因气、液相变化而产生阶跃，在纯组分的临界点附近，黏度曲线变化趋势较复杂。
|
(Tc, pc)—临界点；l—液相区；g—气相区。 图 3 CO2的p-T图 Fig.3 Variation of pressure with temperature for CO2 |
由于第三类纯组分的黏度回归分析篇幅较长，将择篇详述。因此，本文针对前两类天然气纯组分的黏度，分别以CH4和iC5H12为例进行回归分析。
《Matheson气体数据手册》给出的气体黏度计算式将纯组分常压状态的黏度数据拟合成温度为唯一因变量的一元二次回归计算式
| $\mu = A + BT + C{T^2}$ | (1) |
本文提出的新的二元回归法克服常规二元回归模型无法准确描述黏度与温度、压力关系的困难，采用以下步骤进行二元回归分析：
1) 对纯组分分类并制作出纯组分在定温条件下的黏度随压力变化的曲线簇；
2) 将μ-p曲线变化不连续的纯组分分为气相区和液相区处理；
3) 分析曲线簇变化规律，将其拟合为幂次相同以压力为唯一变量的一元非线性方程组；
4) 将一元非线性方程组的各个系数分别拟合为温度的函数。
将上述拟合的计算式整合，则得到纯组分的二元回归计算式。新的二元回归分析法的计算框图如图 4所示。
|
图 4 二元回归分析法计算框图 Fig.4 Calculation block diagram for the binary regression analysis |
在适用状态范围内，CH4的黏度随温度、压力变化具有规律性，将黏度表示成温度和压力的函数
| $\mu = f\left( {T,p} \right)$ | (2) |
在REFPROP中规律提取不同温度、压力状态下CH4的黏度，并将黏度散点绘制成μ-p图(图 5)。
|
图 5 定温条件下CH4黏度随压力变化曲线簇 Fig.5 Variation of viscosity with pressure at different constant temperatures for CH4 |
目前，还没有一种二元非线性回归模型将图 5中的规律散点用函数μ=f(T, p)表示出来，为此，本文提出了新的二元回归分析法对图 5散点进行回归分析(使用Excel软件)。
在同一温度下，黏度随压力变化趋势有明显规律性：当T≥270 K，变化曲线与三次函数曲线贴合；当230 K≤T＜270 K，变化曲线与五次函数曲线贴合。
当T≥270 K时，将曲线簇中所有定温条件下μ-p黏度曲线的回归方程设为三次函数组
| $\mu = a\left( T \right){p^3} + b\left( T \right){p^2} + c\left( T \right)p + d\left( T \right)$ | (3) |
其中，a(T)、b(T)、c(T)、d(T)为只与温度有关的系数(表 1)，绘制a(T)-T散点图(图 6)，并进行数据回归，回归方程如式(4)，相关系数R2=0.9995。
|
下载CSV 表 1 CH4黏度的二元回归方程系数 Table 1 Coefficients of regression equations for CH4 |
|
图 6 CH4黏度的拟合系数a(T)-T散点图 Fig.6 a(T)-T scatter diagram of coefficients for CH4 |
| $\begin{array}{l} \quad \quad a\left( T \right) = - 3.756188 \times {10^{ - 16}}{T^6} + 8.537517 \times {\rm{ }}\\ {10^{ - 13}}{T^5} - 8.0629185 \times {10^{ - 10}}{T^4} + 4.0505275 \times {10^{ - 7}}{T^3} - {\rm{ }}\\ 1.14186674 \times {10^{ - 4}}{T^2} + 1.713282 \times {10^{ - 2}}T - \\ 1.069494 \end{array}$ | (4) |
同理得到b(T)、c(T)、d(T)的回归方程
| $\begin{array}{l} \quad \quad b\left( T \right) = - 1.2894206 \times {10^{ - 12}}{T^5} + 2.4727304 \times \\ {10^{ - 9}}{T^4} - 1.893064156 \times {10^{ - 6}}{T^3} + 7.237689 \times {10^{ - 4}}{T^2} - \\ 0.13836997T + 1.0608058 \times 10 \end{array}$ | (5) |
| $\begin{array}{l} \quad \quad c\left( T \right) = 6.737514 \times {10^{ - 12}}{T^5} - 1.2881283 \times {10^{ - 8}}\\ {T^4} + 9.82659 \times {10^{ - 6}}{T^3} - 3.740359 \times {10^{ - 3}}{T^2} + \\ 0.710612T - 5.3813975 \times 10 \end{array}$ | (6) |
| $d\left( T \right) = 0.02982731T + 2.21397$ | (7) |
当230 K≤T＜270 K时，将曲线簇中所有定温条件下μ-p黏度曲线的回归方程设为五次函数组
| $\begin{array}{l} \quad \quad \mu = {a_1}\left( T \right){p^5} + {b_1}\left( T \right){p^4} + {c_1}\left( T \right){p^3} + {d_1}\left( T \right){p^2} + \\ {e_1}\left( T \right)p + {f_1}\left( T \right) \end{array}$ | (8) |
同理可得系数回归方程
| $\begin{array}{l} \quad \quad {a_1}\left( T \right) = - 7.881 \times {10^{ - 10}}{T^3} + 6.3542 \times {10^{ - 7}}{T^2} - \\ 1.70976 \times {10^{ - 4}}T + 1.5354703 \times {10^{ - 2}} \end{array}$ | (9) |
| $\begin{array}{l} \quad \quad {b_1}\left( T \right) = 3.158926 \times {10^{ - 8}}{T^3} - 2.57598 \times {10^{ - 5}}{T^2} + \\ 7.012673 \times {10^{ - 3}}T - 6.37434 \times {10^{ - 1}} \end{array}$ | (10) |
| $\begin{array}{l} \quad \quad {c_1}\left( T \right) = - 3.5061024 \times {10^{ - 7}}{T^3} + 2.933751 \times \\ {10^{ - 4}}{T^2} - 8.191852 \times {10^{ - 2}}T + 7.6353283 \end{array}$ | (11) |
| $\begin{array}{l} \quad \quad {d_1}\left( T \right) = 5.4082757 \times {10^{ - 7}}{T^3} - 5.4539805 \times \\ {10^{ - 4}}{T^2} + 1.76634 \times {10^{ - 1}}T - 1.857905 \times 10 \end{array}$ | (12) |
| $\begin{array}{l} \quad \quad {e_1}\left( T \right) = 2.2967767 \times {10^{ - 6}}{T^3} - 1.528944 \times {10^{ - 3}}\\ {T^2} + 3.2302563 \times {10^{ - 1}}T - 2.079767 \times 10 \end{array}$ | (13) |
| ${f_1}\left( T \right) = 0.0374889T + 0.0172501$ | (14) |
式(3)~(14) 为CH4的黏度回归分析计算式，适用于温度、压力范围在230 K≤T≤450 K，p≤20 MPa内任意状态点黏度的预测。
1.3 iC5H12黏度回归分析 1.3.1 iC5H12黏度随压力变化曲线簇在REFPROP中规律提取不同状态点下iC5H12的黏度，并将黏度与温度、压力的关系绘制成μ-p散点图(图 7)。
|
图 7 定温条件下iC5H12的μ-p散点图 Fig.7 Scatter diagrams of viscosity with pressure at different constant temperatures for iC5H12 |
iC5H12的临界温度为460.35 K，在230~450 K温度范围内不会进入超临界状态。假设温度不变，单质相态随压力升高会由气相变为液相，产生图 7中的阶跃。如何界定气液两相？如何将阶跃的散点拟合成二元非线性回归方程？这都将成为二元回归分析的重点。
在REFPROP中规律提取iC5H12在不同压力下对应的露点温度，如表 2所示，将表 2中数据绘制成p-T散点图(图 8)。
|
|
下载CSV 表 2 iC5H12不同压力下对应的露点温度 Table 2 Dew point temperatures for different pressures of iC5H12 |
|
图 8 iC5H12气液相变曲线 Fig.8 Phase transition curve of iC5H12 |
iC5H12相变曲线的回归方程为
| $\begin{array}{l} \quad \quad {p_{\rm{k}}}\left( T \right) = 9.676740213 \times {10^{ - 13}}{T^5} - 9.559293079 \times \\ {10^{ - 10}}{T^4} + 5.3759899 \times {10^{ - 7}}{T^3} - 1.81907976 \times {10^{ - 4}}{T^2} + \\ 3.152401235 \times {10^{ - 2}}T - 2.111413 \end{array}$ | (15) |
式(15) 可以精确区分iC5H12在给定状态点(T0, p0)的相态：p0＜pk(T0)，为气相；p0≥pk(T0)，则为液相。iC5H12为气相时，绘制在定温条件下黏度随压力变化的μ-p散点图，如图 9所示。
|
图 9 定温条件下iC5H12黏度随压力变化散点图 Fig.9 Scatter diagrams of viscosity with pressure at different temperatures of iC5H12 |
iC5H12为气相时，按照1.2节中黏度回归分析方法进行拟合，则iC5H12气相黏度拟合计算式为
| $\mu = {a_2}\left( T \right){p^2} + {b_2}\left( T \right)p + {c_2}\left( T \right)$ | (16) |
| $\begin{array}{l} \quad \quad {a_2}\left( T \right) = - 2.9313902 \times {10^{ - 9}}{T^4} + 4.2638864 \times \\ {10^{ - 6}}{T^3} - 2.302761 \times {10^{ - 3}}{T^2} + 5.4701302 \times {10^{ - 1}}T - \\ 4.80281 \times 10 \end{array}$ | (17) |
| $\begin{array}{l} \quad \quad {b_2}\left( T \right) = 3.4774952 \times {10^{ - 9}}{T^4} - 5.1822019 \times \\ {10^{ - 6}}{T^3} + 2.8739075 \times {10^{ - 3}}{T^2} - 7.019954 \times {10^{ - 1}}T + \\ 6.38183314 \times 10 \end{array}$ | (18) |
| ${c_2}\left( T \right) = 2.4392 \times {10^{ - 2}}T - 4.77 \times {10^{ - 3}}$ | (19) |
由图 7可以看出：当T＜370 K时，黏度随压力变化趋近于线性；当370≤T≤450 K时，黏度随压力变化曲线为非线性。按1.2节黏度回归分析方法进行拟合，当T＜370 K时，线性回归分析为
| $\mu = {a_3}\left( T \right)p + {b_3}\left( T \right)$ | (20) |
| $\begin{array}{l} \quad \quad {a_3}\left( T \right) = - 3.6384523 \times {10^{ - 7}}{T^3} + 4.1929373 \times \\ {10^{ - 4}}{T^2} - 1.58136816 \times {10^{ - 1}}T + 2.1937288 \times 10 \end{array}$ | (21) |
| $\begin{array}{l} \quad \quad {b_3}\left( T \right) = 8.633715 \times {10^{ - 7}}{T^4} - 1.18338615 \times \\ {10^{ - 3}}{T^3} + 6.137439 \times {10^{ - 1}}{T^2} - 1.4403343 \times {10^2}T + \\ 1.3143756 \times {10^4} \end{array}$ | (22) |
当T≥370 K时，二次回归分析
| $\mu = {a_4}\left( T \right){p^2} + {b_4}\left( T \right)p + {c_4}\left( T \right)$ | (23) |
| $\begin{array}{l} \quad \quad {a_4}\left( T \right) = - 1.9852273 \times {10^{ - 8}}{T^4} + 3.2790985 \times \\ {10^{ - 5}}{T^3} - 2.0305526 \times {10^{ - 2}}{T^2} + 5.5860807T - \\ 5.759877 \times {10^2} \end{array}$ | (24) |
| $\begin{array}{l} \quad \quad {b_4}\left( T \right) = 5.2696212 \times {10^{ - 7}}{T^4} - 8.700847 \times {10^{ - 4}}\\ {T^3} + 5.38654 \times {10^{ - 1}}{T^2} - 1.4815692 \times {10^2}T + \\ 1.5276611 \times {10^4} \end{array}$ | (25) |
| $\begin{array}{l} \quad \quad {c_4}\left( T \right) = - 1.0181722 \times {10^{ - 4}}{T^3} + 1.263649 \times \\ {10^{ - 1}}{T^2} - 5.31594637 \times 10T + 7.643228 \times {10^3} \end{array}$ | (26) |
式(15)~(26) 为iC5H12黏度在适用范围(230 K≤T≤450 K，p≤20 MPa)内的二元回归计算式。
2 黏度拟合数据分析 2.1 与数据源REFPROP对比误差分别把以CH4、iC5H12为例的黏度拟合计算公式编入Fortran程序，在适用温度压力范围内，随机取17组状态点(T0，p0)输入程序，将输出结果与NIST数据软件REFPROP结果进行对比，如表 3所示。
|
|
下载CSV 表 3 CH4、iC5H12黏度拟合计算结果与数据源REFPROP对比 Table 3 Accuracy in predicting the viscosity of CH4 and iC5H12 compared with REFPROP |
以相对误差δ和平均相对误差R(AARE)检验计算结果精确度，具体计算方法为
| $\delta = \frac{{|{\mu _{{\rm{cal}}}} - {\mu _{{\rm{exp}}}}|}}{{{\mu _{{\rm{exp}}}}}} \times 100\% $ | (27) |
| $R = \left( {100/n} \right)\sum\limits_j {\frac{{|{\mu _{{\rm{cal}}}} - {\mu _{{\rm{exp}}}}|}}{{{\mu _{{\rm{exp}}}}}}} $ | (28) |
式中，μcal表示黏度拟合计算值，μexp表示数据源REFPROP提取值或实验值，n表示状态点数。
由表 3可以看出，新的二元回归分析计算式的拟合误差结果为：CH4的相对误差最小0.03%，最大0.68%，平均0.20%；iC5H12的相对误差最小0.04%，最大0.84%，平均0.29%。两者误差都在1%以内，随着压力的升高，预测精度稳定，而且纯组分iC5H12的气液相处理十分准确。
表 4为基于新的二元回归分析方法，针对17组随机状态点的分析结果。所获得的9种天然气纯组分计算黏度与REFPROP软件数据平均误差均小于1%，表明二元回归计算式计算结果精确。
|
|
下载CSV 表 4 黏度预测与数据源REFPROP对比 Table 4 Accuracy in predicting the viscosity compared with REFPROP |
对比5种纯组分黏度的实验数据[15]，应用回归分析法拟合的黏度预测精度见表 5。由表 5可以看出，在新方法适用范围内，天然气纯组分(H2，N2，AR，CO，CH4)的二元回归分析计算式平均相对误差均小于1%，预测结果精确，表明采用新的二元回归分析方法能够满足工程精度要求。
|
|
下载CSV 表 5 黏度预测与实验数据对比 Table 5 Accuracy in predicting the viscosity compared with experimental data |
(1) 基于多元函数微积分计算原理，提出了针对天然气纯组分黏度预测的二元回归分析法，将天然气纯组分的黏度表示成温度、压力的函数，使得计算式在高压状态也适用，适用范围230 K≤T≤450 K，p≤20 MPa。
(2) 新的二元回归分析法预测结果与数据源REFPROP和实验值的平均相对误差均小于1%，精确度较高。该方法易于编程应用实现天然气纯组分黏度的快速计算，并可进行二次开发。若通过应用合理混合法则，可以实现天然气黏度的准确计算，这对于天然气压缩机的正确设计和三维流场的准确数值分析都是非常重要的。
致谢: 衷心感谢辽宁重大装备制造协同创新中心对本文研究工作的支持！| [1] |
刘夏天. 天然气压缩机多组分混合工质物性计算研究[D]. 辽宁大连: 大连理工大学, 2014. Liu X T. The study of multi-component mixture physical properties calculation for natural gas compressor[D]. Dalian, Liaoning:Dalian University of Technology, 2014.(in Chinese) |
| [2] |
Carr N L, Kobayashi R, Burrows D B. Viscosity of hydrocarbon gases under pressure[J]. Journal of Petroleum Technology, 1954, 6(10): 47-55. DOI:10.2118/297-G |
| [3] |
Dean D E, Stiel L I. The viscosity of nonpolar gas mixtures at moderate and high pressures[J]. AIChE Journal, 1965, 11(3): 526-532. DOI:10.1002/(ISSN)1547-5905 |
| [4] |
Lohrenz J, Bray B G, Clark C R. Calculating viscosities of reservoir fluids from their compositions[J]. Journal of Petroleum Technology, 1964, 16(10): 1171-1176. DOI:10.2118/915-PA |
| [5] |
Huber M L, Hanley H J M. The corresponding-states principle:dense fluids[M]//Millat J, Dymond J H, Nieto de Castro C A. Transport properties of fluids. Their correlation, prediction and estimation. Cambridge, UK:Cambridge University Press, 1996.
|
| [6] |
Pedersen K S, Fredenslund A. An improved corresponding states model for the prediction of oil and gas viscosities and thermal conductivities[J]. Chemical Engineering Science, 1987, 42(1): 182-186. DOI:10.1016/0009-2509(87)80225-7 |
| [7] |
王利生, 郭天民. 基于Patel-Teja状态方程的统一粘度模型(Ⅱ)应用于油气藏流体粘度的预测[J]. 化工学报, 1993, 44(6): 685-691. Wang L S, Guo T M. A unified viscosity model based on Patel-Teja equation of state(Ⅱ) Application in predicting viscosity of reservoir fluid[J]. Journal of Chemical Industry and Engineering(China), 1993, 44(6): 685-691. (in Chinese) |
| [8] |
荣淑霞, 郭绪强, 杨继涛, 等. 烃类混合物及油气藏流体的粘度计算模型[J]. 中国石油大学学报:自然科学版, 1999, 23(1): 79-82. Rong S X, Guo X Q, Yang J T, et al. The viscosity model of hydrocarbon mixture and reservoir fluid[J]. Journal of China University of Petroleum:Natural Sciences, 1999, 23(1): 79-82. (in Chinese) |
| [9] |
魏凯丰, 宋少英, 张作群. 天然气混合气体黏度和雷诺数计算研究[J]. 计量学报, 2008, 29(3): 248-250. Wei K F, Song S Y, Zhang Z Q. Research on calculation of natural mix gas viscosity and Reynolds number[J]. Acta Metrologica Sinica, 2008, 29(3): 248-250. (in Chinese) |
| [10] |
王霞, 卢坤. 天然气混合气体粘度计算方法[J]. 延安大学学报:自然科学版, 2011, 30(2): 54-55. Wang X, Lu K. The calculation method of natural gas mixture viscosity[J]. Journal of Yanan University:Natural Science, 2011, 30(2): 54-55. (in Chinese) |
| [11] |
董情, 蔡博. 基于REFPROP高温热泵工质性能研究[J]. 日用电器, 2015(8): 127-130. Dong Q, Cai B. Performance study for high temperature heat pump working fluid based on REFPROP[J]. Electrical Appliances, 2015(8): 127-130. (in Chinese) |
| [12] |
卡尔 L约斯. Matheson气体数据手册[M]. 7版. 北京: 化学工业出版社, 2003. Yaws C L. Matheson gas data book[M]. 7th ed. Beijing: Chemical Industry Press, 2003. (in Chinese) |
| [13] |
陈绍凯, 李自力, 雷思罗, 等. 高压天然气压力能的回收利用技术[J]. 煤气与热力, 2008, 28(4): 1-6. Chen S K, Li Z L, Lei S L, et al. Recovery and utilization of pressure energy from high-pressure natural gas[J]. Gas & Heat, 2008, 28(4): 1-6. (in Chinese) |
| [14] |
Sanjari E, Nemati Lay E, Peymani M, et al. An efficient reliable method for calculation of natural gas viscosity[J]. Petroleum Science & Technology, 2014, 32(11): 1300-1308. |
| [15] |
石油化学工业部化工设计院. 氮肥工艺设计手册:理化数据[M]. 北京: 石油化学工业出版社, 1977.
|