基于Matlab的C2,C3類(lèi)氣體壓縮機(jī)功率求解

沈順成，何 博，黎章杰

(武漢理工大學(xué) 機(jī)電工程學(xué)院，湖北 武漢 430070)

?

基于Matlab的C2,C3類(lèi)氣體壓縮機(jī)功率求解

沈順成，何 博，黎章杰

(武漢理工大學(xué) 機(jī)電工程學(xué)院，湖北 武漢 430070)

對(duì)于傳統(tǒng)功率計(jì)算不適合C2，C3類(lèi)壓縮機(jī)在臨界參數(shù)附近工作時(shí)的問(wèn)題，通過(guò)對(duì)比態(tài)原理和一系列實(shí)際氣體狀態(tài)方程，推導(dǎo)出絕熱指數(shù)k值的公式并推出更具普遍性的功率計(jì)算關(guān)系式，并通過(guò)Matlab編寫(xiě)出計(jì)算程序，成功運(yùn)用到C2，C3類(lèi)壓縮機(jī)的功率求解上。

C2和C3類(lèi)壓縮機(jī);絕熱指數(shù)；Matlab；功率求解

某石油化工廠的乙烯增壓壓縮機(jī)投入實(shí)際生產(chǎn)后出現(xiàn)超負(fù)荷運(yùn)行狀況，壓縮機(jī)額定功率為870 kW，驅(qū)動(dòng)機(jī)額定功率為1 000 kW，然而壓縮機(jī)實(shí)際運(yùn)行功率為1 149 kW，出現(xiàn)32%誤差。基于此，筆者用對(duì)比態(tài)原理、實(shí)際氣體Lee-Kesler狀態(tài)方程、starling狀態(tài)方程、Soave-Redlish-Kwang (RKS)方程和維里方程等基本狀態(tài)方程推導(dǎo)出功率計(jì)算中關(guān)鍵性的溫度絕熱指數(shù)的求解方法，這不僅可以解決查表法誤差大、應(yīng)用范圍小的問(wèn)題，而且由此可以推出適用更廣的功率計(jì)算關(guān)系式。借助Matlab畫(huà)出壓縮機(jī)運(yùn)行中溫度、壓縮性系數(shù)和絕熱指數(shù)的變化圖，并求解出理論總功率，不僅解決C2，C3類(lèi)壓縮機(jī)在功率求解上的空白，而且能形象直觀地了解和研究壓縮機(jī)運(yùn)行中不同氣體、不同工況下各指數(shù)的變化關(guān)系。

1 絕熱指數(shù)k的公式推導(dǎo)

經(jīng)過(guò)一系列演算[1]可得出過(guò)程指數(shù)與氣體常數(shù)、氣體比熱、氣體壓縮性系數(shù)、溫度的關(guān)系式：

(1)

令：

(2)

則式(1)可變?yōu)椋?/p>

(3)

由式(3)可知，要求出溫度絕熱，就必須先求出壓縮性系數(shù)Z和X之值。目前求Z值方法有很多，其中Lee-Kesler方程在0.3≤Tr≤4,0≤Pr≤10，范圍內(nèi)可使用,基本滿足所有工程實(shí)際，尤其適用于非極性和輕微極性氣體及其混合物。對(duì) C1到C6的烷烴及乙烯、丙烯、氮、二氧化碳和苯氣體的計(jì)算，在寬溫度和壓力范圍內(nèi)其計(jì)算出的壓縮因子Z的偏差為：過(guò)熱汽≤0.94%，飽和蒸汽≤1.02%[2]。L-K方程[3]如式(4)所示：

(4)

式中:Z為待求氣體的壓縮性系數(shù)；Z0為簡(jiǎn)單流體的壓縮性系數(shù)；ZR為參考流體的壓縮性系數(shù)；ω為待求氣體的偏心因子；ωR為參考流體的偏心因子。簡(jiǎn)單流體是指Ar等球形結(jié)構(gòu)分子的氣體，參考流體用正辛烷(C8H18)，其ωR=0.397 8。

在一定溫度和壓強(qiáng)下，ZR和Z0可以通過(guò)維里方程迭代求得，如下：

(5)

(6)

(7)

(8)

