李春奇，袁獻(xiàn)忠，李 俊，郭 喬

(1.中國石油化工股份有限公司天然氣分公司，北京 100029； 2.國家管網(wǎng)集團(tuán)川氣東送天然氣管道有限公司，湖北 武漢 430223； 3.西安石油大學(xué) 石油工程學(xué)院，陜西 西安 710065)

引 言

作為一種清潔、高效的化石能源，天然氣的利用受到世界各國的重視，長距離管道運輸無疑是最方便、經(jīng)濟(jì)的天然氣運輸方式。由于季節(jié)、天氣、用戶數(shù)量、突發(fā)事件等因素的影響，管道內(nèi)的天然氣流動一般是不穩(wěn)定和瞬態(tài)的，對于愈加大型化的天然氣管網(wǎng)，管道內(nèi)的流動更加復(fù)雜。對流動參數(shù)的準(zhǔn)確預(yù)測是保證天然氣管道安全運行的前提。在這種情況下，管網(wǎng)仿真的作用就顯得尤為重要[1]。

輸氣管道流動狀態(tài)可由質(zhì)量守恒方程、動量守恒方程和能量守恒方程描述，有多種方法可求解三大守恒方程。由于特征線法[2]、MacCormack[3]法均受制于穩(wěn)定條件的約束，時間步長只能取得很小，如果仿真過程較長，則需要消耗較長的計算時間。而隱式差分法計算步長比較方便，能夠選取較大的時間步長，減少計算次數(shù)和時間[4]。

對于隱式差分法，存在多種差分格式。國內(nèi)眾多學(xué)者[5-12]將不穩(wěn)定流方程以有限差形式應(yīng)用于位居四點網(wǎng)格中心的點子。Kiuchi[13]，Abbaspour和Chapman[14]，對時間項采用向前差分進(jìn)行離散，對流項采用中心差分進(jìn)行離散，并稱這種差分格式為Implicit Cell Centered Method；鄭建國[15]，王鵬[16-18]，向月[19-20]在Implicit Cell Centered Method的基礎(chǔ)上，對偏微分方程組的非線性項進(jìn)行泰勒展開，保留一階項忽略高階項。不同的差分格式直接影響管網(wǎng)仿真的計算精度，不同的差分格式適用于不同的求解方式，進(jìn)而影響仿真計算速度，因此，尋找方程組合適的差分格式至關(guān)重要。

1 仿真理論

1.1 管道模型

管道模型用于描述管道內(nèi)的流動狀態(tài)，對于管道橫截面積不發(fā)生變化的管道，可建立質(zhì)量守恒方程、動量守恒方程和能量守恒方程[7]，并將其統(tǒng)一寫為向量形式[16]

(1)

其中各參數(shù)見表1。

表1 管道數(shù)學(xué)模型相應(yīng)變量的含義Tab.1 Meaning of variables in pipeline mathematical model

其中，p為流體壓力，Pa；m為質(zhì)量流量，kg/s；T為流體溫度，K；Tg為環(huán)境溫度，K；A為管道截面積，m2；ρ為流體密度，kg/m3；t為時間變量，s；θ為管段傾角，rad；w為流體流速，m/s；x為沿管道變量，m；g為重力加速度，m/s2；d為管道內(nèi)徑，m；D為管段外徑，m；z為高程，m；cv為流體比熱容，J/(kg·K)；K為總傳熱系數(shù)，W/(m2·K)；λ為摩阻系數(shù)，無因次，可由柯爾布魯克公式計算，即

(2)

式中：Re為雷諾數(shù)，無因次。狀態(tài)方程選擇BWRS狀態(tài)方程。

1.2 隱式差分法

為求解式進(jìn)行離散化，使用不同的差分格式得到不同的離散方程，進(jìn)而得到不同的計算結(jié)果，因此研究中心有限差分法、Implicit Cell Centered Method以及Implicit Cell Centered Method線性化法3種方法的差分格式以及相容性。

1.2.1 差分格式

將管道在空間上劃分為N等分，每一等分長度Δx，在時間上每一時間步長為Δt，縱坐標(biāo)表示時間，橫坐標(biāo)表示空間，如圖1所示。

(1)方法1：中心有限差分法

中心有限差分法是將不穩(wěn)定流方程以有限差形式應(yīng)用于位居四點網(wǎng)格中心的點子，離散格式為

(3)

(4)

(5)

圖1 隱式差分法網(wǎng)格圖Fig.1 Grid diagram of implicit difference method

(2)方法2：Implicit Cell Centered Method

Implicit Cell Centered Method是對時間項采用向前差分進(jìn)行離散，對流項采用中心差分進(jìn)行離散，離散格式為

(6)

(7)

(8)

(3)方法3：Implicit Cell Centered線性化

Implicit Cell Centered線性化是在Implicit Cell Centered方法離散格式的基礎(chǔ)上，對式中的B、F在上一時間層處進(jìn)行泰勒展開，并忽略2階無窮小項，具體格式為

(9)

1.2.2 相容性分析

