基于隱式法的天然氣管道動(dòng)態(tài)仿真

董 斌

（西安石油大學(xué)，陜西 西安 710065）

伴隨著天然氣生產(chǎn)以及應(yīng)用規(guī)模的不斷擴(kuò)大，輸送天然氣的管網(wǎng)系統(tǒng)也變得越來(lái)越龐大和復(fù)雜。在天然氣管道系統(tǒng)運(yùn)行的過(guò)程中，由于供氣和用氣量的變化，閥門的開(kāi)關(guān)，壓縮機(jī)的啟停，都會(huì)使管網(wǎng)系統(tǒng)處于不穩(wěn)定的運(yùn)行狀態(tài)，為了便于合理的生產(chǎn)和輸送，實(shí)現(xiàn)有效的調(diào)度和管理，對(duì)于管網(wǎng)系統(tǒng)的工況模擬是很有必要的。

天然氣在管網(wǎng)系統(tǒng)中不穩(wěn)定流動(dòng)的問(wèn)題可以采用一些數(shù)值解法并借助計(jì)算機(jī)來(lái)進(jìn)行有效的分析，目前最廣泛采用的數(shù)值解法有特征線法、隱式法[2]等。隱式法的優(yōu)點(diǎn)[1]是它可以采用較長(zhǎng)的時(shí)間步長(zhǎng)，同時(shí)還能夠保證計(jì)算結(jié)果的精確性和穩(wěn)定性，這就可以減少相應(yīng)的計(jì)算次數(shù)。本文主要研究基于隱式法的單根天然氣管道動(dòng)態(tài)仿真，通過(guò)Python語(yǔ)言進(jìn)行計(jì)算機(jī)編程，將其編程結(jié)果與PNS仿真軟件的結(jié)果進(jìn)行對(duì)比，判斷所編程序的可行性，從而為天然氣管道在生產(chǎn)運(yùn)行中提供可靠的數(shù)據(jù)，以便進(jìn)行合理調(diào)度和管理。

1 數(shù)學(xué)微分方程

對(duì)于輸氣管道來(lái)說(shuō)，在不用考慮溫度的變化情況下，流動(dòng)狀態(tài)可以用連續(xù)性方程和運(yùn)動(dòng)方程來(lái)描述。

連續(xù)性方程：

運(yùn)動(dòng)方程：

引入質(zhì)量流量M=ρvA，由于運(yùn)動(dòng)方程中慣性項(xiàng)影響較小，可以忽略，則數(shù)學(xué)微分方程如下：

由熱力學(xué)知?dú)怏w絕熱傳播音速為：

對(duì)等溫流動(dòng)，考慮狀態(tài)方程：

則式 （3）、 （4） 變?yōu)椋?/p>

在方程（7）、（8） 中有兩個(gè)變量P和M，可以聯(lián)立進(jìn)行求解。

2 隱式法

2.1 差分方程

如圖1為隱式法的xt圖。橫坐標(biāo)表示為管道長(zhǎng)度，根據(jù)需要可以劃分為步長(zhǎng)為Δx的網(wǎng)格。縱坐標(biāo)表示時(shí)間，在計(jì)算過(guò)程中也可以離散成時(shí)長(zhǎng)為Δt的網(wǎng)格。對(duì)于每一個(gè)網(wǎng)格節(jié)點(diǎn)來(lái)說(shuō)，可以用兩個(gè)數(shù)去描述它，其中第一個(gè)數(shù)為管段的網(wǎng)格數(shù)，第二個(gè)數(shù)為時(shí)步數(shù)。在圖1中的四個(gè)網(wǎng)格節(jié)點(diǎn)上，未知變量包括 Pi，j+1，Pi+1，j+1，Mi，j+1，Mi+1，j+1。

圖1 隱式法網(wǎng)格圖

而已知變量包括 Pi，j，Pi+1，j，Mi，j，Mi+1，j。

中心有限差分方法是把不穩(wěn)定流的方程以有限差形式應(yīng)用在位居四點(diǎn)網(wǎng)格中心的點(diǎn)子上。因此，方程（7）、（8） 的差分格式可以寫為：