通過(guò)查表可查出實(shí)際待求氣體的偏心因子ω，由于Z計(jì)算十分復(fù)雜，筆者采用迭代解法，為了表述方便，不妨將該氣體在某一溫度T、壓強(qiáng)P下壓縮性系數(shù)Z的迭代解設(shè)為Z(T,P)，當(dāng)T，P已知時(shí)其為一確定值。

由式(2)可知，要求X，必須得到一個(gè)實(shí)際氣體Z與T的關(guān)系式，筆者采用對(duì)輕烴類(lèi)適用性強(qiáng)、精度高的RKS狀態(tài)方程進(jìn)行X的計(jì)算。RKS狀態(tài)方程[4]為：

(9)

(10)

(11)

(12)

m=f(w)=0.480+1.57w-0.176w2(13)

其中,Ωa=0.427 48,Ωb=0.086 64。

基本實(shí)際氣體狀態(tài)方程為：

Pv=ZRT

(14)

(15)

將X，Z代入式(3)即可求出某一溫度、壓強(qiáng)下的溫度絕熱指數(shù)：

(16)

2 功率公式推導(dǎo)

圖1 壓縮機(jī)工作時(shí)實(shí)際循環(huán)圖

要精確地表示出絕熱指數(shù)隨時(shí)間的變化關(guān)系是極其困難甚至是不可能的，因?yàn)槟承?shí)驗(yàn)數(shù)據(jù)[6]表明，在極其微小的時(shí)間段內(nèi)，絕熱指數(shù)甚至出現(xiàn)上下波動(dòng)的情況。但是宏觀上看，絕熱指數(shù)k的變化規(guī)律[7]還是可知的，為了畫(huà)出C2，C3類(lèi)氣體壓縮時(shí)的ab和cd段實(shí)際變化曲線，筆者作出以下假定：在一極小段的壓強(qiáng)變化內(nèi)氣體的絕熱指數(shù)是不變的，以此求解出功率p。即當(dāng)n無(wú)窮大時(shí)，p～p+p/n的壓強(qiáng)變化范圍內(nèi)，絕熱指數(shù)k是不變的。

已知進(jìn)氣壓力為Ps，排氣壓力為Pd，將壓強(qiáng)變化分為相等的n個(gè)區(qū)間，假定每個(gè)壓力區(qū)間內(nèi)絕熱指數(shù)是不變的，則每個(gè)區(qū)間的壓力變化量為：

(17)

那么第i個(gè)區(qū)間的起始點(diǎn)壓力為:

Pi=Ps+(i-1)ΔP

(18)

當(dāng)i=1時(shí)，即起始點(diǎn)處，此時(shí)P1=Ps，T1=Ts，則第一個(gè)區(qū)間的絕熱指數(shù)k1可按式(16)求出：

(19)

已知溫度變化公式為：

(20)

第一個(gè)區(qū)間結(jié)束處溫度即為第二個(gè)區(qū)間起始溫度，即：

(21)

則第二個(gè)區(qū)間的絕熱指數(shù)為：

(22)

通過(guò)每個(gè)區(qū)間起始點(diǎn)壓力算出每個(gè)區(qū)間起始溫度，然后依次算出每個(gè)區(qū)間的絕熱指數(shù):

(23)

(24)

對(duì)于每個(gè)區(qū)間，再根據(jù)式(25)和第一個(gè)區(qū)間起始點(diǎn)體積V1，可以依次算出每個(gè)區(qū)間的起始點(diǎn)體積Vi,如式(26)所示。

Pi-1Vi-1ki-1=PiViki-1=const

(25)

(26)

已知每個(gè)區(qū)間起始點(diǎn)壓強(qiáng)和體積，根據(jù)式(25)可以把每一個(gè)區(qū)間在壓縮機(jī)工作時(shí)的P-V循環(huán)圖表示出來(lái)，這樣就可以畫(huà)出n個(gè)小段連接出來(lái)的ab曲線，同理可以畫(huà)出cd曲線，然后通過(guò)Matlab就可以計(jì)算出abcd包圍的面積，從而求解出總功率。

3 Matlab編程求解

