基于瞬態(tài)仿真的天然氣管道動(dòng)態(tài)管存分析

溫凱，焦健豐，袁運(yùn)棟，董楠，殷雄，樊迪，宮敬

1.中國(guó)石油大學(xué)（北京）城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室；2.中油國(guó)際管道有限公司；3.中國(guó)石油工程建設(shè)有限公司

0 引言

隨著“雙碳”目標(biāo)帶動(dòng)的能源結(jié)構(gòu)轉(zhuǎn)型持續(xù)推進(jìn)，天然氣作為最清潔的化石資源，在中國(guó)實(shí)現(xiàn)“碳達(dá)峰、碳中和”任務(wù)目標(biāo)過(guò)程中起到橋梁和支撐作用［1-3］，其在中國(guó)能源消費(fèi)占比從 2010年的 4%，提升至2020年的8.2%［4-5］。國(guó)家石油天然氣管網(wǎng)集團(tuán)有限公司（簡(jiǎn)稱國(guó)家管網(wǎng)集團(tuán)）的成立標(biāo)志著整個(gè)天然氣市場(chǎng)向著上游多主體供應(yīng)、中游單管網(wǎng)集輸、下游多市場(chǎng)競(jìng)爭(zhēng)的“X+1+X”格局加速推進(jìn)［6-7］。隨著天然氣在能源消費(fèi)中的占比日漸提高、交易環(huán)節(jié)和主體的增多［8-10］，輸差作為天然氣管輸過(guò)程中一個(gè)極其重要的指標(biāo)越發(fā)引起重視。輸差與天然氣各環(huán)節(jié)主體的經(jīng)濟(jì)效益直接掛鉤，受多種實(shí)際因素影響。為提高天然氣貿(mào)易的精細(xì)化管理，保障天然氣管網(wǎng)高效運(yùn)行，有必要研究天然氣管道輸送的輸差問(wèn)題。天然氣管道輸差的產(chǎn)生原因分為計(jì)量偏差和計(jì)算不精確兩大類［11-13］。近年來(lái)，研究人員在傳統(tǒng)分析方法的基礎(chǔ)上，借助人工智能、大數(shù)據(jù)分析等新算法深入分析輸差的影響因素，最終，都將問(wèn)題的關(guān)鍵點(diǎn)指向管存計(jì)算［14-19］。

目前，管存的計(jì)算方法主要有理論計(jì)算、軟件仿真及在線數(shù)據(jù)分析等。研究人員基于實(shí)際氣體的狀態(tài)方程，采用理論公式計(jì)算管存的過(guò)程中發(fā)現(xiàn)公式法耗時(shí)費(fèi)力，同時(shí)存在壓縮因子計(jì)算不準(zhǔn)確的問(wèn)題［20-27］。林棋等［28］采用微粒群進(jìn)階算法改進(jìn)了上述問(wèn)題，同時(shí)制訂了在不同輸量下的管存控制原則。林楊等［29］在西氣東輸天然氣管道的管存計(jì)算中設(shè)計(jì)了SCADA（數(shù)據(jù)采集與監(jiān)視控制）數(shù)據(jù)質(zhì)量補(bǔ)救方法，有效提高了計(jì)算準(zhǔn)確度，并使用蘇霍夫公式結(jié)合相關(guān)數(shù)據(jù)，得出實(shí)際的沿線溫度，一定程度上提高了數(shù)據(jù)的可靠性。上述管存計(jì)算均基于穩(wěn)態(tài)算法實(shí)現(xiàn)，但在管輸過(guò)程中天然氣往往處于非穩(wěn)態(tài)。李猷嘉［30］比較了穩(wěn)態(tài)與非穩(wěn)態(tài)兩種計(jì)算方法，結(jié)果表明穩(wěn)態(tài)計(jì)算方法的結(jié)果偏于安全。

自20世紀(jì)90年代以來(lái)，國(guó)外各水力仿真模擬軟件的引入為管存計(jì)算帶來(lái)了幫助。邢曉凱等［31］將TGNET（氣體管道瞬態(tài)和穩(wěn)態(tài)模擬計(jì)算軟件）模擬的管存值同穩(wěn)態(tài)計(jì)算公式進(jìn)行對(duì)比，結(jié)果表明穩(wěn)態(tài)計(jì)算公式有所偏差；但朱瑞華等［32］進(jìn)一步加入氣體熱力學(xué)過(guò)程，并采用蘇霍夫公式進(jìn)行溫度修正，計(jì)算結(jié)果與TGNET最接近。李欣澤［33］使用TGNET仿真軟件計(jì)算出某管道的管存并確定了一些經(jīng)典輸量下管存的合理控制區(qū)間。丁媛等［34］以陜京天然氣管道為例，使用 SPS（石油天然氣長(zhǎng)輸管道模擬計(jì)算軟件）進(jìn)行瞬態(tài)模擬，由于瞬態(tài)模擬相較于穩(wěn)態(tài)模擬考慮到了流量隨時(shí)間的變化，減小了流動(dòng)參數(shù)隨時(shí)間變化的誤差。使用仿真軟件雖然能夠較快地計(jì)算出管存，但若該仿真系統(tǒng)不是在線仿真系統(tǒng)，那么工況的變化將會(huì)對(duì)管存計(jì)算產(chǎn)生較大影響［33］，同時(shí)仿真軟件需要人為設(shè)置管道的基礎(chǔ)數(shù)據(jù)及邊界條件，而往往設(shè)定值并不等于現(xiàn)場(chǎng)所采集的數(shù)據(jù)。

