李 希，叶 晔，刘 恒，康国强，苏万兵，姚 帅

（1.西安石油大学，陕西西安 710065；2.中国石油长庆油田分公司第三采气厂，内蒙古乌审旗 017300；3.中国石油长庆油田分公司苏里格南作业分公司，陕西西安 710018）

1 清管现状

清管器在管道中的定位显得尤为重要，以下为某采气厂某天然气处理厂的清管速度计算如公式（1）所示：

清管器运行时的瞬时速度计算公式为：

式中：v-清管器瞬时速度，m/s；Qn-始发站流量，104m3/d；d-输气管线内径，m；p-清管器后压力，MPa。

此公式经一次计算便得出全程运行速度和预计行进时间，而真实情况是清管器在管道内情况多变，管道部分部位结垢或沉积物较多，会造成清管器行进速度变慢，从而影响实际上的到达时间。且始发站流量Qn与压力p均不能保持绝对恒定，也会对计算结果造成误差，从而导致计算值与实际运行时间相差较大，往往并不具有参考意义（见表1）。

由表1总结可知，其理论耗时10.3 h与实际耗时29.5 h相差近三倍。

目前某采气厂各处理厂均无法在清管过程中确定清管器的实时位置，所以导致收球工作过早或过晚，耗费人力物力。

综合以上清管球球速计算与清管球运行里程计算的各项不足，现提出利用现有的数字化管理平台PKS系统进行改进，可较为精确的预估清管器在管道中的实时位置及实时运行速度[1-3]。

2 理想气体的状态方程

理想气体状态方程是描述理想气体在处于平衡态时，压强、体积、物质的量、温度间关系的状态方程，公式（2）：

表1 X干线收球总结

式中：P-气体绝对压力，Pa；T-气体绝对温度，K；n-气体物质的量，mol；R-理想气体常数。

由此公式变换可得公式（3）：

假设很短一段时间区间t内P、T变化忽略不计，则如公式（4）：

根据公式（3）可知，当已知工况温度T1、标况温度T2、工况压力P1、标况压力P2与标况体积流量V2时即可得出工况条件下体积V1。根据此原理得出初步方案：如公式（4）所示利用PKS系统上干管的实时流量Q2（标况）、压力P2（标准大气压）与平均地温值 T2（即气体温度）。计算出当压力为干管压力P1时的工况流量Q1。则Q1体积便为干管压力下时间区间t内流经管道的体积V。

则根据公式（5）已知管道横截面积S便可计算出清管器在管道内行进的速度V。

3 实际气体状态方程

气体状态方程是表示气体温度、压力和比容之间相互联系的函数。理想气体被看作分子之间没有力的相互作用且气体分子本身没有体积。输气管道压力高达几兆帕甚至十几兆帕，此时天然气的性质与理想气体的性质有很大差异，单纯利用理想气体状态方程计算球速时，误差依然较大。为进一步减少计算误差，精确实时定位球速与位置，现引入实际气体的状态方程，用压缩因子来表示实际气体与理想气体的差异，实际气体的状态方程表示为公式（6）：

式中：p-气体绝对压力，Pa；T-气体绝对温度，K；V-气体的比容，m3；n-气体物质的量，mol；Rg-气体常数，8.314 J（/mol·K）。

Z值与气体成分相关，而且随温度和压力的改变而变化，故不同种类气体的Z=（fP，T）曲线是不同的。

由于压缩因子Z在任意温度和压力下均不相同，则为精确计算球速，需要求出每一时间区间内的Z值，现引入BWRS方程如公式（7）：

式中：p-系统压力，kPa；T-系统温度，K；ρ-气相或液相密度，kmol/m3；R-气体常数，8.314 3 kJ（/kmol·K）。

公式中 A0、B0、C0、D0、E0、a、b、c、d、α、γ 是 BWRS状态方程的11个参数。对于某个纯组分i的各项参数均可由它的临界参数Tci，pci和偏心因子wi由以下公式（8）求得：

