殷雄，溫凱*，謝萍，閃向營(yíng)，夏子杰，宮敬*

1 中國(guó)石油大學(xué)(北京)城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 102249 2 中石油管道有限責(zé)任公司西部分公司，烏魯木齊 830013*通信作者, ydgj@cup.edu.cn, kewin1983@126.com

1 天然氣管道流動(dòng)計(jì)算模型介紹

管道是長(zhǎng)距離運(yùn)輸天然氣的主要設(shè)備，天然氣管道系統(tǒng)是由多個(gè)壓氣站及其調(diào)節(jié)單元組成的高度非線性、強(qiáng)耦合、大時(shí)滯分布參數(shù)系統(tǒng)。天然氣管道內(nèi)氣體的流動(dòng)可以用連續(xù)性方程、運(yùn)動(dòng)方程、能量方程、氣體狀態(tài)方程等非線性偏微分方程組來(lái)描述。要求出該方程組的解析解是非常困難，工程應(yīng)用中一般只求得方程組的近似解或數(shù)值解。數(shù)值解法主要是通過(guò)特征線法、隱式中心法等將偏微分方程轉(zhuǎn)化成差分方程求解管內(nèi)氣體流動(dòng)狀態(tài)分布和動(dòng)態(tài)變化過(guò)程。這些方法可以直接服務(wù)于天然氣管道的設(shè)計(jì)計(jì)算、流動(dòng)分析等各個(gè)方面。但是這些方法的計(jì)算效率較低，并不適合狀態(tài)估計(jì)、實(shí)時(shí)優(yōu)化之類的天然氣管道在線應(yīng)用，無(wú)法滿足管道實(shí)時(shí)控制的要求。

1984年，Kralik[1]等提出一種簡(jiǎn)化的天然氣管道傳遞函數(shù)模型，探討了管道控制偏微分方程通過(guò)隱式中心差分空間離散而引入的誤差，研究結(jié)果可以提高復(fù)雜天然氣管道的計(jì)算效率。雖然Kralik僅僅把這種方法用于天然氣管道離線仿真，但是卻為傳遞函數(shù)的在線應(yīng)用提供了新思路。其他學(xué)者提出液體管道傳遞函數(shù)模型并在頻域進(jìn)行分析，隨后傳遞函數(shù)模型被用于水管道的泄漏檢測(cè)。2006年，Reddy[2]等引用Kralik提出的天然氣管道傳遞函數(shù)模型用于分輸量發(fā)生變化的在線狀態(tài)估計(jì)和管道破裂的在線泄漏檢測(cè)。2010年，Behbahani-Nejad[3]等采用二階傳遞函數(shù)在MATLAB-Simulink中構(gòu)建了用于模擬天然氣管道中流體運(yùn)動(dòng)的模型。由于傳遞函數(shù)模型不適用于具有復(fù)雜拓?fù)浣Y(jié)構(gòu)的天然氣管道，而狀態(tài)空間模型可以很好地解決這一問(wèn)題，Behbahani-Nejad[4]等進(jìn)一步通過(guò)傳遞函數(shù)得到天然氣管道的狀態(tài)空間模型。2012年，Alamian[5]等推導(dǎo)了不同邊界條件下天然氣管道傳遞函數(shù)，并且得到不同邊界條件約束的狀態(tài)空間模型。2015年，K.A.Pambour[6]等針對(duì)天然氣管網(wǎng)的仿真控制提出了一種集總參數(shù)瞬態(tài)模型。2015年，R.Whalley[7]等建立了天然氣長(zhǎng)輸管道的分布式參數(shù)模型，通過(guò)魯棒控制算法提出了一種閉環(huán)優(yōu)化調(diào)控策略。同在2015年，S.E.Mohammed[8]等提出了一種基于網(wǎng)絡(luò)分析與仿真的管道狀態(tài)評(píng)價(jià)模型。2016年，劉亞龍[9]等提出了一種基于自適應(yīng)漸進(jìn)約束參量和控制目標(biāo)函數(shù)的控制模型。Wen[10]的研究表明集總參數(shù)狀態(tài)空間模型是對(duì)傳遞函數(shù)模型的改進(jìn)，去掉了控制直通項(xiàng)，使方程的響應(yīng)過(guò)程更接近于實(shí)際管道。