目前，國(guó)內(nèi)外大型管道公司在管存運(yùn)行管理方面也紛紛采用動(dòng)態(tài)管存，如：國(guó)家管網(wǎng)集團(tuán)北京油氣調(diào)控中心將管道生產(chǎn)系統(tǒng)應(yīng)用于實(shí)際運(yùn)營(yíng)中，通過(guò) SCADA系統(tǒng)采集的實(shí)時(shí)數(shù)據(jù)計(jì)算管存以指導(dǎo)管道運(yùn)行管理；哥倫比亞管道集團(tuán)利用 IPS（基于工業(yè)互聯(lián)網(wǎng)的智能化管道系統(tǒng)）將內(nèi)外部數(shù)據(jù)充分整合，在此基礎(chǔ)上搭建了工作管理系統(tǒng)、商業(yè)調(diào)度系統(tǒng)用以實(shí)時(shí)計(jì)算管存；TransCanada（加拿大橫加管道公司）基于建設(shè)運(yùn)營(yíng)期的數(shù)據(jù)及跨業(yè)務(wù)信息的共享，實(shí)現(xiàn)管存計(jì)算及預(yù)測(cè)性分析功能。

由于管存計(jì)算的精確程度與所使用的計(jì)算方法密切相關(guān)，同時(shí)管道輸氣工況在氣源供氣量變化與用戶需求波動(dòng)的影響下常處于非穩(wěn)定狀態(tài)，而目前現(xiàn)場(chǎng)缺少國(guó)產(chǎn)自主研發(fā)的動(dòng)態(tài)管存計(jì)算方法，因此，為適應(yīng)天然氣管道貿(mào)易精細(xì)化管理的要求，本文基于有限容積法開(kāi)發(fā)一套具有自主知識(shí)產(chǎn)權(quán)的仿真算法，并以此為基礎(chǔ)，進(jìn)行天然氣管道非穩(wěn)定流動(dòng)工況下的管存分析。在分析天然氣管道穩(wěn)態(tài)計(jì)算方法及常用的AGA8氣體狀態(tài)方程的基礎(chǔ)上，推導(dǎo)了基于有限容積法離散后的天然氣管道瞬態(tài)仿真模型，從理論上分析了穩(wěn)態(tài)管存與動(dòng)態(tài)管存計(jì)算方法的差異；使用現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)驗(yàn)證了仿真模型的準(zhǔn)確性，并對(duì)比分析了不同工況下管道穩(wěn)態(tài)與動(dòng)態(tài)管存差異，以指導(dǎo)管道運(yùn)行方進(jìn)行管存計(jì)算與輸差分析。

1 穩(wěn)態(tài)管存計(jì)算

1.1 穩(wěn)態(tài)管存

管道某一時(shí)刻實(shí)際容納的氣體數(shù)量稱為管存，通常用工程標(biāo)準(zhǔn)狀況下的體積表示。在計(jì)算上有兩種方法：穩(wěn)態(tài)管存和動(dòng)態(tài)管存。穩(wěn)態(tài)管存即認(rèn)為管道運(yùn)行工況為穩(wěn)定工況，根據(jù)QSY 05197—2019《油氣管道輸送損耗計(jì)算方法》中的天然氣管道管存量計(jì)算方法，同時(shí)參考某國(guó)際天然氣管道實(shí)際運(yùn)行情況，可以用平均壓力和平均溫度表征管內(nèi)氣體當(dāng)前狀態(tài)，并以此為依據(jù)計(jì)算管存，具體計(jì)算方法如式（1）和式（2）。

式中：V—管道管存量，m3；p1—管道入口壓力，Pa；p2—管道出口壓力，Pa；pav—管道平均壓力，Pa；V1—管道幾何容積，m3；Z1—管道平均壓力、平均溫度下的氣體壓縮因子；T1—管道平均溫度，K；Z0—工程標(biāo)準(zhǔn)狀況下的氣體壓縮因子；T0—工程標(biāo)準(zhǔn)狀況下的溫度，取293.15 K；p0—工程標(biāo)準(zhǔn)狀況下的壓力，取101.3 kPa。

