有限容積法燃?xì)夤艿婪€(wěn)態(tài)模擬及分析

摘 要

采用可壓縮流體的高階有限容積法(FVM)及其系列改進(jìn)算法對(duì)燃?xì)夤艿肋M(jìn)行了穩(wěn)態(tài)分析，探討了數(shù)學(xué)模型的建立和數(shù)學(xué)模型的求解。
摘要：采用可壓縮流體的高階有限容積法(FVM)及其系列改進(jìn)算法對(duì)燃?xì)夤艿肋M(jìn)行了穩(wěn)態(tài)分析，探討了數(shù)學(xué)模型的建立和數(shù)學(xué)模型的求解。通過實(shí)例對(duì)比可知，F(xiàn)VM方法具有較快的運(yùn)算速度和較好的收斂性，計(jì)算結(jié)果較準(zhǔn)確。
關(guān)鍵詞：燃?xì)夤艿?；穩(wěn)態(tài)分析；水力計(jì)算；有限容積法
Steady-state Simulation and Analysis of Gas Pipeline Based on Finite Volume Method
LIU Xiaojing，ZHOU Weiguo，WANG Hai
Abstract：The steady-state analysis of gas pipeline is conducted using compressible fluid high-order finite volume method(FVM)and its improved algorithms.The establishment and solution of mathe-matical model are discussed.The ease comparison shows that FVM has rapider operation speed，better convergence and exact calculation result.
Key words：gas pipeline；steady-state analysis；hydraulic calculation；finite volume method(FVM)
1 概述
   天然氣管網(wǎng)穩(wěn)態(tài)分析是天然氣管網(wǎng)設(shè)計(jì)和改造的依據(jù)，是進(jìn)行管網(wǎng)動(dòng)態(tài)模擬和分析的基礎(chǔ)，也是加強(qiáng)天然氣管網(wǎng)系統(tǒng)優(yōu)化運(yùn)行以及確定最優(yōu)改擴(kuò)建方案的基礎(chǔ)。管道模擬準(zhǔn)確與否關(guān)鍵在于所建立的數(shù)學(xué)模型能否準(zhǔn)確、全面地描述管內(nèi)流動(dòng)過程以及能否找到快速求解模型且收斂性好的方法。
   在HARDY等人提出迭代法解非線性方程后，目前常用的穩(wěn)態(tài)分析方法有牛頓-拉夫遜(Newton-Raphson)法、線性逼近法、流體網(wǎng)絡(luò)理論。管延文、段常貴等都曾用牛頓-拉夫遜法進(jìn)行過計(jì)算[1～2]。牛頓-拉夫遜法是求解天然氣管網(wǎng)穩(wěn)態(tài)分析節(jié)點(diǎn)法數(shù)學(xué)模型的常用算法，但它有兩個(gè)主要缺點(diǎn)：①每次迭代都要計(jì)算雅可比矩陣及其逆矩陣，計(jì)算量大；②對(duì)迭代初始值要求比較高，如不能提供較好的初始值，則收斂較慢，甚至計(jì)算過程發(fā)散[3]。國(guó)內(nèi)的姜東琪等[4]和國(guó)外的KE[5]也采用線性逼近法對(duì)燃?xì)夤芫W(wǎng)進(jìn)行過水力計(jì)算。KE和TA0等人利用流體網(wǎng)絡(luò)理論對(duì)管網(wǎng)進(jìn)行了穩(wěn)態(tài)模擬[6～7]。西南石油大學(xué)對(duì)燃?xì)夤芫W(wǎng)的穩(wěn)態(tài)計(jì)算方法進(jìn)行了研究，分別應(yīng)用節(jié)點(diǎn)法和環(huán)路法進(jìn)行求解，但在計(jì)算中忽略了壓縮因子對(duì)氣體流動(dòng)的影響[8]。進(jìn)行管網(wǎng)水力計(jì)算時(shí)，必須考慮管道中燃?xì)獾目蓧嚎s性。
    本文采用SIMPLE改進(jìn)算法，用有限容積法(FVM)導(dǎo)出的離散方程可以保證具有守恒特性，而且離散方程系數(shù)的物理意義明確。這克服了采用有限差分法求解管網(wǎng)流動(dòng)偏微分方程的多種困難；也可以克服管網(wǎng)計(jì)算中各管段管徑不同，以及在連接處流量突變的困難。對(duì)不同的管段進(jìn)行長(zhǎng)度不同的離散容積劃分，不會(huì)影響計(jì)算結(jié)果的連續(xù)性和迭代的收斂穩(wěn)定性。而傳統(tǒng)方法離散時(shí)，為了保證計(jì)算結(jié)果收斂，必須對(duì)長(zhǎng)管段進(jìn)行很細(xì)密的網(wǎng)格劃分，導(dǎo)致計(jì)算時(shí)間很長(zhǎng)且無(wú)法保證數(shù)值解的守恒性。根據(jù)有限容積法離散方程的特點(diǎn)，可以采用高效的TDMA算法，存儲(chǔ)量很小[9]。