本文采用機(jī)理建模的方法，通過(guò)天然氣管道流動(dòng)的基本方程推導(dǎo)建立狀態(tài)空間模型。一般仿真模型都是從管道的偏微分方程組通過(guò)數(shù)值解法求解，本文提出用狀態(tài)空間模型對(duì)管道水力進(jìn)行建模分析。狀態(tài)空間模型是一種集總參數(shù)模型，在機(jī)理建模的基礎(chǔ)上將系統(tǒng)簡(jiǎn)化為狀態(tài)量、輸入和輸出3類變量。通過(guò)綜合空間變量個(gè)數(shù)，實(shí)現(xiàn)在保證計(jì)算精度的同時(shí)提高計(jì)算效率。狀態(tài)空間模型不僅可用于實(shí)時(shí)控制，還能夠良好的體現(xiàn)管道本身的非線性、大時(shí)滯等特性。對(duì)于仿真模型，時(shí)間上和空間上的劃分精細(xì)程度，通常對(duì)仿真結(jié)果有很大影響。同時(shí)由于在建模過(guò)程中進(jìn)行了一定的近似和線性替換，所以需要對(duì)模型的適用范圍進(jìn)行分析。為展示模型的實(shí)用性，本文選定管徑1016 mm的天然氣管道作為研究對(duì)象。

2 天然氣管道的狀態(tài)空間模型

狀態(tài)空間模型依據(jù)管內(nèi)氣體動(dòng)力學(xué)方程進(jìn)行推導(dǎo)，其中有2處簡(jiǎn)化假設(shè)：(1)輸送過(guò)程中氣體溫度場(chǎng)變化忽略不計(jì)；(2)忽略動(dòng)量方程中對(duì)流項(xiàng)的影響。以下推導(dǎo)以剛性管道的一元流動(dòng)為基礎(chǔ)，針對(duì)控制體的動(dòng)力學(xué)方程，即連續(xù)性方程、運(yùn)動(dòng)方程和能量方程[11]，加上氣體狀態(tài)方程構(gòu)成方程組(1)。

式中：m為氣體的質(zhì)量流量，kg/s；p為氣體的壓力，Pa；A為管道流通截面面積，m2；ρ為氣體的密度，kg/m3；D為管道內(nèi)徑，m；g為重力加速度，m/s2；Δh為管道與水平面的垂直高度，m；L為管道的長(zhǎng)度，m；R為氣體常數(shù)，kJ/(kg·K)；λ為管道水力摩阻系數(shù)；Z為氣體壓縮因子；c為氣體的波速，m/s；x為管道位置變量，m；t為時(shí)間變量，s。

天然氣長(zhǎng)輸管道中氣體運(yùn)行大多在阻力平方區(qū)，摩阻系數(shù)λ可以通過(guò)科爾布魯克公式計(jì)算，也可以用由Hofer(Hofer，1973)[12]得到的適用于阻力平方區(qū)的近似顯式公式(2)計(jì)算；當(dāng)氣體確定，壓縮因子Z是氣體壓力和溫度的函數(shù)，可通過(guò)公式(3)計(jì)算。

式中：Re為輸氣管道的雷諾數(shù)；r為管壁的絕對(duì)當(dāng)量粗糙度，m；pc氣體的臨界壓力，Pa；Tc氣體的臨界溫度，K。

管段的傳遞函數(shù)則是在方程組(1)的基礎(chǔ)上推導(dǎo)得到，是關(guān)于上游的壓力、流量和下游的壓力、流量的頻域表達(dá)式。對(duì)管段的動(dòng)力學(xué)方程線性化并進(jìn)行拉普拉斯變換獲得頻域上的常微分方程組。假定管道的進(jìn)口壓力和出口流量被指定為邊界條件，可獲得高階的傳遞函數(shù)表達(dá)式。但通過(guò)這種方式獲得的傳遞函數(shù)很難求得時(shí)域解析解，因此需要采用泰勒展開得到簡(jiǎn)化的傳遞函數(shù)表達(dá)式。