1.2 物性計(jì)算

天然氣在不同壓力及溫度下的密度計(jì)算精度直接影響到天然氣管道水力、熱力計(jì)算及管存計(jì)算的準(zhǔn)確性，因此確定合適的物性計(jì)算方法至關(guān)重要。參照GB/T 17747.2—2011《天然氣壓縮因子的計(jì)算第2部分：用摩爾組成進(jìn)行計(jì)算》的推薦算法，同時(shí)，為了與某國(guó)際天然氣管道（下文算例）物性計(jì)算方法保持一致，本文采用AGA8—92DC方程（簡(jiǎn)稱AGA8方程）進(jìn)行物性計(jì)算，在輸入的參數(shù)誤差滿足標(biāo)準(zhǔn)且氣體的熱力狀態(tài)在公式應(yīng)用范圍之內(nèi)時(shí)，AGA8方程的不確定度小于 0.1%。AGA8方程于 1992年發(fā)表在 AGA（美國(guó)煤氣協(xié)會(huì)）的 AGA8報(bào)告上，主要適用于絕對(duì)壓力 0～12 MPa、溫度263～338 K的管輸天然氣。該方程可輸入多達(dá) 21種組分的摩爾比例進(jìn)行計(jì)算，同時(shí)要求組分的摩爾分?jǐn)?shù)大于0.000 05；對(duì)于H2S和CO2，適用的摩爾分?jǐn)?shù)范圍分別為0～0.02%，0～30%。

2 基于瞬態(tài)仿真的動(dòng)態(tài)管存計(jì)算

2.1 動(dòng)態(tài)管存

在天然氣管道的運(yùn)行中，常常伴隨著運(yùn)行工況的調(diào)整，如管道的進(jìn)氣量發(fā)生變化、壓縮機(jī)的啟停、調(diào)節(jié)閥的動(dòng)作、沿線的下載氣量發(fā)生變化等，這些都會(huì)導(dǎo)致管道處于非穩(wěn)態(tài)，在這種情況下沿線的壓力分布與穩(wěn)態(tài)時(shí)的并不一致，因而不能采用平均壓力公式（1）計(jì)算管道的壓力。動(dòng)態(tài)管存計(jì)算則是根據(jù)管道的基礎(chǔ)數(shù)據(jù)，以管道兩端的實(shí)測(cè)壓力、溫度和流量為邊界條件，實(shí)時(shí)計(jì)算出沿線各個(gè)位置的壓力、溫度、流速和密度等參數(shù)，再根據(jù)管存計(jì)算公式（3），計(jì)算出管道的管存。管道各個(gè)位置的壓力、溫度、流速和密度等氣體狀態(tài)參數(shù)需要通過(guò)流動(dòng)仿真獲得。

式中：ρ—管段實(shí)際壓力、溫度下天然氣的密度，kg/m3；0ρ—工程標(biāo)準(zhǔn)狀況下天然氣的密度，kg/m3；Vd—管段的設(shè)計(jì)管容量，m3；Vs—管段在標(biāo)準(zhǔn)參比條件下的管存量，m3。

2.2 基本模型

天然氣管道流動(dòng)仿真就是采用離散手段將無(wú)法直接得到解析解的描述天然氣管道流動(dòng)的偏微分方程組轉(zhuǎn)化為一系列離散代數(shù)方程組，借助計(jì)算機(jī)完成數(shù)值求解，再通過(guò)對(duì)數(shù)值結(jié)果的分析再現(xiàn)管道內(nèi)天然氣的不穩(wěn)定流動(dòng)狀態(tài)。用來(lái)描述天然氣在管道內(nèi)流動(dòng)狀態(tài)的控制方程稱為管道的數(shù)學(xué)模型，是由連續(xù)性方程、動(dòng)量方程和能量方程構(gòu)成守恒型的天然氣管流的基本控制方程。其中前兩者可稱為水力方程，后者稱為熱力方程。

式中：A—管道橫截面積，m2；t—時(shí)間，s；w—?dú)怏w流速，m/s；x—管道軸向長(zhǎng)度，m；p—管道中氣體的壓力，Pa；D—管道內(nèi)徑，m；θ—管道傾角，°；λ—摩阻系數(shù)，無(wú)量綱；g—重力加速度，m/s2；e—?dú)怏w的內(nèi)能，J/kg；h—?dú)怏w的焓值，J/（kg·K）；s—管道高程，m；Qq—天然氣與環(huán)境之間的換熱量，W/m3。

2.3 基于有限容積法的方程離散

天然氣管道瞬態(tài)仿真的主流數(shù)值算法主要有特征線法、有限差分法、有限元法、狀態(tài)空間模型法及有限容積法等。