2 數(shù)學(xué)模型的建立
   ① 管道數(shù)學(xué)模型
   管道內(nèi)氣體流動(dòng)數(shù)學(xué)模型由動(dòng)量方程、連續(xù)性方程、狀態(tài)方程組成[10]，見式(1)。因?yàn)楣艿纼?nèi)溫度變化很小，所以不考慮能量方程。
 
式中ρ——燃?xì)饷芏?，kg/m3
    u——燃?xì)馑俣龋?span>m/s
    x——管道軸向坐標(biāo)，m
    p——燃?xì)鈮毫Γ?span>Pa
    g——重力加速度，m/s2
    ρa——空氣密度，kg/m3
    α——管道與水平面夾角，rad
    λ——管道摩擦阻力系數(shù)
    d——管道內(nèi)徑，m
    z——壓縮因子
    R——氣體常數(shù)，J/(kg·K)
    T——燃?xì)鉁囟龋?span>K
   ② 模型的簡(jiǎn)化
   a. 在動(dòng)量方程中，對(duì)流項(xiàng)只在燃?xì)饬魉贅O大(接近聲速)時(shí)才有意義，而通常管道中燃?xì)饬魉俨淮笥?0m/s。此外，在高壓管網(wǎng)瞬態(tài)變化中，對(duì)流項(xiàng)相對(duì)于其他項(xiàng)數(shù)值較小，可以忽略[11]。b.在城市燃?xì)夤芫W(wǎng)中，管道的標(biāo)高差值不太大，因此，動(dòng)量方程中的重力項(xiàng)g(ρ-ρa)sinα一般忽略不計(jì)。
    通過簡(jiǎn)化，可以得到仿真數(shù)學(xué)模型，見式(2)。
 
式中C——常量
   ③ 摩擦阻力系數(shù)
   管道摩擦阻力系數(shù)除與燃?xì)庑再|(zhì)、管道材質(zhì)等因素有關(guān)外，還與燃?xì)獾牧鲃?dòng)狀態(tài)有關(guān)。不分低壓和高中壓管道，統(tǒng)一按流動(dòng)狀態(tài)和管道材質(zhì)來分類選取計(jì)算公式，進(jìn)行燃?xì)夤艿滥Σ磷枇ο禂?shù)的計(jì)算[12]。
3 數(shù)學(xué)模型的求解
   ① 數(shù)學(xué)模型的求解
   在計(jì)算氣體穩(wěn)定流動(dòng)時(shí)，利用有限容積法一維交錯(cuò)網(wǎng)格將壓力-速度耦合方程中不同的變量離散式存儲(chǔ)在不同的網(wǎng)格系統(tǒng)，將式(2)中的動(dòng)量方程進(jìn)行離散，見式(3)。i節(jié)點(diǎn)是速度控制容積的中心節(jié)點(diǎn)，I節(jié)點(diǎn)是壓力控制容積的中心節(jié)點(diǎn)。
 
式中ai——i節(jié)點(diǎn)的離散方程系數(shù)，kg/s
    ui——i節(jié)點(diǎn)的速度，m/s
    n——i節(jié)點(diǎn)周圍的節(jié)點(diǎn)數(shù)
    aj——i節(jié)點(diǎn)周圍各節(jié)點(diǎn)的離散方程系數(shù)，kg/s
    uj——i節(jié)點(diǎn)周圍各節(jié)點(diǎn)的速度，m/s
    pI-1——I-1節(jié)點(diǎn)的壓力，Pa
    pI——I節(jié)點(diǎn)的壓力，Pa
    Ai——速度控制容積的界面面積，m2
    bi——方程源項(xiàng)，(kg·m)/s2
   在初始?jí)毫ο拢?/div>
 