相容性意味著當(dāng)時、空步長均趨近于零時，離散方程逼近微分方程，也就是說當(dāng)時、空步長趨近于零時，如果離散方程的截斷誤差趨近于零，則稱此離散方程與微分方程相容[21]。對于管道模型，3種不同方法的相容性如下：

定義微分算子：

(10)

對于方法1，定義差分算子

(11)

因此，方法1的截斷誤差為

(12)

同理，方法2的截斷誤差為

(13)

方法3的截斷誤差為

(14)

對于3種差分格式，當(dāng)時、空步長趨近于零時，其對應(yīng)的截斷誤差也趨近于零，因此3種格式都具有相容性。但是由于方法3采用了泰勒展開，這就導(dǎo)致該方法要求相鄰時間層的B、F、U相差不能過大，否則將存在較大的截斷誤差。

1.3 方程組求解方法

將管道劃分為n個管段后，針對每個管段建立質(zhì)量守恒、動量守恒以及能量守恒方程，并聯(lián)立邊界條件，即可求解管道中每個節(jié)點的壓力、流量和溫度。一般而言，求解非線性方程組要比求解相同維度的線性方程組的計算速度慢，方法1和方法2所建立的方程組是非線性方程組，可以使用牛頓拉夫遜法求解，對于方法3則是將非線性項線性化后，可建立稀疏線性方程組，求解該稀疏線性方程組，因此方法3的計算速度較快。

2 結(jié)果對比

2.1 算例描述

為了對比仿真結(jié)果，分別以3種差分格式應(yīng)用于同一個仿真案例。水平管道內(nèi)徑600 mm，管道長度為5 km，粗糙度為0.025 mm，標(biāo)準(zhǔn)狀態(tài)：壓力為101.325 kPa，溫度為20 ℃。天然氣組分見表2。

表2 天然氣主要組分Tab.2 Main composition of natural gas

初始時刻管道處于穩(wěn)定狀態(tài)：進(jìn)口壓力為6 MPa，溫度為30 ℃，流量為0 Nm3/h。零時刻以后，進(jìn)口恒定壓力6 MPa，在第6 s時刻出口流量變?yōu)?×105Nm3/h，保持到第18 s，并瞬變?yōu)? Nm3/h，直至36 s仿真結(jié)束，如圖2所示。

圖2 邊界條件Fig.2 Boundary conditions

2.2 計算精度

為了對比計算精度，分別以3種差分格式求解上述仿真案例，分別設(shè)置空間步長、時間步長為0.2 km、0.6 s，并且在相同的條件下與PNS油氣管網(wǎng)仿真軟件[22]計算結(jié)果對比，分析第12 s、24 s、30 s以及36 s時刻的流量、壓力分布(圖3)。

圖3 不同時刻的流量、壓力分布圖Fig.3 Flow rate and pressure distribution at different times

從圖3可以看出，在第12 s、24 s以及36 s時刻，方法2、方法3兩種方法的流量分布以及壓力分布基本重合，且與PNS仿真結(jié)果吻合度較高。而方法1計算的流量和壓力分布與PNS計算結(jié)果存在較大偏差，并且無論在哪個時刻，方法1所計算的流量分布和壓力分布均存在不符合客觀規(guī)律的波動。

2.3 計算速度

在求解管道動態(tài)仿真時，以第i時層的水力熱力參數(shù)作為初始值，求解第i+1時層的，而計算速度又與方程組初始值有關(guān)，為了避免初始值不同對求解速度的影響，因此只觀察管道流動突變時刻時的單層計算速度，此時方程組的初始值均相同，即上一時層的計算結(jié)果。

不同的分段數(shù)意味著構(gòu)建不同規(guī)模的方程組，分段數(shù)越大，方程組維數(shù)越多。為了對比不同分段數(shù)下的計算速度，分別對比將管道劃分為5段、10段、20段和25段時的仿真耗時(圖4)。

圖4 3種方法仿真耗時Fig.4 Time consumption of three simulation methods

從圖4中可以看出，無論在哪一時刻，方法3計算耗時基本為0，計算速度遠(yuǎn)高于另外兩種，而方法2與方法1計算速度相差無幾。這是由于方法1與方法2需要通過迭代求解非線性方程組，并在迭代的過程中，需要更新流體物性，這就導(dǎo)致每迭代一次的耗時不僅包含非線性方程組的求解耗時，也包括每次迭代時的流體物性更新耗時，仿真耗時較大；而方法3直接求解線性方程組，無需迭代也無需更新物性，故而計算速度較快。

3 結(jié) 論

(1)相容性：3種方法均具有相容性，即當(dāng)時間步長和空間步長趨近于零時，離散截斷誤差也趨近于零。

(2)計算精度：中心差分法計算精確性較差，無法正確描述天然氣管道的動態(tài)仿真過程，而Implicit Cell Centered和Implicit Cell Centered線性化法計算精度相當(dāng)。

(3)計算速度：由于中心差分法與Implicit Cell Centered法計算速度較慢，而Implicit Cell Centered線性化法直接求解線性方程組，并且計算過程中不需要更新物性，計算速度非?？欤试诤笃诘奶烊粴夤芫W(wǎng)動態(tài)仿真研究中推薦使用Implicit Cell Centered線性化法。