特征線法其特點(diǎn)為顯式求解，同時(shí)穩(wěn)定性、適應(yīng)性較高；但求解步長(zhǎng)受條件限制，在快瞬變流研究中較為廣泛［35］?，F(xiàn)較常用的有限差分法為隱式有限差分法，該方法雖然求解速度很慢，但求解步長(zhǎng)不受條件限制，同時(shí)不受物性及狀態(tài)的影響，在長(zhǎng)時(shí)間瞬變流的求解中應(yīng)用較為廣泛［36］。有限元法相較于有限差分法計(jì)算精度有所提高，同時(shí)求解速度也較快，在求解大型復(fù)雜管網(wǎng)時(shí)也具有一定優(yōu)勢(shì)［37］。狀態(tài)空間模型法計(jì)算速度相較于有限差分法大幅提高，這是由于該方法基于傳遞函數(shù)直接表征了管道出入口的流動(dòng)關(guān)系，但對(duì)于管道瞬時(shí)沿線壓力流量難以求解［38］。有限容積法因其離散格式具有穩(wěn)定、精準(zhǔn)等特點(diǎn)，同時(shí)求解步長(zhǎng)不受限制，因此在快、慢瞬變流動(dòng)方面都有廣泛的應(yīng)用［39］。

基于以上對(duì)天然氣瞬態(tài)仿真方法的調(diào)研，本文采用交錯(cuò)網(wǎng)格下的有限容積法進(jìn)行流動(dòng)仿真。有限容積法的基本思想是將問(wèn)題的求解域離散化成一系列小區(qū)域，對(duì)每個(gè)小區(qū)域進(jìn)行積分，得到關(guān)于未知函數(shù)的代數(shù)方程，聯(lián)立求解這些代數(shù)方程就可得到未知函數(shù)在節(jié)點(diǎn)的參數(shù)值。

為了方便積分，需要對(duì)能量方程（6）進(jìn)行轉(zhuǎn)換。根據(jù)熱力學(xué)定律有如下兩式：

式中：u—?dú)怏w黏度，Pa·s；T—管道內(nèi)氣體溫度，K；Cp—?dú)怏w的定壓比熱容J/（kg·K）。

將式（4）、（5）、（7）和（8）帶入式（6），整理可得：

式（4）、式（5）和式（9）組成了天然氣管道仿真的控制方程組。在交錯(cuò)網(wǎng)格下采用有限容積法對(duì)控制方程組進(jìn)行離散（見(jiàn)圖1）。主變量中速度變量位于網(wǎng)格的邊緣，用虛直線標(biāo)記；壓力和溫度變量位于網(wǎng)格的中心，用黑點(diǎn)標(biāo)記。同時(shí)，密度、熱容、聲速等物性參數(shù)也位于網(wǎng)格的中心。

圖1 管道網(wǎng)格劃分示意圖

對(duì)于計(jì)算所需，但該位置并未定義的變量，采用一階迎風(fēng)格式—公式（10）進(jìn)行確定，如網(wǎng)格邊緣的壓力變量或者網(wǎng)格中心的速度變量。式（10）中，i表示第i個(gè)網(wǎng)格節(jié)點(diǎn)。

采用上述變量定義方法，動(dòng)量方程式（5）可在動(dòng)量控制體中離散成公式（11），式（11）中：上標(biāo)“#”代表變量的當(dāng)前值，“*”代表該變量下一次迭代中待求解的值。同理，連續(xù)性方程和能量方程可以在質(zhì)量控制體中分別離散成公式（12）、公式（13）。

該節(jié)中，采用SIMPLE（壓力耦合方程組的半隱式方法）求解天然氣管道瞬態(tài)仿真中的水力、熱力狀態(tài)［40］。將下一時(shí)步速度、壓力和密度的值寫(xiě)成當(dāng)前值加修正量的形式，如式（14）～式（16）。

式中：w′，p′和ρ′分別為速度、壓力和密度的修正量。將式（14）～式（16）帶入式（11），并忽略二階修正項(xiàng)，可得壓力修正方程：

該方法的整體仿真流程如圖2所示。仿真開(kāi)始時(shí)，根據(jù)穩(wěn)態(tài)計(jì)算結(jié)果對(duì)網(wǎng)格上各點(diǎn)參數(shù)賦值。接著根據(jù)離散后的動(dòng)量方程式（11）求得一個(gè)新的速度值w。這個(gè)速度值不能匹配當(dāng)前的壓力值，因此根據(jù)式（17）求解壓力修正量p′。繼而根據(jù)式（16）求解密度修正量ρ′。至此，壓力、速度和密度均進(jìn)行了更新。這個(gè)循環(huán)修正的步驟迭代到壓力修正量足夠小才會(huì)結(jié)束。之后根據(jù)式（13）求解各網(wǎng)格點(diǎn)的溫度值，這樣這一時(shí)步的各點(diǎn)的壓力、流量和溫度均可獲得，接著求解下一步，直到整個(gè)仿真計(jì)算完成。