由上述壓縮因子Z、過(guò)程量X和總功率N的求法可知，其計(jì)算量很大，且部分運(yùn)用到數(shù)值分析、多元隱函數(shù)，必須借助計(jì)算機(jī)來(lái)進(jìn)行計(jì)算。Matlab是一套高性能的數(shù)值計(jì)算和可視化軟件，因此筆者采取Matlab編程來(lái)進(jìn)行求解和仿真驗(yàn)算，簡(jiǎn)化的計(jì)算流程圖如圖2所示。

圖2 簡(jiǎn)化的計(jì)算流程圖

以乙烯壓縮機(jī)進(jìn)行驗(yàn)算，乙烯壓縮機(jī)全負(fù)荷狀態(tài)下的基本數(shù)據(jù)如表1所示。

傳統(tǒng)功率計(jì)算式[8]為：

表1 乙烯壓縮機(jī)全負(fù)荷狀態(tài)下的基本數(shù)據(jù)

(27)

運(yùn)行程序，輸入初始參數(shù)，壓強(qiáng)等分?jǐn)?shù)n取100時(shí)運(yùn)行結(jié)果如圖3所示。程序計(jì)算的功率值N=1 072 917.207,即算出的設(shè)計(jì)功率約為1 073 kW，與1 149 kW的誤差只有6.6%，可見(jiàn)筆者提出的功率求解公式相比于傳統(tǒng)計(jì)算公式精確度提高了很多。但n取值較大時(shí)，程序運(yùn)行時(shí)間較長(zhǎng)。當(dāng)n=100時(shí)，運(yùn)行時(shí)間已超過(guò)1 min。

圖3 乙烯各參數(shù)變化圖

由于總功率是由各小區(qū)間功率相加得到，需要求n次積分，其計(jì)算量大，故筆者利用Matlab對(duì)折線ab進(jìn)行擬合，得到一條光滑曲線，那么只需求一次積分即可求出總功率。修改程序后再次運(yùn)行，得運(yùn)行結(jié)果N=1 072 917.497，并仿真出循環(huán)示意圖，如圖4所示。其中ab、cd表示擬合后的曲線，可見(jiàn)擬合后的功率基本沒(méi)變，且大大縮短運(yùn)算周期。

圖4 乙烯壓縮循環(huán)示意圖

由式(9)～式(15)可知，臨界參數(shù)、對(duì)比態(tài)參數(shù)(Tr，Pr)對(duì)X有著決定影響，乙烯壓縮機(jī)運(yùn)行時(shí)恰好是在乙烯的臨界溫度和臨界壓強(qiáng)附近工作，X值會(huì)大于0，從而導(dǎo)致絕熱指數(shù)k比零度時(shí)的絕熱指數(shù)要大，最終導(dǎo)致實(shí)際功率比按傳統(tǒng)功率計(jì)算的設(shè)計(jì)功率要大[9]。

為了驗(yàn)證筆者所提功率求解方法的正確性，把同是臨界參數(shù)附近工作的C2、C3類(lèi)氣體(如乙烷)帶入程序，輸入相同的初始參數(shù)運(yùn)行，得到圖5所示乙烷參數(shù)變化圖和圖6所示功率計(jì)算指示圖，功率N=1 075 719.662，而按傳統(tǒng)功率計(jì)算公式計(jì)算出的功率約為850 kW，可以得到與前所述相同的結(jié)論。

圖5 乙烷各參數(shù)變化圖

圖6 乙烷功率計(jì)算指示圖

通過(guò)公式推導(dǎo)可知，筆者所提功率公式也適用非C2,C3類(lèi)氣體，以氮?dú)?、氫氣為例，假定入口溫度?0 ℃，入口壓強(qiáng)為2.5 MPa，出口壓強(qiáng)為5.0 MPa，用傳統(tǒng)公式計(jì)算的功率和筆者方法解出的指示功率如表2所示。

由表2可知，對(duì)于介質(zhì)不在臨界參數(shù)附近工作的壓縮機(jī)，筆者方法解出的功率與傳統(tǒng)功率解相對(duì)誤差的絕對(duì)值不到1%，可見(jiàn)筆者所提出的功率求解法更具普適性。同時(shí)氫氣和氮?dú)飧鲄?shù)變化也完全符合實(shí)際情況，如圖7和圖8所示。

表2 傳統(tǒng)公式與筆者方法解出指示功率對(duì)比

圖7 氫氣各參數(shù)變化圖