式中：

狀態(tài)空間模型表達(dá)式如下

3 狀態(tài)空間模型的適用性分析

狀態(tài)空間模型是一個(gè)集總參數(shù)模型，即將一整段管道作為一個(gè)對(duì)象進(jìn)行建模和分析。由于管道的分布參數(shù)特性，當(dāng)我們將管道簡(jiǎn)化為某一點(diǎn)的狀態(tài)集合進(jìn)行分析時(shí)，管道的長(zhǎng)度將對(duì)模型的精度和適用性產(chǎn)生直接影響。根據(jù)狀態(tài)空間模型表達(dá)式(6)、(7)，編程實(shí)現(xiàn)模型參數(shù)的計(jì)算。

管道基礎(chǔ)參數(shù)設(shè)定為：管道管徑為1016 mm，粗糙度為0.04 mm，管道整體絕熱，進(jìn)口規(guī)定壓力為10 MPa，出口規(guī)定流量為131 m3/s。管段按30、40、60、80、100、120、140、160、200 km共9種長(zhǎng)度情況進(jìn)行分析，對(duì)于不同的管道長(zhǎng)度分別代入方程進(jìn)行參數(shù)計(jì)算。

我們首先根據(jù)上述參數(shù)在SPS中搭建對(duì)應(yīng)的仿真模型，進(jìn)口規(guī)定壓力為10 MPa，出口規(guī)定流量為131 m3/s，模型如圖1。

在上文中推導(dǎo)的狀態(tài)空間模型的基礎(chǔ)上，根據(jù)管道對(duì)應(yīng)狀態(tài)參數(shù)和狀態(tài)方程系數(shù)矩陣，計(jì)算得到描述管道系統(tǒng)特性的對(duì)應(yīng)狀態(tài)空間模型的參數(shù)。在MATLAB-Simulink中建立管道的可視化模型，包括描述管段特性的狀態(tài)空間模型模塊以及對(duì)應(yīng)連接的輸入輸出模塊。將計(jì)算得到狀態(tài)空間模型的參數(shù)矩陣輸入到對(duì)應(yīng)的狀態(tài)空間模塊中，狀態(tài)空間模型的輸入是管段的進(jìn)口壓力和出口流量，輸出是進(jìn)口流量和出口壓力。建立的Simulink模型如圖2所示。

在分析中，我們對(duì)兩種輸入變化展開研究：一種是輸入的進(jìn)口壓力產(chǎn)生階躍擾動(dòng)(0.2 MPa)而出口流量不變，而另一種是壓力、流量都產(chǎn)生階躍擾動(dòng)(壓力0.2 MPa、流量4 m3/s)。在這兩種情況下，分別修改管段為30、40、60、80、100、120、140、160、200 km共9種長(zhǎng)度進(jìn)行分析。

圖1 SPS模型Fig. 1 SPS model

圖2 MATLAB-Simulink模型Fig. 2 Matlab-simulink model

3.1 輸入變化為壓力階躍

當(dāng)輸入信號(hào)的變化只有壓力階躍時(shí)，不同管段長(zhǎng)度下的仿真對(duì)比結(jié)果如下所述。

從圖3-1和圖3-2中可以明顯看出由于管段長(zhǎng)度增加導(dǎo)致的穩(wěn)定時(shí)間增加，同時(shí)模型與SPS的仿真結(jié)果之間的偏差也在增大。由仿真結(jié)果可知，當(dāng)輸入出現(xiàn)擾動(dòng)時(shí)，輸出的入口流量和出口壓力值也隨之產(chǎn)生變化。對(duì)于入口流量，無(wú)論管段長(zhǎng)度如何，其輸出流量都會(huì)在擾動(dòng)出現(xiàn)的時(shí)間點(diǎn)瞬時(shí)增加。由于出口流量不變，輸出的流量值也會(huì)逐漸趨于穩(wěn)定值，即設(shè)定的流量值。對(duì)于出口壓力，在擾動(dòng)出現(xiàn)后，輸出的壓力也會(huì)隨之變化，然后達(dá)到新的穩(wěn)態(tài)值。