式中ui*——初始?jí)毫ο耰節(jié)點(diǎn)的速度，m/s
    uj*——初始?jí)毫ο拢?span>i節(jié)點(diǎn)周圍各節(jié)點(diǎn)的速度，m/s
    pI-1*——I-1節(jié)點(diǎn)的初始?jí)毫Γ?span>Pa
    pI*——I節(jié)點(diǎn)的初始?jí)毫Γ?span>Pa
    設(shè)修正的壓力方程為：
    p=p*+p′    (5)
式中p′——壓力修正值，Pa
    設(shè)修正的速度方程為：
    u=u*++u′    (6)
式中u′——速度修正值，m/s
    由式(4)～(6)得：
 
式中ui′——i節(jié)點(diǎn)速度修正值，m/s
    uj′——i節(jié)點(diǎn)周圍各節(jié)點(diǎn)速度修正值，m/s
    pI-1′——I-1節(jié)點(diǎn)的壓力修正值，Pa
    pI′——I節(jié)點(diǎn)的壓力修正值，Pa
由于式(7)第一項(xiàng)是周圍節(jié)點(diǎn)速度引起的修正值，值很?。欢诙?xiàng)是同一方向相鄰節(jié)點(diǎn)壓力差引起的修正值，影響較大。故為簡(jiǎn)化計(jì)算，略去第一項(xiàng)，并由式(6)、(7)得到速度的改進(jìn)值：
 
   將速度改進(jìn)值代入方程組(2)的連續(xù)性方程中，并將壓力修正值p′的系數(shù)歸一化處理，得到壓力修正方程：
    aIpI′=aI+1pI+1′+aI-1pI-1′+bI′    (9)
式中aI——I節(jié)點(diǎn)離散方程系數(shù)，kg/m
    aI+1——I+1節(jié)點(diǎn)離散方程系數(shù)，kg/m
    pI+1′——I+1節(jié)點(diǎn)的壓力修正值，Pa
    aI-1——I-1節(jié)點(diǎn)離散方程系數(shù)，kg/m
    bI′——由于速度場(chǎng)的不正確引起的不平衡量，kg/s
   經(jīng)過多次迭代修正，最終bI′是否趨于0可以作為判斷迭代過程是否滿足要求的判據(jù)。
   在計(jì)算時(shí)，管道內(nèi)氣體壓力和流速的分布采用雙精度Simple方法的PIS0模型進(jìn)行求解。即通過壓力預(yù)測(cè)一修正方法，包括1個(gè)預(yù)測(cè)步驟和2個(gè)校正步驟，不斷地修正計(jì)算結(jié)果，反復(fù)迭代，最后求出壓力和速度的收斂解。
   ② 程序編制
   采用VB編寫計(jì)算程序，對(duì)網(wǎng)格參數(shù)、壓力參數(shù)內(nèi)循環(huán)迭代的計(jì)算均采用穩(wěn)定且快速收斂的TDMA方法。
4 實(shí)際應(yīng)用與分析
   為驗(yàn)證燃?xì)夤芫W(wǎng)仿真軟件是否準(zhǔn)確，英國(guó)氣體公司倫敦研究站(LRS)采用了一組不同管徑、流量下符合實(shí)際的燃?xì)獾姆€(wěn)態(tài)數(shù)據(jù)進(jìn)行驗(yàn)證。在相同的管徑、流量條件下，本文采用FVM方法進(jìn)行穩(wěn)態(tài)仿真，并將仿真結(jié)果與新加坡國(guó)立大學(xué)采用電網(wǎng)絡(luò)傳輸線理論的模擬結(jié)果進(jìn)行對(duì)比。
   p1，p2，p3，p4，p5分別表示管道被平均分為4段后，5個(gè)節(jié)點(diǎn)處的壓力值。FVM方法的仿真結(jié)果與LRS的數(shù)據(jù)對(duì)比見表1，新加坡國(guó)立大學(xué)的電網(wǎng)絡(luò)傳輸線理論模擬結(jié)果與LRS的數(shù)據(jù)對(duì)比見表2。