式中参数A1～A11和B1～B11的值（见表2），这是Starling通过正构烷烃采用多种热力性质分析（PVT，焓和蒸气压）关联得到的。

表2 通用参数Ai、Bi的值

表2 通用参数Ai、Bi的值（续表）

部分天然气组分的 Tci、ρci和 ωi值（见表3）。

表3 部分天然气组分的Tci、ρci和ωi值

计算气体混合物性质时的混合规则都是半经验性的，是实验数据与理论计算的结果反复比较试探取得的。BWRS方程应用于混合物计算时，11个参数按下列公式（9）混合规则计算：

以上各式中xi，xj分别为气相混合物中i、j组分的摩尔分数；Kij为i、j组分间的交互作用系数，它表示i、j组分混合物和理想混合物所发生的偏差，Kij越大说明偏离理想状态越远。Kij的值（见表4）。

BWRS方程的 11个参数计算：首先利用式（8）计算流体中各个纯组分对应的11个参数，然后再利用式（9）的混合规则求出所求混合物对应的 A0、B0、C0、D0、E0、a、b、c、d、α、γ 等 11 个参数。

4 压缩因子Z的计算

利用BWRS公式计算压缩因子。

压缩因子Z利用公式（10）求出：

其中：A0、B0、C0、D0、E0、a、b、c、d、α、γ 等 11 个参数的计算方法与公式（9）相同。

计算气体密度，BWRS方程可看作压力关于密度的函数P=F（ρ），可选用抛物线法求解，迭代方程如下：

式中：下标j表示迭代序号，先设两个密度的初值，再用正割法求解。对于气体初值的选取可假设为理想气体，设迭代计算到为止。当收敛指标为ερ=10-4时，一般ρ的迭代次数只需要3～6次（注意：此时求出的密度单位为kmol/m3）。求出密度ρ之后，代入公式即可求出压缩因子Z。

表4 组分间的交互作用系数Kij

则将每一时段求得的Z=Z1代入公式（6）中，且设Z2为大气压下的压缩因子，近似取值为1，即可得公式（12）：

式中：Q1-工况下的流量，m3/h；Q2-标况下的流量，m3/h；Z1-工况下的压缩因子；P1-工况下的压力；P2-标况下的压力；T1-工况下的温度，K；T2-标况下的温度，K。

在输入界面中输入标况下的体积流量，管道内径，温度，压力等参数，点击计算，则可在输出界面中得出精确的工况下体积流量、清管器滑行距离与压缩因子Z，其中压缩因子Z作为计算的中间量代入公式继续计算得出前两个结果。

表5 X干线参数

根据X干线参数（见表5），试求某一时刻清管器的运行速度（见图1）。

从编程运行界面可以看出，当输入标况下体积流量45 000 m3/h、管道内径693.6 mm、温度15℃、压力2.3 MPa时。运算得出结果为工况下体积流量1 857.0 m3/h、清管器速度4 905.78 m/h、压缩因子Z=0.953。能够较为精确的描绘清管球在管道内的状态。由于Z=0.953，可见理想气体与高压的实际气体之间的确存在着较大差距，所以为了提高计算精确度，引入压缩因子Z是十分必要的。

图1 清管球速计算结果

5 总结

（1）相较于人工粗略计算清管器运行时间，软件实时精确计算提高了精度，减少了出错几率，提高了清管工作的效率。

（2）可实时监控球速，球速变化时还能记录下球速变化所在管道的位置。帮助清管人员对管道全程状态进行大致的了解。例如通过定位球速降低的某个位置，大致了解该位置管道的积液或腐蚀情况。

（3）通过多次清管并记录后，可较精确的描绘出整个管段上可能存在较为严重的腐蚀或结垢。在通过低洼管段时可大致了解该管段内的积液情况，从而更好的维护天然气管道。