圖2 交錯(cuò)網(wǎng)格下的有限容積法計(jì)算過(guò)程

無(wú)論是穩(wěn)態(tài)管存算法還是動(dòng)態(tài)算法，其計(jì)算過(guò)程都是先設(shè)法獲得管段中天然氣的壓力和溫度，然后再根據(jù)該壓力和溫度進(jìn)行其他計(jì)算，因此，確定壓力和溫度這一步驟十分關(guān)鍵。穩(wěn)態(tài)算法利用天然氣管道處于穩(wěn)定流動(dòng)時(shí)壓力、溫度的分布特點(diǎn)，給出了一種基于管道兩端的壓力、溫度計(jì)算平均壓力、溫度的方法；動(dòng)態(tài)管存算法則摒棄了這種思路，將管段劃分為長(zhǎng)度為1 km左右的小管段，再根據(jù)管道兩段實(shí)測(cè)的壓力、溫度，仿真出每一小管段的壓力、溫度，最后再進(jìn)行管存計(jì)算。顯然后者管道劃分更加精細(xì)，且對(duì)于非穩(wěn)定工況下的管存描述更加準(zhǔn)確。

3 應(yīng)用實(shí)例

以某國(guó)際天然氣管道 C1至 C2站間管道為例（管道參數(shù)如表1所示，氣體組成如表2所示）。通過(guò)對(duì)比SCADA系統(tǒng)測(cè)量數(shù)據(jù)與仿真數(shù)據(jù)，驗(yàn)證基于有限容積法的天然氣管網(wǎng)仿真方法的準(zhǔn)確性。同時(shí)，針對(duì)此管道進(jìn)行不同壓力波動(dòng)下的管存分析，用于分析管內(nèi)流動(dòng)的不穩(wěn)定性對(duì)于管存計(jì)算的影響。

表1 管道參數(shù)

表2 氣體組成

3.1 模型驗(yàn)證

通過(guò)SCADA系統(tǒng)獲取該管道3 d的運(yùn)行數(shù)據(jù)，包括C1站出站壓力、流量與溫度，以及C2站進(jìn)站壓力與溫度。將C1站出站壓力、C2站進(jìn)站壓力，C1站出站溫度分別作為仿真的管道進(jìn)口壓力、出口壓力及入口溫度，仿真距離步長(zhǎng)為1 km，時(shí)間步長(zhǎng)為60 s，仿真160 h管內(nèi)流動(dòng)情況（見(jiàn)圖3），共耗時(shí)3 min 44 s。將仿真所得C1站出站流量即管道入口流量與 C1站實(shí)測(cè)流量進(jìn)行對(duì)比，平均相對(duì)誤差為1.85%；同理，將仿真所得C2站進(jìn)站溫度即管道出口溫度與 C2站實(shí)測(cè)溫度進(jìn)行對(duì)比，平均相對(duì)誤差為0.92%（見(jiàn)圖4）。以上結(jié)果驗(yàn)證了基于有限容積法的天然氣管網(wǎng)仿真方法計(jì)算速度快且精度高，可以滿足工業(yè)需求。

圖3 仿真水力、熱力邊界

圖4 仿真結(jié)果對(duì)比

3.2 動(dòng)態(tài)管存分析

天然氣管道內(nèi)流動(dòng)的不穩(wěn)定性可以用壓力變化程度進(jìn)行表征。為了定量分析壓力變化程度對(duì)于管存計(jì)算的影響，本文采用管道出口為定壓，入口壓力呈階躍變化或周期變化兩類算例進(jìn)行仿真實(shí)驗(yàn)。對(duì)比入口壓力波動(dòng)，出口壓力保持定值下，不同階躍時(shí)長(zhǎng)或不同周期下的動(dòng)態(tài)管存與穩(wěn)態(tài)管存差異。初始狀態(tài)為入口壓力8.5 MPa，入口溫度45 ℃，出口壓力6.5 MPa的穩(wěn)定狀態(tài)。

3.2.1 入口壓力階躍變化

壓力階躍變化即壓力從一個(gè)穩(wěn)定值線性增長(zhǎng)到另一個(gè)穩(wěn)定值，本算例中其數(shù)學(xué)模型為：

式中：Astep—階躍振幅，MPa；Tstep—階躍時(shí)長(zhǎng)，s。Astep與Tstep共同決定了階躍變化的形狀，因此各設(shè)計(jì)3種階躍幅度與時(shí)長(zhǎng)，組合9種階躍變化（見(jiàn)圖5～圖7）進(jìn)行仿真實(shí)驗(yàn)，每次仿真時(shí)間為1 000 min，得到9種入口壓力階躍變化情況下的穩(wěn)態(tài)與動(dòng)態(tài)管存計(jì)算結(jié)果（見(jiàn)圖8～圖9）。