表1 FVM方法仿真結(jié)果與LRS數(shù)據(jù)對(duì)比
管道長(zhǎng)度/km
管道內(nèi)徑/m
流量/(m3·s-1)
相對(duì)誤差/%
p1
p2
p3
p4
p5
128.747
0.457
0.408
0.000
0.057
0.074
0.088
0.004
128.747
0.457
0.612
0.000
0.100
0.182
0.436
0.009
32.187
0.305
0.218
0.000
0.206
0.105
0.237
0.O01
32.187
0.305
0.327
0.000
0.355
0.960
2.175
0.045
9.656
0.356
0.082
0.000
0.072
0.109
0.166
0.000
9.656
0.356
0.122
0.000
0.123
0.387
1.069
0.000
4.828
0.356
0.022
0.000
0.000
0.117
0.322
0.000
4.828
0.356
0.033
0.000
0.059
0.220
0.388
0.000
0.805
0.102
0.001
0.000
0.000
0.152
0.000
0.000
0.805
0.102
0.002
0.000
0.153
0.000
0.196
0.000
表2 新加坡國(guó)立大學(xué)的電網(wǎng)絡(luò)傳輸線理論模擬結(jié)果與LRS的數(shù)據(jù)對(duì)比
管道長(zhǎng)度/km
管道內(nèi)徑/m
流量/(m3·s-1)
相對(duì)誤差/%
p1
p2
p3
p4
p5
128.747
0.457
0.408
0.000
0.000
0.180
0.250
0.170
128.747
0.457
0.612
0.000
0.060
0.130
0.040
0.170
32.187
0.305
0.218
0.000
0.410
0.560
0.940
0.720
32.187
0.305
0.327
0.000
0.380
0.890
1.510
0.300
9.656
0.356
0.082
0.000
0.850
1.360
2.450
3.210
9.656
0.356
0.122
0.000
1.830
3.280
4.170
3.460
4.828
0.356
0.022
0.000
0.000
0.590
0.650
1.870
4.828
0.356
0.033
0.000
0.600
1.470
1.940
0.290
0.805
0.102
0.001
0.000
0.000
0.000
0.580
0.610
0.805
0.102
0.002
0.000
0.000
0.620
0.710
1.640
    由表1可知，運(yùn)用FVM方法的仿真結(jié)果與LRS的數(shù)據(jù)非常接近，平均相對(duì)誤差小于1.000%。表1與表2對(duì)比可知，F(xiàn)VM方法的仿真結(jié)果更準(zhǔn)確。采用FVM方法具有較快的運(yùn)算速度和較好的收斂性，為管網(wǎng)的動(dòng)態(tài)分析提供了基礎(chǔ)，可以滿足工程的需要，對(duì)于管網(wǎng)的設(shè)計(jì)和安全控制都有重要的參考作用。
參考文獻(xiàn)：
[1] 管延文，王永剛，李統(tǒng)強(qiáng)，等.牛頓管段方程法的燃?xì)夤芫W(wǎng)穩(wěn)態(tài)分析[J].煤氣與熱力，2006，26(6)：10-12.
[2] 段常貴，徐彥峰，呂文哲，等.燃?xì)夤芫W(wǎng)的穩(wěn)態(tài)分析[J].煤氣與熱力，2000，20(2)：95-97.
[3] 白建輝，汪玉春，郜峰，等.天然氣管網(wǎng)穩(wěn)態(tài)分析綜合方法[J].油氣儲(chǔ)運(yùn)，2009，28(2)：37-39.
[4] 姜東琪，杜建梅，崔建華，等.燃?xì)夤芫W(wǎng)水力計(jì)算及水力計(jì)算圖的繪制[J].煤氣與熱力，2001，21(5)：453-455.
[5] KE S L.Transient analysis of pipeline network(master’s theses)[D].Singapore：National University of Singapore.1999：54-60.
[6] KE S L.TI H C.Transient analysis of isothermal gas flow in pipeline network[J].Chemical Engineering Journal，2000，76(2)：169-177.
[7] TAO W Q，TI H C.Transient analysis of gas pipeline network[J].Chemical Engineering Journal，1998，69(3)：47-52.
[8] 李長(zhǎng)明.燃?xì)夤芫W(wǎng)水力計(jì)算程序設(shè)計(jì)基礎(chǔ)[M].北京：煤炭工業(yè)出版社，1997：10-80.
[9] 李人憲.有限體積法基礎(chǔ)[M].北京：國(guó)防工業(yè)出版社，2001：81-83.
[10] OSIADACZ A J.Different transient models——limitations.advantages and disadvantages[C]∥.28th Annual Meeting PSIGs.San Francisco，California：Pipe]ine Simulation Interest Group，1996：1-26.
[11] 段常貴.燃?xì)廨斉鋄M].3版.北京：中國(guó)建筑工業(yè)出版社，2001：87-88.
[12] 蔣祥龍.燃?xì)夤艿滥Σ磷枇τ?jì)算的探討[J].煤氣與熱力，2004，24(1)：37-40.
 
(本文作者：劉曉婧 周偉國(guó) 王海 同濟(jì)大學(xué)機(jī)械工程學(xué)院 上海 201804)