殷雄，溫凱*，謝萍，閃向營，夏子杰，宮敬*

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

1 天然氣管道流動計算模型介紹

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

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

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

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

狀態(tài)空間模型依據(jù)管內(nèi)氣體動力學(xué)方程進行推導(dǎo)，其中有2處簡化假設(shè)：(1)輸送過程中氣體溫度場變化忽略不計；(2)忽略動量方程中對流項的影響。以下推導(dǎo)以剛性管道的一元流動為基礎(chǔ)，針對控制體的動力學(xué)方程，即連續(xù)性方程、運動方程和能量方程[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為管道的長度，m；R為氣體常數(shù)，kJ/(kg·K)；λ為管道水力摩阻系數(shù)；Z為氣體壓縮因子；c為氣體的波速，m/s；x為管道位置變量，m；t為時間變量，s。

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

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

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

式中：

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

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

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

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

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

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

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

圖1 SPS模型Fig. 1 SPS model

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

3.1 輸入變化為壓力階躍

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

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

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

3.2 輸入擾動為壓力和流量同時存在

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

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

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

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

圖3-2 只有壓力階躍時，不同管道長度下的輸出壓力對比Fig. 3-2 Comparison of output pressure with different pipeline length

圖4-1 即有壓力階躍，又有流量階躍時，不同管道長度下的輸出流量對比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)，管道長度較長時，狀態(tài)空間模型仿真穩(wěn)定時間較長且與SPS仿真結(jié)果會有較大偏差。管段狀態(tài)空間模型在一階條件下其實是把管段作為存儲可壓縮氣體的圓筒看待。當(dāng)管段較短時，模型特性無法得到滿足，則會出現(xiàn)振蕩現(xiàn)象，而當(dāng)管段長度較長時，一階模型并不能良好的描述系統(tǒng)，會出現(xiàn)較大的偏差。故而對于管徑為1016的天然氣管道，狀態(tài)空間模型的適用范圍為80～120 km長管段，其中以100 km為最佳。

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

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

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

在對連續(xù)性方程和運動方程進行線性化處理時，不同的穩(wěn)態(tài)可以得到不同的傳遞函數(shù)。對于長度為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和進口壓力為7.6 MPa、8.0 MPa、8.4 MPa、8.8 MPa條件下采用Bode分析展開結(jié)果。

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

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

為了驗證模型的魯棒性還需要對模型在單壓力變化和單流量變化下的穩(wě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和進口壓力8.0 MPa為標(biāo)稱狀態(tài)。圖6-1中假定穩(wěn)態(tài)流量不變，穩(wěn)態(tài)進口壓力為7.2 MPa、7.6 MPa、8.0 MPa、8.4 MPa、8.8 MPa；而 在 圖6-2中 假 定穩(wěn)態(tài)進口壓力不變，穩(wěn)態(tài)流量為70 m3/s、85 m3/s、100 m3/s、115 m3/s、130 m3/s。當(dāng)壓力、流量分別發(fā)生改變時，系統(tǒng)的頻率特性并無較大偏移，狀態(tài)空間模型的動態(tài)性能仍然非常接近。因此，我們認為狀態(tài)空間模型在描述適用長度的管段時具有較強的魯棒性。

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

5 結(jié)論

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

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