圖8 氮?dú)飧鲄?shù)變化圖

4 結(jié)論

傳統(tǒng)功率計(jì)算認(rèn)定氣體的絕熱指數(shù)是不變的，如氫氣、氮?dú)狻⒖諝獾确荂2，C3類(lèi)氣體取0 ℃測(cè)得的絕熱指數(shù)計(jì)算功率的誤差較小[10]，但傳統(tǒng)功率公式并不適合所有氣體，C2，C3類(lèi)氣體壓縮機(jī)在臨界參數(shù)附近工作時(shí)，工作溫度和壓強(qiáng)皆在臨界參數(shù)附近，當(dāng)工作溫度高于0 ℃時(shí)，絕熱指數(shù)較傳統(tǒng)公式在0 ℃測(cè)得的絕熱指數(shù)會(huì)大很多，而且有一定的波動(dòng)，仍以0 ℃時(shí)的絕熱指數(shù)來(lái)進(jìn)行計(jì)算，會(huì)導(dǎo)致設(shè)計(jì)功率偏小。筆者推導(dǎo)出的功率解法，并不只適用于C2，C3類(lèi)氣體，而是對(duì)所有氣體都適用。不僅能很較地解決傳統(tǒng)功率公式不適合C2、C3類(lèi)壓縮機(jī)的難題，而且可以避免傳統(tǒng)公式計(jì)算中查圖表選取Z、k值的誤差，其精度也比傳統(tǒng)功率計(jì)算公式高。

[1] 郁永章.活塞式壓縮機(jī)[M].北京：機(jī)械工業(yè)出版社,1982：12-57.

[2] 李珍華，繆玲梅.實(shí)際氣體壓縮因子的應(yīng)用[J].上?；?2014，39(3)：7-10.

[3] 蘇洛養(yǎng).三參數(shù)法求氣體壓縮性系數(shù)的計(jì)算機(jī)解[J].壓縮機(jī)技術(shù),1990,100(2)：21-25.

[4] 王東宇,李岳,王淑蘭.R-K-S方程在壓縮因子計(jì)算中的應(yīng)用[J].化學(xué)工業(yè)與工程技術(shù),2001,22(5)：32-36.

[5] 林梅,孫嗣瑩.活塞式壓縮機(jī)原理[M].北京：機(jī)械工業(yè)出版社,1987:16-34.

[6] 張勇，陳浩，楊金.天然氣壓縮機(jī)多變過(guò)程指數(shù)的研究[J].機(jī)械,2008，35(1):27-33.

[7] 王君,李強(qiáng),李娟娟.往復(fù)壓縮機(jī)氣體壓縮熱力過(guò)程研究[J].流體機(jī)械,2008(1)：22-25.

[8] 活塞式壓縮機(jī)編寫(xiě)組.活塞式壓縮機(jī)設(shè)計(jì)[M].北京：機(jī)械工業(yè)出版社,1974：22-53.

[9] 石玉美，孫嗣瑩.超高壓乙烯的熱力學(xué)參數(shù)計(jì)算[J].流體機(jī)械,1994(11):22-29.

[10] 馮全科.往復(fù)式壓縮機(jī)指示功率計(jì)算公式的精度分析[J].壓縮機(jī)技術(shù),1989，95(3)：13-15.

SHEN Shuncheng:Assoc. Prof.; School of Mechanical and Electric Engineering, WUT, Wuhan 430070, China.

[編輯：王志全]

C2and C3Gas Compressor Power Solution Based on Matlab

SHENShuncheng,HEBo,LIZhangjie

In consideration of that the traditional power calculation is not suitable for the C2，C3 compressor near the critical parameters problems, we introduced a formula of adiabatic index k and a more general power calculation formula through a series of real gas state equation and the theorem of corresponding states. And the power of C2 C3 compressor can be successfully solved out by the Matlab calculation program.

C2C3compressor; adiabatic index; Matlab; power solution

2015-07-02.

沈順成(1959-)，男，湖北武漢人，武漢理工大學(xué)機(jī)電工程學(xué)院副教授.

武漢高新技術(shù)開(kāi)發(fā)基金資助項(xiàng)目(2014010101010010).

2095-3852(2015)06-0832-05

A

TH457

10.3963/j.issn.2095-3852.2015.06.037