對(duì)于管道系統(tǒng)中的每一個(gè)管段單元，可以寫出類似式（9） 和（10） 的代數(shù)方程組。這樣每個(gè)管段就可以提供包括四個(gè)未知量（兩個(gè)壓力和兩個(gè)流量）的兩個(gè)方程。需要求出每個(gè)時(shí)間步長(zhǎng)上所有的未知量，就必須要聯(lián)立所有管段上的這些方程。同時(shí)，還要考慮對(duì)應(yīng)的邊界條件，這樣才能使差分方程和邊界方程構(gòu)成封閉的方程組。

2.2 邊界條件

節(jié)點(diǎn)為系統(tǒng)的邊界點(diǎn)，節(jié)點(diǎn)的工藝要求就是動(dòng)態(tài)仿真的邊界條件。對(duì)于存在N個(gè)節(jié)點(diǎn)的天然氣管網(wǎng)來(lái)說(shuō)，每一個(gè)節(jié)點(diǎn)都必須給定一個(gè)對(duì)應(yīng)的工藝參數(shù)來(lái)作為邊界條件，即

式中：j為時(shí)間層，j=0,1,2,3……；Mi,j或Pi,j為 i節(jié)點(diǎn)j時(shí)刻的流量或壓力；Ci,j為i節(jié)點(diǎn)j時(shí)刻的邊界條件值。

2.3 初始條件

在管網(wǎng)進(jìn)行不同時(shí)間的工況運(yùn)行仿真時(shí)，還首先需要知道管網(wǎng)內(nèi)天然氣的初始時(shí)刻的運(yùn)行參數(shù)，即為初始條件。在大多數(shù)情況下，初始條件所給正確與否可以直接影響計(jì)算結(jié)果的真實(shí)性。本文是通過(guò)穩(wěn)態(tài)計(jì)算，用所得的穩(wěn)態(tài)計(jì)算結(jié)果作為動(dòng)態(tài)仿真的初值。

3 計(jì)算示例

3.1 編程計(jì)算結(jié)果

用隱式法計(jì)算單根天然氣管道的各個(gè)時(shí)間層的壓力和流量（不考慮高程）。單根管道的長(zhǎng)度為1km，管道內(nèi)徑為0.252984m，進(jìn)口壓力設(shè)為3MPa，初始流量為3kg/s，管道步長(zhǎng)設(shè)為200m，時(shí)間步長(zhǎng)設(shè)為1s，由初始時(shí)間各個(gè)點(diǎn)的壓力和流量以及邊界條件（表1），從而進(jìn)行迭代，求每一事件層的各點(diǎn)壓力和流量數(shù)據(jù)。

由所編程序可運(yùn)行出壓力和流量數(shù)據(jù)，如表2和表3。

表1 邊界條件

表2 壓力數(shù)據(jù) Pa

表3 流量數(shù)據(jù) kg/s

由以上數(shù)據(jù)對(duì)比可發(fā)現(xiàn)，每一時(shí)間層的壓力由初始?jí)毫?MPa開(kāi)始隨管道長(zhǎng)度的增加而遞減，每一時(shí)間層的壓力遞減程度都會(huì)越來(lái)越平緩，最后趨近于收斂層的趨勢(shì)；每一層的流量圍繞初始流量3kg/s波動(dòng)，在最后逐漸收斂于初始流量3kg/s。

3.2 PNS單根天然氣管道建立與運(yùn)行

由PNS管網(wǎng)仿真軟件建立單根的天然氣管道模型，數(shù)據(jù)和本文對(duì)應(yīng)的實(shí)例數(shù)據(jù)一致，詳見(jiàn)圖2、表 4。

圖2 管道模型

表4 數(shù)據(jù)對(duì)比

3.3 編程程序結(jié)果與PNS結(jié)果對(duì)比

由所編程序運(yùn)行的結(jié)果與PNS結(jié)果對(duì)比可以發(fā)現(xiàn)，其中流量和壓力所得的趨勢(shì)圖的趨勢(shì)一致，所編程序流量結(jié)果與PNS管道的運(yùn)行結(jié)果一致；所編程序壓力結(jié)果與PNS管道的運(yùn)行結(jié)果有所偏差，相對(duì)誤差最大為0.0591443%，誤差主要來(lái)自于對(duì)摩阻系數(shù)的計(jì)算。

4 結(jié)論

管道無(wú)高程最后所得的相對(duì)誤差都很小，證明該編程可以用于運(yùn)行，所得數(shù)據(jù)能夠提供精度高的數(shù)據(jù)進(jìn)行生產(chǎn)調(diào)度和管理。