對(duì)于只有壓力變化的情況，在100 km以上長(zhǎng)度的管道，流量經(jīng)歷短期的波動(dòng)后，最后都會(huì)回到原來(lái)的穩(wěn)定值。這說(shuō)明狀態(tài)空間模型在流量的關(guān)系上可以達(dá)到與SPS模型相同的結(jié)果。圖3-1(a)和圖3-2(a)可以看到狀態(tài)空間模型結(jié)果會(huì)出現(xiàn)振蕩，這是由于管道長(zhǎng)度過(guò)短時(shí)，管道本身的壓力波傳播特性占主要因素，會(huì)表現(xiàn)出一定的波動(dòng)。在實(shí)際管道中由于摩阻和氣體傳熱特性，并不會(huì)出現(xiàn)較大的波動(dòng)。因此，狀態(tài)空間模型并不適用分段長(zhǎng)度過(guò)小的情況。由圖3-1(c)和圖3-2(c)可以看出對(duì)于過(guò)大的管段分段情形，狀態(tài)空間模型的代數(shù)特性使它表現(xiàn)出與SPS模型在反應(yīng)時(shí)間上較大的偏差。圖3-1(b)和圖3-2(b)反映了狀態(tài)空間模型的仿真結(jié)果與SPS的仿真結(jié)果較為符合。其中100 km長(zhǎng)管段條件下振蕩完全消失，且在仿真反應(yīng)時(shí)間偏差較小的情況下與SPS仿真結(jié)果較為接近。因此，對(duì)于管徑為1016的氣體管道，將其切分為100 km左右長(zhǎng)度的管段采用狀態(tài)空間模型進(jìn)行分析比較合理。

3.2 輸入擾動(dòng)為壓力和流量同時(shí)存在

當(dāng)輸入的變化既有壓力階躍，又有流量階躍(4 m3/s)時(shí)，不同管段長(zhǎng)度下的仿真結(jié)果如下所述。

在壓力和流量同時(shí)出現(xiàn)階躍擾動(dòng)的情況下，我們可以得到與單獨(dú)存在壓力擾動(dòng)情況下相同的結(jié)論。

由圖3-1(a)、圖3-2(a)、圖4-1(a)和圖4-2(a)中發(fā)現(xiàn)無(wú)論輸入是否存在流量階躍擾動(dòng)，當(dāng)管道長(zhǎng)度較小時(shí)，模型仿真都會(huì)存在振蕩現(xiàn)象。由3-1(b)、圖3-2(b)、圖4-1(b)和圖4-2(b)可以發(fā)現(xiàn)振蕩現(xiàn)象直到管道長(zhǎng)度增加到80 km時(shí)才消失。理論上，隨著管段的增長(zhǎng)，模型的精度會(huì)逐漸降低。而在本算例中，由于振蕩的存在導(dǎo)致在管段長(zhǎng)度較小時(shí)，精度不是很高；當(dāng)振蕩消失時(shí)由于長(zhǎng)度較大，精度也不是很高。狀態(tài)空間模型是由傳遞函數(shù)模型變形得到的，使用這類模型進(jìn)行仿真時(shí)，當(dāng)管段長(zhǎng)度足夠小時(shí)均會(huì)出現(xiàn)振蕩現(xiàn)象。隨著管段長(zhǎng)度的增加，振蕩則會(huì)消失。管徑越大，振蕩現(xiàn)象消失時(shí)對(duì)應(yīng)的管段長(zhǎng)度越長(zhǎng)。

圖3-1 只有壓力階躍時(shí)，不同管道長(zhǎng)度下的輸出流量對(duì)比Fig. 3-1 Comparison of output flow with different pipeline length

