活塞压缩机气流脉动数值模拟和实验验证.docx
《活塞压缩机气流脉动数值模拟和实验验证.docx》由会员分享，可在线阅读，更多相关《活塞压缩机气流脉动数值模拟和实验验证.docx（54页珍藏版）》请在冰豆网上搜索。
活塞压缩机气流脉动数值模拟和实验验证
活塞压缩机气流脉动数值模拟及实验验证
1、绪论
1.1研究背景及意义
活塞式压缩机广泛应用于石油、化工、冶金、天然气行业，作为一种重要的气体增压设备，在一些工艺流程中发挥着关键作用，这些设备能否正常运行直接关系到企业的生产能力[1]。
在持续安全生产中威胁最大的是管道振动，而管道振动的最大诱因就是气流脉动。
由于活塞式压缩机吸、排气的非连续性，不可避免使管道内气体压力出现周期性的波动，这就是气流脉动[1,2];活塞式压缩机管道系统都存在一定程度的气流脉动，这种脉动的压力在管道的突变截面、弯头、盲管、阀门等处产生交变的激振力，进而引发振动，工业现场经常出现剧烈的管道振动导致管路焊接处或法兰联接处振断，造成生产事故。
控制管道振动首先应准确掌握管道系统的气流脉动情况，尤其是管道系统中关键节点如气缸连接法兰、弯头、阀门等处的压力脉动幅值。
分析气流脉动的方法主要有两种，一种是平面波动理论，另一种是一维非定常可压缩流体流动理论[3]。
平面波动理论是研究气流脉动现象时最早发展起来的理论，这种方法做了几个方面的重要假定：
压力脉动值相对管道气流的平均压力值很小[4,5];气体遵守理想气体的性质;认为管道中气体流速相对声速小到可以忽略不计的程度[6]。
因此波动理论建立气体脉动的控制方程时能做线性化处理，最终得出能求解析解的波动方程。
在符合假定的条件下，波动理论能预测出符合实际的压力脉动幅值。
波动理论作出的假定在数学模型上就决定了它不能完整描述管道内压力波和非稳态流动耦合的复杂现象。
一般认为波动理论对气体与管道壁面摩擦考虑不足，导致其在脉动幅值较大尤其共振状态下计算值偏大。
此外波动理论在实际求解过程中将整个管道元件中的气流参数平均值取作气流参数值进行计算，这就决定了管道内气流参数值是常数而不是随实际状态变化的值，这降低了波动理论的模拟压力脉动的准确度。
非定常可压缩流动理论在建立描述管道内气流脉动现象的控制方程时，没有忽略非线性因素，综合考虑了气体与管道壁面的摩擦问题，实际气体性质的问题[2]。
而且多认为非定常可压缩流动理论在摩擦问题上处理的更符合实际，因而在脉动幅值较大的情况下计算值比波动理论更符合实测值。
但是摩擦阻尼能否显著抑制脉动幅值还有待进一步验证，其它影响气流脉动的因素还有哪些?
哪一个因素起了重要作用?
如何定量分析它们的影响?
这些问题目前研究的还不够。
此外，用非定常方法建立的双曲型控制方程组需要用数值方法求解，双曲型方程应用在压力脉动上会有哪些特性，数值求解的特点、如何获得较准确的收敛解，这些问题都有待进一步分析。
1.2气流脉动研究现状
气流脉动的研究是随压缩机工业的建立开始的，工程师很早就认识到这种现象对压缩机管道系统的重要影响，美国西南研究院自20世纪50年代已经展开气流脉动的理论和实验研究[7]。
1962年，Kinsl和Kfrey[8]最早提出经典的平面波动理论，至今仍是气流脉动研究的基础性理论之一[9]，波动理论不考虑管道内气流流速和气体实际性质，并忽略非线性因素，最终得出波动方程，从而用声波传播的原理很好的揭示了气流脉动的机理，对加深认识气流脉动的本质有重要意义。
气流脉动研究的两大任务是压力脉动幅值和气柱固有频率的计算，60年代后期有学者开始对压力脉动幅值计算进行初步探索[10,11]。
1970年，日本学者Toru等[12]提出转移系数法，用结构离散化的思想，将通常复杂的管道系统分割成不同的元件，分别计算。
这样处理的优点是易于实现数字计算机编程，因而得到了广泛应用，至今仍是脉动计算的主流方法之一。
70年代初山田荣[13]、野田桂一郎[14]提出刚度矩阵法，克服了转移系数法对分支管路处理繁琐的缺陷。
1973年酒井敏之等[15]提出计算复杂管系气柱固有频率的转移矩阵法，仍然借助结构离散化思想，首先计算每个管道元件的转移矩阵，再进行总装配，最后用计算机求解出各阶气柱固有频率，这种方法同样易于编程计算，因而应用非常广泛。
同年，美国的Sodel教授引入经典的亥姆霍兹共鸣器法，开始了压缩机消声器研究[16]。
以上几种方法都是基于波动理论发展起来的，而波动理论在阻尼因素上作了线性化处理即认为阻尼与速度成正比，当阻尼超出线性范围时，计算值比实际值偏大，因此限制了它的应用范围。
后来有研究人员[17]对波动理论进行改进，认为速度的平方决定摩擦力的大小，使波动理论能计算脉动幅值较大的情况，拓展了它的应用范围。
与此同时，不作简化直接用数值计算手段求解管道内非定常气流流动控制方程组的方法从70年代初开始，1972年Benson[18]总结了数值模拟方法的一些进展，提出可处理管道边界的匀熵特征线法。
自1974年起，在美国普渡大学历届召开的国际压缩机会议，都会讨论气流脉动项目，大大推动了此项研究。
这一年的会议上Singh和Sodel[19]教授共同发表一篇综述，全面总结了压力脉动和气柱固有频率计算的各种方法，制订出衰减压力脉动的评价标准。
同年，Elson[20]首次考虑了气阀阀片运动和管路压力波动的相互影响，为精确模拟压缩机吸、排气口处压力脉动情况打下基础。
随着计算机技术的进步，数值模拟的手段越来越受重视，1976年在普渡大学召开的国际压缩机会议上，Maclaren[21]等基于一维非定常流动理论，提出了较为完善的数学模型，建立的非线性双曲型方程组中考虑了气体与管道壁面的非线性摩擦问题以及管道截面变化的影响，得到与实测波形吻合程度较高的计算结果，验证了一维非定常可压缩流动数学模型应用于气流脉动模拟的可行性;文中对比了特征线法、Lax-Wendroff格式和Leap-Frog格式三种算法的数值计算结果，指出特征线法比后两种算法计算精度低，而且更容易衰减压力波的高频成分，但也指出特征线法是计算边界节点信息必不可缺的方法，文中还首次采用非匀熵特征线法计算边界节点，精度比匀熵特征线法高，该文对数值模拟气流脉动有巨大的指导意义。
此后，以Sodel[22]、Singh[23]为代表的研究人员在前人研究成果的基础上进一步取得进展，不断完善气流脉动的数学模型，将已经取得的成果推广到结构更复杂的多气缸大型压缩机上。
随着理论的不断成熟，20世纪80年代以后工程界侧重控制技术的研究[24-27]，并逐步形成了在石化、天然气工业界广泛认可的API618标准[28]，该标准由美国石油协会联合会员单位共同制订，详细规定了石化与天然气行业用压缩机气流压力脉动幅值上限和管道振幅允许值，并约定了分析气流脉动和管道振动的三种方法。
此标准的广泛认可也使压缩机制造商和用户越来越重视气流脉动问题，并积极开发控制技术。
美国西南研究院自2007年起，展开以声学衰减器为突破点的新一代压力脉动控制技术[29-31]，目前已经取得阶段性的成果。
国内是西安交通大学的党锡淇和陈守五教授等人最早发起气流脉动的研究。
从1974年开始着手，他们借鉴了国外转移矩阵法、转移系数法和刚度矩阵法的研究成果，并进一步发展：
推导出各种典型管道元件的转移矩阵[32]，在转移系数法中引入线性摩擦阻尼[33];对一维非定常流动也作了一定研究，推导出等截面管内气流的非稳态流动控制方程组，用匀熵特征线法处理容器、突变截面、汇流点等元件联接处，使数值计算得到简化[34,35];在理论分析的基础上进行了大量实验研究[36];在深入理论研究和大量工程实践的基础上总结出压力脉动的控制措施[37,38]，他们的研究成果集中体现在一本关于活塞式压缩机管道气流脉动与振动的专著上[39]。
近年来，国内学者进一步取得进展：
2001年，西安交通大学的彭学院教授基于平面波动理论开发出气流脉动分析软件，该软件能够计算任意复杂管系的气柱固有频率及各节点处压力脉动幅值，为快速分析压缩机管道系统声学特性提供了有效工具;2003年，李志博通过大量的实验验证了该软件计算结果的可靠性[40]。
近年来气流脉动的研究趋势表现在：
以美国西南研究院为代表侧重使用纳维斯托克斯方程一维流动模型建立描述管道内非稳态气流流动的控制方程，引入因粘性产生的气体与管道壁面的摩擦力，改变了以往一维非定常气流方程中摩擦力靠经验公式计算的方式[21]，方程同样需要有限元或有限差分的数值方法求解[41]，并将这种数值解法定义为时域分析法，将波动理论的解析解法定义为频域分析法，认为时域法比频域法作的假设更少，计算结果更符合实际，借助时域法还可以计算出因压力脉动造成的动态压力损失，进而帮助设计者改进压缩机整体性能。
因此认为时域法更有价值，投入了大量精力研究它的计算特性，探讨提高计算精度的方法。
另外也有研究人员[42,43]使用CFD软件运用三维流动理论模拟管道内气体的压力脉动，一般认为缓冲罐、气液分离器等三维结构特征明显的元件以及压缩机吸、排气口等复杂流道处三维方法的结果更准确，西安交通大学的徐斌[44]用Fluent软件在大脉动情况下获得了比一维方法更准确的结果，但也指出一维流动理论在小脉动时精度仍然很高。
以上研究现状的分析表明，基于一维流动的理论仍是分析压缩机管道气流脉动的有效方法，一维非定常流动理论是较为完善的数学模型，随着计算科学的进步，用数值解法精确模拟管道内流体运动越来越重要，但其计算特性如何;如何准确、可靠的得出结果;怎样用数值方法定量分析影响气流脉动的各种因素;摩擦阻尼是否有显著的影响;如何分析非定常方法和波动理论计算差异。
这些问题有待进一步探索，本文将在这些方面进行研究。
1.3本文所做工作
为了深入研究活塞式压缩机管道内气流脉动的机理，探索更加精确的模拟方法，在前人研究的基础上进一步认识气流脉动的内在规律，本文拟做以下几个方面的研究：
1）基于一维非定常可压缩流动理论建立描述活塞式压缩机管道内气流脉动现象的控制方程组，分析差分方程的稳定性条件，在用特征线法建立差分格式的过程中分析稳定性条件的物理意义。
2）编写一维非定常方法数值计算程序，通过大量的计算分析双曲型方程数值计算特性和程序的准确度、可靠性。
讨论影响计算结果准确度的主要因素，尤其是网格长度的影响。
3）搭建专门研究活塞式压缩机管道内气流脉动的实验台，测量管道不同位置处的压力脉动值。
通过与实验测量值对比，分析导致计算和实测差异的原因，指出数学模型上可改进之处;定量分析影响压力脉动波形和幅值的因素，尤其是摩擦阻尼的影响;分析导致波动理论方法和一维非定常方法计算差异的原因;定量评价局部阻力在变截面处抑制气流脉动的作用。
2、气流脉动的数学模型及求解
平面波动理论分析气流脉动时作了理想气体、等熵流动等假设，并且基本方程忽略了非线性项、气流平均流速的影响[45]。
为了在数学模型上更完整准确的描述脉动现象，本章建立一维非定常气流流动方程，着重考虑管路中的摩擦、实际气体性质等问题。
2.1一维非定常气流的守恒型方程组
由于实际输气管路管径与管长之比一般非常小，流体在同一截面上的各参数如压力、密度、速度等可以认为相等[39]，所以能够从一维的角度分析气流脉动现象。
在管道内取相邻两个截面形成的微团作为研究对象，推导连续方程、运动方程和能量方程。
得出一组描述一维非定常可压缩气流运动的偏微分方程。
这组方程可以表示成守恒型和非守恒型的形式[46]，在空气动力学数值计算上守恒型方程更受重视[47]。
气流在管路内作一维流动，则压力、速度、密度分别为坐标x和时间t的函数，即
2.1.1连续方程
1）通过控制面净流出控制体的流体质量
如图2-1所示取等截面管左侧I截面及相邻右侧II截面包围的空间为控制体，轴向长度取为dx。
在dt时间内由x截面气流流进的质量为
。
由x+dx截面气流流出的质量为：
则时间内通过I、II控制面净流出控制体的流体质量为：
2）控制体内流体质量的变化
在dt时间内控制体内流体质量的变化为：
3）流体流动的连续方程
根据质量守恒定律，可以得出以下关系式：
等截面管中横截面积s是常数，于是可以从上式消去sdxdt，则得等截面管内流体流动的连续方程：
它表示了对于非定常流动，单位时间净流出控制体的质量等于微元控制体内密度的变化。
2.1.2动量方程
如图2-2，在管道内仍取I、II截面内控制体为研究对象，控制体内流体的动量在t瞬时为
，在dt时间内的变化量为：
在
时间内，通过控制体的动量净流出量为：
另外，作用在截面I和截面II上的瞬时压力冲量代数和为：
考虑气体与管道壁面的摩擦，假设单位流体质量受到的摩擦力与瞬时速度的平方成正比，即：
式中：
——摩擦系数，
=
为壁面对气体的切应力[48];
D——管道内径/m。
由以上分析可知管内流体所受摩擦力主要与流体速度的二次方成正比，其摩擦系数与雷诺数有关[49]。
则可以写出控制体表面摩擦力的冲量：
根据动量定理可以得到以下数学表示式：
消去公因子
后得：
上式即为等截面管内气流的动量方程。
2.1.3能量方程
能量方程是对流动流体运用能量守恒定律得出的数学表达式。
在等截面管内任取一封闭控制面，其所包围的空间为控制体。
根据能量守恒定律，单位时间控制体内能量的变化量与控制体能量净流出量之和等于热交换的能量加上表面力所做的功。
所以能量方程的建立要考虑到以下因素[39]：
1）控制体能量净流出量
控制体内流体的能量由两部分组成：
宏观流体运动的动能和微观分子运动的动能（内能），对单位质量流体分别为
和，则控制体能量净流出量为：
2）控制体内能量的变化量
控制体内流体具有的能量为
。
在单位时间
内变化量为：
3）表面力所做功
表面力即瞬时压力，在单位时间
内，压力作的功为
，所以压力在控制面I、II上所作功的代数和为：
4）热交换
在单位时间内，设单位质量流体与外界的热交换量为
，输入热量取正值，输出时取负值。
假定气体流动过程中，气流与管道壁面的摩擦力作负功，全部转变为热量并且都被气体吸收，则摩擦力作功与气体吸收的热量平衡，在
时间内，控制体与外界的热交换量为：
综合以上分析，根据能量守恒定律得出等截面管内气流的能量方程为：
化简后得：
上式中
表示单位长度管道气流与壁面的换热量，它可由下式计算出：
式中：
——换热系数/W·（m2·K）-1;
——管道内径/m;
——管内流体温度/K;
——管外环境温度/K。
换热系数
的计算公式如下：
式中：
——流体的导热系数/W·kg·K）-1;
——努塞尔数。
它可用雷诺数Re和普朗特数Pr求得[49]：
方程（2-1）、（2-3）、（2-4）独立的未知物理量有
、
四个，而方程组只有三个，为使方程封闭，需要补充内能表达式，根据热力学理论，气体内能的计算式为：
式中：
上式即为控制体内气流的能量方程。
将连续方程（2-1）、动量方程（2-3）和能量方程（2-9）联立写成矩阵形式为：
上式为一维非定常气流的守恒型方程组，为将其无因次化，引入以下参数：
——参考压力/Pa;
——参考密度/kg·m-3;
——参考长度/m;
——参考声速/m·s-1。
则各参数的无因次表示式各为：
——无因次声速;
——无因次速度;
——无因次压力;
——无因次密度;
——无因次坐标;
——无因次时间。
由于
则连续方程的无因次表达式为
同理，动量方程的无因次表达式为：
能量方程（2-9）为：
将（2-11）、（2-12）、（2-13）仍然写成矩阵形式为：
引入符号
于是方程组（2-14）化为：
上式是一个非线性的一阶双曲型偏微分方程组，它的解要用近似的数值方法求得。
2.2一维非定常守恒型方程组的数值解法
2.2.1有限差分法
方程组无法获得解析解，必须采用近似的数值方法。
有限差分法是求解偏微分方程最常用的数值解法之一，其基本原理是：
在积分域内用有限的数值差商代替极限形式的微商，将连续问题离散化，最终化成有限形式的线性代数方程组[50]。
用差分法将偏微分方程组离散化的步骤是，首先在求解区域作网格划分，对于一维情形是把x区间分成一些等距或不等距的小区间即空间步长，用有限数目的网格节点代替连续的求解区域[51];然后将原微分方程组转化成差分形式的方程组;最后从已知的初始值开始，按照一定时间步长沿时间轴逐步推算，直至符合设定的精度[39]。
2.2.2一阶双曲型方程组差分格式
差分格式的构造与偏微分方程的特征及解的性质有关，由于特征型方程的两大优点：
（1）便于反映物理意义
（2）便于边界处理[50]。
所以在得出双曲型方程组的差分格式之前先引入关于特征的一些概念。
将一维非定常气流的守恒型方程组写成如下形式：
方程称为一维非定常气流的特征型方程组，下面对方程中的第一式进行分析：
由式（2-22）知
，即u沿直线值L保持不变，这种直线是特征线[52]。
图2—3是a>0和a<0时的特征线示意图。
沿特征线，方程可以化为常微分形式，而且波则沿着特征线以有限速度a传播。
因为波速是有限值，所以存在依赖区域和影响区域，这些特点对双曲型方程的数值求解很重要[50]。
上式中C为常数。
为确保时间增加时，解
有界，也必须使V有界，即Re（a）<0,则a>0，因此对方程,为使计算稳定，若a>0,则用空间向后差分近似
,反之a<0,则空间导数应向前差分，否则不稳定。
前的系数a表示波运动的速度，a>0表明波是沿x轴正方向运动，这时要用向后差分的格式来近似空间一阶导数才能保证差分格式条件稳定。
由于差分指向与波前进方向刚好相反，所以称迎风或逆风，如图24所示，可见迎风格式与特征线的方向相关[52]。
当k>0时，对式构造右偏心的迎风差格式：
下面用特征线方法构造本文所用的Lax-Wenrodff差分格式，令a>0，特征线方向和网格如图26所示，假定第n时间层值
已知，要计算第n+1时间层p点的（m,n+1）值
。
过P点作特征线与n时间层相交于Q点，若CFL条件成立，即点在线段BC上。
根据特征线上参数值保持不变的特点可知
都已得到，因此可以在B、C、D三点作抛物型插值来求出的值，从而得到的值[52]。
将非线性双曲方程的两步Lax-Wendroff差分格式应用在无因次化后的一维非定常气流守恒型方程组：
2.2.3边界节点的特征线法计算
图27展示了Lax-Wendroff两步法用到的网格。
第一步，从Z时刻节点1和节点2的信息计算出Z+1/2△Z时刻节点4的信息，同样从节点2和节点3计算出节点5;第二步，根据节点4和节点5的信息，计算出Z+△Z时刻节点6的结果。
但此差分格式仅适用于计算管道内部节点，不能计算管道端点即边界点，因为用Z时刻各节点的信息去计算Z+△Z时刻的信息时，需要用到相邻节点的信息，而对边界节点缺乏相邻点的信息，所以不能计算[2]。
边界节点和的计算需要借助特征线法[21]。
将方程式（2-20）的方程式单独列出：
根据特征线的性质，沿特征线方程化为常微分关系式。
方程
于是特征型方程组
设满足稳定条件的时间步长是△Z，它与网络均分距离△X构成计算网络。
在管路的始端和末端两个边界节点需要用特征线法处理。
将式和中的特征线视为直线，其斜率由所在点Z时刻的U和A值来确定。
沿第一特征线的黎曼变量的变化量的求法是：
综上所述，在等截面管的内部节点使用精度较高的两步Lax-Wendroff法，在边界节点上使用匀熵修正理论给出的特征线法计算[34]。
2.2.4边界条件
管道端点分内外两种。
与外界相连的端点称为外端点，主要有三类：
开口端、闭口端和压缩机端，管道不同单元间连接的端点称为内端点，如突变截面联接点、容器联接点等[39]。
下面用特征线法说明外端点与内端点的处理方法。
用U1、U2分别表示节点1、2在t时刻的无因次速度，
表示节点1、2在t+△t时刻的无因次速度，密度、压力、黎曼变量等参数的表示方式与此相同。
首先介绍两个重要的关系式：
根据管内气流方向或指定方向，左端点用第二特征线求β;右端点用第一特征线求λ。
管道左、右端点的特征线如图2-9所示。
1）外端点
（1）闭口端
盲管、关闭的阀门处、压缩机气阀关闭时都是闭口边界，闭口端速度为0，U=0;根据和求出密度和压力：
（2）开口端
管道端部与大气连通或者压力为定值都是开口端，端点处压力是常数即常数;其它参数求法如下：
（3）压缩机端
与气缸相连的管道端点认为是此边界条件[54]。
气阀开启后，受气阀阀片运动和活塞运动的影响，压缩机端的气流运动是非常复杂的，为简化求得与气缸相连的管道端点的气流速度，作出以下假定[55]：
a）气阀的关闭与打开都是瞬间进行，于是忽略掉阀片运动对气流流动的影响;
b）气阀打开后，吸、排气口处气流速度与活塞速度成相关，其相关性系数是活塞面积与管道横截面积的比值[56]。
基于以上假定可以得出气阀开启的时间内，压缩机端气流的无因次速度为[39]：
2.3数值计算程序设计
本文数学模型是以等截面管道内气流为分析对象建立起来的，所以只能对等截面管道划分网格进行计算。
而实际管道系统还有突变截面、容器、孔板、阀门等元件，这些非等截面管道元件将整个管路分割成不同长度的管段，管段内部都是用二阶精度的差分格式计算，管段与非等截面管道元件的联接处，流动条件发生变化，必须用特征线法处理才能计算出联接处的参数，使下游管路的计算得以继续进行。
实际管道系统要分段计算，各元件有独自的结构参数，又必须单独处理它们，所以设计程序时是用结构离散化的思路，对管道系统
升级会员 