圖5 入口壓力階躍變化下（階躍時(shí)長(zhǎng)6 h）穩(wěn)態(tài)管存與動(dòng)態(tài)管存對(duì)比

圖6 入口壓力階躍變化下（階躍時(shí)長(zhǎng)2 h）穩(wěn)態(tài)管存與動(dòng)態(tài)管存對(duì)比

圖7 入口壓力階躍變化下（階躍時(shí)長(zhǎng)0.5 h）穩(wěn)態(tài)管存與動(dòng)態(tài)管存對(duì)比

圖8 入口壓力階躍變化下管存偏差對(duì)比

圖9 入口壓力階躍變化下管存相對(duì)管存差對(duì)比

由圖5可以看出，兩類管存均是先上升，后平穩(wěn)，總體趨勢(shì)一致，但是穩(wěn)態(tài)管存上升速度比動(dòng)態(tài)管存快，且流動(dòng)穩(wěn)定后，穩(wěn)態(tài)管存均較動(dòng)態(tài)管存低。同時(shí)，穩(wěn)態(tài)管存變化過(guò)程在t=Tstep時(shí)刻存在一個(gè)明顯的拐點(diǎn)，該拐點(diǎn)之前由于入口壓力快速增長(zhǎng)，因此平均壓力也快速增長(zhǎng)，穩(wěn)態(tài)管存也變化迅速；拐點(diǎn)之后，進(jìn)出口壓力均不變，但溫度場(chǎng)并未達(dá)到穩(wěn)定，出口溫度緩慢下降，造成平均溫度減小，穩(wěn)態(tài)管存緩慢增長(zhǎng)直至穩(wěn)定。由此可知，由于穩(wěn)態(tài)管存計(jì)算方法僅僅與進(jìn)出口兩個(gè)點(diǎn)壓力與溫度有關(guān)，因此階躍變化造成的入口壓力變化會(huì)立刻影響到平均壓力的計(jì)算，而實(shí)際情況是入口壓力的上升是管道充裝氣體造成的，而數(shù)百千米的管道充裝是需要時(shí)間的，因此實(shí)際管存變化應(yīng)該如動(dòng)態(tài)管存所示，平滑變化，且當(dāng)入口壓力階躍結(jié)束后，仍維持一段時(shí)間的管存上升狀態(tài)。此外，流動(dòng)穩(wěn)定后，高壓狀態(tài)下的穩(wěn)態(tài)管存與動(dòng)態(tài)管存偏差變小，說(shuō)明高壓穩(wěn)定狀態(tài)下的平均法計(jì)算管存是相對(duì)有效的；相同階躍振幅下，階躍時(shí)長(zhǎng)越短，管存上升階段（0～Tstep）穩(wěn)態(tài)管存變化斜率較動(dòng)態(tài)管存變化斜率更大，說(shuō)明在管內(nèi)不穩(wěn)定流動(dòng)剛發(fā)生時(shí)段進(jìn)行管存穩(wěn)態(tài)計(jì)算會(huì)造成較大的誤差，且振幅越大，誤差越大。

在輸差分析過(guò)程中，管存所計(jì)算的往往是相對(duì)量，即，終了狀態(tài)管內(nèi)管存量與初始時(shí)刻之差。因此本文對(duì)比了穩(wěn)態(tài)與動(dòng)態(tài)管存差之間的偏差，即：

式中：ΔLP—穩(wěn)態(tài)管存差與動(dòng)態(tài)管存差之間的偏差（簡(jiǎn)稱管存偏差），104m3；dLsP—當(dāng)前時(shí)刻穩(wěn)態(tài)管存與初始時(shí)刻偏差，104m3；dLPd—當(dāng)前時(shí)刻動(dòng)態(tài)管存與初始時(shí)刻偏差，104m3。

由圖6（a）可以看出，管存偏差變化趨勢(shì)均為先增加后減小再平穩(wěn)，且增加減小的速率均是逐漸變緩，拐點(diǎn)出現(xiàn)在t=Tstep時(shí)刻。分析原因可知，0～Tstep時(shí)間段內(nèi)，上升速率大于動(dòng)態(tài)管存，但穩(wěn)態(tài)管存接近線性上升，而動(dòng)態(tài)管存呈拋物線增長(zhǎng)，因此管存偏差呈拋物線增長(zhǎng)；當(dāng)t＞Tstep，穩(wěn)態(tài)管存上升速度迅速下降，低于動(dòng)態(tài)管存，因此管存偏差又下降，直到流動(dòng)穩(wěn)定后平穩(wěn)。相同階躍時(shí)長(zhǎng)下，振幅越大，管存偏差越大；相同階躍振幅下，時(shí)長(zhǎng)越長(zhǎng)，最大管存偏差越小，但流動(dòng)穩(wěn)定后管存偏差相同。