圖3-2 只有壓力階躍時(shí)，不同管道長(zhǎng)度下的輸出壓力對(duì)比Fig. 3-2 Comparison of output pressure with different pipeline length

圖4-1 即有壓力階躍，又有流量階躍時(shí)，不同管道長(zhǎng)度下的輸出流量對(duì)比Fig. 4-1 Comparison of output flow with different pipeline length

由圖3-1(c)、圖3-2(c)、圖4-1(c)和圖4-2(c)可以發(fā)現(xiàn)，管道長(zhǎng)度較長(zhǎng)時(shí)，狀態(tài)空間模型仿真穩(wěn)定時(shí)間較長(zhǎng)且與SPS仿真結(jié)果會(huì)有較大偏差。管段狀態(tài)空間模型在一階條件下其實(shí)是把管段作為存儲(chǔ)可壓縮氣體的圓筒看待。當(dāng)管段較短時(shí)，模型特性無(wú)法得到滿足，則會(huì)出現(xiàn)振蕩現(xiàn)象，而當(dāng)管段長(zhǎng)度較長(zhǎng)時(shí)，一階模型并不能良好的描述系統(tǒng)，會(huì)出現(xiàn)較大的偏差。故而對(duì)于管徑為1016的天然氣管道，狀態(tài)空間模型的適用范圍為80～120 km長(zhǎng)管段，其中以100 km為最佳。

通過(guò)本算例我們發(fā)現(xiàn)狀態(tài)空間模型可以簡(jiǎn)化天然氣管道系統(tǒng)的建模，在不影響仿真精度的同時(shí)，減少了建模所需的管段數(shù)量。由于能夠兼顧計(jì)算效率和精度需求，該模型適合于管道的動(dòng)態(tài)控制。

4 狀態(tài)空間模型魯棒性分析

在由傳遞函數(shù)得到狀態(tài)空間方程的過(guò)程中，我們用穩(wěn)態(tài)線性化得到狀態(tài)空間模型用于系統(tǒng)控制特性的分析。如果模型對(duì)于穩(wěn)態(tài)點(diǎn)的選取過(guò)于敏感，則在此條件下做出的分析也只是針對(duì)特定的工況點(diǎn)而失去實(shí)際意義。因此，我們將分析不同工況點(diǎn)得到狀態(tài)空間模型的差異性，從而說(shuō)明狀態(tài)空間模型是否能準(zhǔn)確描述管道在各個(gè)工況點(diǎn)的動(dòng)態(tài)特性。

在對(duì)連續(xù)性方程和運(yùn)動(dòng)方程進(jìn)行線性化處理時(shí)，不同的穩(wěn)態(tài)可以得到不同的傳遞函數(shù)。對(duì)于長(zhǎng)度為60 km直徑為1016 mm的天然氣管道，在穩(wěn)態(tài)流量 為55 m3/s、70 m3/s、85 m3/s、100 m3/s、115 m3/s、130 m3/s和進(jìn)口壓力為7.6 MPa、8.0 MPa、8.4 MPa、8.8 MPa條件下采用Bode分析展開結(jié)果。

線性管段模型的不確定性主要來(lái)源于一些參數(shù)值取進(jìn)出口的平均近似結(jié)果以及參數(shù)由于運(yùn)行點(diǎn)的改變而發(fā)生變化。圖5-1～5-4對(duì)不同輸入輸出組合下的狀態(tài)空間模型的魯棒性進(jìn)行了分析，流量越小阻尼比越小、諧振峰值的幅值越大，進(jìn)口壓力越大阻尼比越小、諧振峰值的幅值越大，諧振峰值的大小決定著時(shí)間響應(yīng)振蕩的強(qiáng)弱。對(duì)于1016管徑管道，建立的狀態(tài)空間模型在不同流量點(diǎn)線性化得到的幅頻響應(yīng)均沒(méi)有出現(xiàn)較大的偏移，模型自身的穩(wěn)定特性較好。這說(shuō)明狀態(tài)空間模型在描述天然氣管道時(shí)，對(duì)于管道不同工況點(diǎn)具有良好的一致性。

圖5-1 p1、p2組合傳遞函數(shù)bode分析Fig. 5-1 Bode analysis of p1 and p2 combined transfer function

為了驗(yàn)證模型的魯棒性還需要對(duì)模型在單壓力變化和單流量變化下的穩(wěn)定性進(jìn)行分析，其bode圖如圖6-1和圖6-2所示。

圖5-2 m2、p2組合傳遞函數(shù)bode分析Fig. 5-2 Bode analysis of m2 and p2 combined transfer function

圖5-3 p1、m1組合傳遞函數(shù)bode分析Fig. 5-3 Bode analysis of p1 and m1 combined transfer function

圖5-4 m2、m1組合傳遞函數(shù)bode分析Fig. 5-4 Bode analysis of m2 and m1 combined transfer function

圖6-1 流量在標(biāo)稱狀態(tài)附近變化Fig. 6-1 Flow changes under standard state

圖6-1和圖6-2為分別選取穩(wěn)態(tài)流量為100 m3/s和進(jìn)口壓力8.0 MPa為標(biāo)稱狀態(tài)。圖6-1中假定穩(wěn)態(tài)流量不變，穩(wěn)態(tài)進(jìn)口壓力為7.2 MPa、7.6 MPa、8.0 MPa、8.4 MPa、8.8 MPa；而 在 圖6-2中 假 定穩(wěn)態(tài)進(jìn)口壓力不變，穩(wěn)態(tài)流量為70 m3/s、85 m3/s、100 m3/s、115 m3/s、130 m3/s。當(dāng)壓力、流量分別發(fā)生改變時(shí)，系統(tǒng)的頻率特性并無(wú)較大偏移，狀態(tài)空間模型的動(dòng)態(tài)性能仍然非常接近。因此，我們認(rèn)為狀態(tài)空間模型在描述適用長(zhǎng)度的管段時(shí)具有較強(qiáng)的魯棒性。

圖6-2 壓力在標(biāo)稱狀態(tài)附近變化Fig. 6-2 Pressure changes under standard state

5 結(jié)論

狀態(tài)方程作為控制理論中常用的建模方法，能夠?qū)⒐艿老到y(tǒng)簡(jiǎn)化為狀態(tài)量、輸入和輸出3個(gè)模塊，兼顧計(jì)算精度和計(jì)算效率。選取管段長(zhǎng)度是建模時(shí)需要考慮的重要因素，選取管段長(zhǎng)度過(guò)短，會(huì)導(dǎo)致整體模型模塊過(guò)多，建模復(fù)雜；而長(zhǎng)度過(guò)長(zhǎng)時(shí)則會(huì)導(dǎo)致模型失真。本文對(duì)不同管段長(zhǎng)度和工況點(diǎn)的模型準(zhǔn)確性進(jìn)行了分析，結(jié)果可以分為3大類：第1類建模管段長(zhǎng)度過(guò)短，模型振蕩嚴(yán)重；第2類長(zhǎng)度適中，仿真結(jié)果良好；第3類長(zhǎng)度過(guò)長(zhǎng)，模型失真。對(duì)于只關(guān)注管道進(jìn)出口壓力、流量的情況下，大口徑天然氣管道可以用狀態(tài)空間方程分段近似分析。該模型可以將管道的分段擴(kuò)大到100 km左右，這對(duì)于提高計(jì)算速度和進(jìn)一步綜合系統(tǒng)整體信息有很大的作用。狀態(tài)空間模型對(duì)于工況點(diǎn)鄰域內(nèi)的波動(dòng)不是非常敏感，具有很好的動(dòng)態(tài)擴(kuò)展性。狀態(tài)空間模型能夠描述管道本身的非線性、大時(shí)滯等特點(diǎn)，在此基礎(chǔ)上可以采用經(jīng)典控制方法進(jìn)行控制器的設(shè)計(jì)。

(責(zé)任編輯 王雨墨 編輯 馬桂霞)