以上為定性分析，以下從定量的角度分析管存計(jì)算方法對(duì)于輸差計(jì)算的影響程度。參考相對(duì)輸差計(jì)算公式［14］，引入相對(duì)管存差參數(shù)RLP，RLP越大，由管存計(jì)算不準(zhǔn)導(dǎo)致的輸差越大，反之亦然。

式中：RLP—相對(duì)管存差，無(wú)量綱；Qall—0～t時(shí)間段內(nèi)的總進(jìn)氣量，104m3；L—輸氣管道總長(zhǎng)度，km；l—某段管道長(zhǎng)度，km。算例中l(wèi)=202.8 km，L=524.4 km，因此RLP=2.59 △RLP/Qall，計(jì)算得到9種階躍變化下的RLP（見(jiàn)圖6（b））。從理論上分析，隨時(shí)間跨度增大，輸氣量增大，由于管存計(jì)算偏差帶來(lái)的輸差會(huì)逐漸減小，這與圖6（b）相符合。除最初幾分鐘外，RLP隨著時(shí)間增加而降低。本文以3%為閾值，當(dāng)RLP＜3%時(shí)，認(rèn)為管存偏差對(duì)輸差影響可以忽略不計(jì)，統(tǒng)計(jì)閾值時(shí)間點(diǎn)Ts，即t＞Ts，RLP＜3%，可以得到入口壓力階躍變化下閾值時(shí)間點(diǎn)（見(jiàn)表3）。此外，將閾值時(shí)間點(diǎn)與階躍時(shí)間相減可得額外穩(wěn)定時(shí)間，用于描述階躍變化停止后到管存偏差對(duì)輸出影響可以忽略不計(jì)的時(shí)間。顯然，階躍時(shí)長(zhǎng)越短波動(dòng)越快，所需額外穩(wěn)定時(shí)間也越長(zhǎng)。如，在階躍振幅都為0.5 MPa時(shí)，階躍時(shí)長(zhǎng)為6 h時(shí)額外穩(wěn)定時(shí)間最短。階躍時(shí)長(zhǎng)為 0.5 h時(shí)穩(wěn)定時(shí)間最長(zhǎng)同時(shí)，階躍振幅越大，所需額外穩(wěn)定時(shí)間也越長(zhǎng)。如，在階躍時(shí)長(zhǎng)都為6 h時(shí)，額外穩(wěn)定時(shí)間隨著振幅的增加而延長(zhǎng)。

表3 入口壓力階躍變化下閾值時(shí)間點(diǎn)

3.2.2 入口壓力周期變化

實(shí)際天然氣管道運(yùn)行工況受上下游供需影響，而供需受溫度等周期變化因素影響而存在一定周期性。因此本文設(shè)計(jì)正弦周期變化的算例，用于分析不同周期與振幅的正弦入口壓力變化下，穩(wěn)態(tài)管存與動(dòng)態(tài)管存計(jì)算結(jié)果偏差，其數(shù)學(xué)模型如下式：

式中：Asin—正弦振幅，MPa；Tsin—正弦周期，s。Asin與Tsin共同決定了正弦變化的形狀，因此各設(shè)計(jì)3種正弦振幅與周期，共組合9種正弦變化（見(jiàn)圖10～圖12），進(jìn)行仿真實(shí)驗(yàn)。

圖10 入口壓力周期變化下（正弦周期24 h）穩(wěn)態(tài)管存與動(dòng)態(tài)管存對(duì)比

圖11 入口壓力周期變化下（正弦周期12 h）穩(wěn)態(tài)管存與動(dòng)態(tài)管存對(duì)比

圖12 入口壓力周期變化下（正弦周期2 h）穩(wěn)態(tài)管存與動(dòng)態(tài)管存對(duì)比

總體來(lái)看，兩類管存變化趨勢(shì)與入口壓力變化相似，均呈現(xiàn)正弦變化，但是動(dòng)態(tài)管存計(jì)算結(jié)果相較于穩(wěn)態(tài)管存存在一定程度的相位差與振幅差。相同Asin、不同Tsin下，穩(wěn)態(tài)管存的正弦變化振幅基本不變，周期也與Tsin保持一致，這是不合理的。在入口壓力持續(xù)變化情況下，管內(nèi)始終處于不穩(wěn)定流動(dòng)狀態(tài)，入口壓力波動(dòng)傳遞到200 km管長(zhǎng)的各處是需要時(shí)間的，造成的管存變化必然滯后。動(dòng)態(tài)管存計(jì)算方法則能較好體現(xiàn)這一特性，管存的變化較之入口壓力變化存在一定程度的時(shí)間滯后。此外，隨著Tsin變小，入口壓力波動(dòng)變大，管內(nèi)氣體無(wú)法充裝到穩(wěn)態(tài)管存量的峰值就下降，因此動(dòng)態(tài)管存的峰值隨著Tsin變小而變小，這與現(xiàn)場(chǎng)實(shí)際情況相符。

類似于階躍變化算例，周期變化算例也對(duì)穩(wěn)態(tài)與動(dòng)態(tài)管存差之間的偏差進(jìn)行了分析（見(jiàn)圖 13）。入口壓力變化振幅越大、周期越短，管存偏差越大，也就是管內(nèi)流動(dòng)波動(dòng)越大，穩(wěn)態(tài)管存計(jì)算方式越不適用，這與直觀理解是相符合的。此外，關(guān)于RLP的結(jié)果（見(jiàn)圖14、表4）所示，其結(jié)果與階躍變化下的結(jié)果略有不同。周期變化下的RLP隨時(shí)間圍繞 0振蕩衰減，且短周期高振幅的周期波動(dòng)所需要的閾值時(shí)間點(diǎn)很長(zhǎng)，如，周期2 h、振幅1.5 MPa算例中，閾值時(shí)間點(diǎn)為4 232 min，需要接近3 d的時(shí)間才可以忽略管存計(jì)算對(duì)于輸差計(jì)算的影響。

表4 入口壓力周期變化下閾值時(shí)間點(diǎn)

圖13 入口壓力周期變化下管存偏差對(duì)比

圖14 入口壓力周期變化下相對(duì)管存差對(duì)比

3.3 小結(jié)

通過(guò)與現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)對(duì)比，驗(yàn)證了基于有限容積法在天然氣管網(wǎng)仿真算法方面的可靠性，同時(shí)通過(guò)對(duì)于入口壓力階躍變化與周期變化兩種情況的穩(wěn)態(tài)管存與動(dòng)態(tài)管存計(jì)算方法對(duì)比可知，動(dòng)態(tài)管存計(jì)算方法的適用性與精確性遠(yuǎn)高于穩(wěn)態(tài)管存計(jì)算方法。對(duì)于管道運(yùn)行方來(lái)說(shuō)，應(yīng)改用動(dòng)態(tài)管存計(jì)算方式進(jìn)行管存管理與輸差分析，以減少由于管存偏差而導(dǎo)致的經(jīng)濟(jì)損失。同時(shí)以下3點(diǎn)還需進(jìn)一步優(yōu)化。

（1）現(xiàn)場(chǎng)數(shù)據(jù)存在大量噪聲與異常，若直接進(jìn)行計(jì)算會(huì)產(chǎn)生計(jì)算偏差，首先應(yīng)進(jìn)行數(shù)據(jù)清洗，通過(guò)提高數(shù)據(jù)可靠性有效降低計(jì)算偏差。

（2）需進(jìn)一步研究SCADA系統(tǒng)和計(jì)算程序之間的通信問(wèn)題，這樣在保證計(jì)算精確度的同時(shí)也提高了運(yùn)行管理的效率。

（3）計(jì)算模型的精確程度直接影響著現(xiàn)場(chǎng)運(yùn)行管理，因此還需要基于運(yùn)行數(shù)據(jù)不斷地修正計(jì)算模型，以提高模型的自適應(yīng)性來(lái)感知現(xiàn)場(chǎng)工況的變化。

4 結(jié)論

精細(xì)化管理是天然氣管道運(yùn)行的重要要求，管存計(jì)算也逐漸成為各管道公司日常運(yùn)行管理內(nèi)容。管存計(jì)算的精確程度不僅涉及輸差問(wèn)題，同時(shí)也會(huì)對(duì)于管道日常運(yùn)行、儲(chǔ)氣調(diào)峰及應(yīng)急供氣產(chǎn)生影響。本文從天然氣管道瞬態(tài)流動(dòng)仿真出發(fā)，提出一套完整的動(dòng)態(tài)管存計(jì)算方法。該方法以有限容積法為核心，求解流動(dòng)三大基本方程，并結(jié)合AGA8氣體狀態(tài)方程，精確計(jì)算天然氣管道非穩(wěn)定流動(dòng)狀態(tài)下沿線壓力溫度密度分布，進(jìn)而獲得管內(nèi)管存精確變化情況。本文用現(xiàn)場(chǎng)數(shù)據(jù)對(duì)該仿真方法進(jìn)行了驗(yàn)證，并對(duì)比了現(xiàn)行穩(wěn)態(tài)管存與本文提出的動(dòng)態(tài)管存計(jì)算方法在不同階躍、周期變化下的偏差，體現(xiàn)了動(dòng)態(tài)管存計(jì)算方法的優(yōu)勢(shì)。此外，隨著天然氣管道運(yùn)行數(shù)字化、智能化，可以預(yù)見(jiàn)將動(dòng)態(tài)管存計(jì)算方法與在線仿真相結(jié)合，可以獲得實(shí)時(shí)管內(nèi)管存變化情況，為運(yùn)行方提供更精確的決策支持信息。