丁延鵬，李玉星，郝宏娜，畢研軍

(1.中國石油大學(xué)儲運與建筑學(xué)院，山東青島 266555;2.中國石油天然氣管道科學(xué)研究院，河北廊坊 065000)

地層蓄熱效應(yīng)對輸氣管道數(shù)值模擬精度的影響

丁延鵬1，李玉星1，郝宏娜1，畢研軍2

(1.中國石油大學(xué)儲運與建筑學(xué)院，山東青島 266555;2.中國石油天然氣管道科學(xué)研究院，河北廊坊 065000)

采用一維非等溫天然氣流動模型對高壓天然氣輸送管線中的快瞬變流動進行數(shù)值模擬，研究地層蓄熱效應(yīng)對高壓天然氣快瞬變管流數(shù)值模擬精度的影響。從工程實際角度對質(zhì)量、動量和能量守恒方程的基本變量進行等價變換，建立天然氣管流與外界地層之間一維軸對稱的非穩(wěn)態(tài)傳熱模型，以BWRS狀態(tài)方程為基礎(chǔ)，采用特征線法對天然氣管流的快瞬變過程數(shù)學(xué)模型進行數(shù)值求解，并與簡化穩(wěn)態(tài)模型數(shù)值結(jié)果進行對比。結(jié)果表明：國內(nèi)常用的簡化穩(wěn)態(tài)傳熱模型過高估計了管線中流量的波動幅度，工程應(yīng)用中應(yīng)采用考慮地層蓄熱效應(yīng)的非穩(wěn)態(tài)傳熱模型對高壓天然氣管流進行仿真模擬;非穩(wěn)態(tài)模型和算法具有較高的準確性。

輸氣管道;快瞬變流;數(shù)值模擬;地層蓄熱效應(yīng);特征線法

由于輸氣管道運行中存在著天然氣的需求波動，需要進行閥門開關(guān)、壓縮機啟停、調(diào)壓器調(diào)節(jié)等操作，因而輸氣管道往往處于不穩(wěn)定的運行狀態(tài)。投產(chǎn)時間長、腐蝕磨損嚴重的管網(wǎng)在閥門開關(guān)等操作過程中就存在著事故隱患，因此有必要對輸氣管道在閥門開關(guān)過程中的運行狀況進行瞬態(tài)模擬分析。天然氣管流的瞬態(tài)模擬是指對描述管內(nèi)流體質(zhì)量、動量和能量守恒以及狀態(tài)變化的偏微分方程組進行數(shù)值求解。目前，國內(nèi)外對天然氣管流瞬態(tài)模擬的研究大都采用一維非穩(wěn)態(tài)數(shù)學(xué)模型，國外的研究重點在于數(shù)學(xué)模型以及求解方法上的創(chuàng)新與改進，國內(nèi)在此領(lǐng)域的研究則側(cè)重于模擬仿真技術(shù)在輸氣管道設(shè)計運行中的應(yīng)用。管流與地層之間的傳熱通過能量方程中的傳熱項進行表述，天然氣管流模型中傳熱項的處理存在著兩種極端情況：等溫模型和絕熱模型［1］。目前，介于等溫與絕熱模型之間的非等溫傳熱模型應(yīng)用比較廣泛［2－4］，依據(jù)傳熱過程是否隨時間變化分為穩(wěn)態(tài)和非穩(wěn)態(tài)兩種。穩(wěn)態(tài)非等溫傳熱模型在我國的輸氣管道模擬仿真研究中處于主導(dǎo)地位。筆者在一維非等溫天然氣管流模型中加入地層蓄熱效應(yīng)，建立非穩(wěn)態(tài)非等溫傳熱模型，對閥門瞬時開關(guān)產(chǎn)生的快瞬變流進行數(shù)值模擬研究。

1 數(shù)學(xué)模型

1.1 氣體管流瞬變模型

由于天然氣管輸系統(tǒng)中經(jīng)常存在著由于需求波動以及系統(tǒng)控制裝置(如閥門、壓縮機和調(diào)壓器等)調(diào)節(jié)等原因形成的動態(tài)管流，因此實際模擬須采用考慮到壓力、溫度和流量等突變的雙曲型非等溫天然氣管流模型。包含質(zhì)量、動量和能量守恒定律的描述天然氣管道內(nèi)一維可壓縮流動的數(shù)學(xué)模型［5－6］如下：

式中，p為管道壓力，Pa;w為氣體流速，m/s;ρ為氣體密度，kg/m3;x為管內(nèi)軸向長度，m;u為比內(nèi)能，J/kgh;h為氣體比焓，J/kg。

上述天然氣管流模型是一組以密度、壓力和流速為基本變量的雙曲型偏微分方程組，比較適用于管線參數(shù)發(fā)生突然變化的情況［1，7］。天然氣工業(yè)中，氣體流速的測量數(shù)據(jù)并不容易得到，而且密閉輸送中標(biāo)準工況下的體積流量或質(zhì)量流量才是天然氣輸送過程常用的參數(shù)。因此，將質(zhì)量、動量和能量守恒方程變換成以壓力、標(biāo)況體積流量和溫度為基本變量的等價形式［9－10］：

式中，A為管道橫截面積，m2;D為管內(nèi)徑，m;Qm為體積流量，m3/s;s為管道沿線高程，m;g為重力加速度，m/s2;λ為摩阻系數(shù);t為時間，s;q為單位時間內(nèi)單位質(zhì)量的氣體與外界環(huán)境間的傳熱量，kJ/(kg·s)。

方程(6)為管線內(nèi)控制體內(nèi)單位質(zhì)量氣體的能量守恒定律，其最后一項是管內(nèi)氣體與管道外界環(huán)境間的傳熱項，該項對上述模型求解氣體參數(shù)影響較大［1，11］。

1.2 狀態(tài)方程

考慮到方程組(4)～(6)的封閉性及內(nèi)在變量間的水力熱力關(guān)聯(lián)，需加入狀態(tài)方程。BWRS狀態(tài)方程屬于多參數(shù)狀態(tài)方程，具有很寬的操作工況，對比溫度可低至Tr=0.3，對比密度則可高達ρr=0.3，并可以用于計算含有CO2、H2S、N2等氣體的體系，是目前最佳的天然氣狀態(tài)方程之一。BWRS狀態(tài)方程可以將壓力表達成溫度和密度的顯示函數(shù)［11］，即

式中，T為系統(tǒng)溫度，K;R為氣體常數(shù)，J/(kg·K);A0、B0、C0、D0、E0、a、b、c、d、α、γ 為狀態(tài)方程中的參數(shù)。

1.3 傳熱模型

在能量方程(6)中，傳熱項q代表著單位時間內(nèi)單位質(zhì)量的氣體與外界環(huán)境間的傳熱量。對離散的管元，采用傅里葉法則計算氣體和地層間的總傳熱量：

式中，K為傳熱系數(shù)，W/(m·K);Tamb為管道軸線兩側(cè)足夠遠處環(huán)境溫度，K。

導(dǎo)熱是天然氣快瞬變管流中起主導(dǎo)作用的傳熱機制［12］，因此管道內(nèi)沿軸向和環(huán)向的熱傳遞可以忽略不計，只須考慮管壁和地層中的導(dǎo)熱。通常的工程實踐中，描述輸氣管道與外界地層間傳熱的模型通常是一維的，因為缺少描述管線沿線附近地層構(gòu)造以及地層性質(zhì)(溫度、傳熱特性等)的準確數(shù)據(jù)，因此大都假定管道與外界地層間的熱阻恒定，并進行簡化穩(wěn)態(tài)傳熱處理：假定管道周圍地層的環(huán)境溫度恒定Tamb不變。

將輸氣管道與外界地層間的傳熱項進行非穩(wěn)態(tài)處理：由于壓縮機出口管線的溫度較高，對周圍地層有一定加熱效果，考慮到管道周圍地層的蓄熱效應(yīng)，周圍地層中會形成一個以管道為軸線的徑向溫度變化，如圖1所示。

圖1 非穩(wěn)態(tài)傳熱模型離散示意圖Fig.1 Dispersion diagram of unsteady heat transfer model

非穩(wěn)態(tài)一維軸對稱的徑向傳熱過程用方程組表述為

式中，n為管道周圍地層傳熱區(qū)的離散元數(shù)目;mi為第 i層離散元的質(zhì)量(i=1，…，n)，kg;cp，i為第 i單元的比定壓熱容，kJ/(kg·K);dx為管元長度，m;Ti為第i地層單元的溫度，K;ki為第i－1層與第i層地層單元之間的傳熱系數(shù)，W/(m·K)。

在采用一維模型進行數(shù)值模擬的過程中，傳熱模型最少應(yīng)含兩層圓柱型傳熱層。假定不同地層離散單元的蓄熱特性差別很大，時間常數(shù)也就不同，管道附近和遠處的環(huán)境對溫度變化的反應(yīng)也分別有快有慢。非穩(wěn)態(tài)傳熱過程中，考慮地層蓄熱的初始條件設(shè)定為管道附近相鄰地層單元之間的溫度分布與簡化穩(wěn)態(tài)情況一致。

2 數(shù)值求解

采用的管流模型為一組雙曲型偏微分方程組，因而采用特征線法進行數(shù)值求解。特征線法是求解二元一階雙曲型偏微分方程(組)的一類經(jīng)典方法，其基本思想是先將基本方程組變換為沿特征線的常微分方程組，然后對常微分方程組結(jié)合初始條件和邊界條件進行數(shù)值解。特征線法的優(yōu)勢是它可以處理模擬中的不連續(xù)性，比較方便地對邊界條件進行處理，物理概念明確、數(shù)學(xué)分析嚴謹、計算精度較高［13］。

將方程組(4)～(6)采用特征線法進行基本變換之后得到沿軌線、順特征線和逆特征線的3個特征方程：

沿軌線的特征方程為

式中，T0為管壁溫度，K;as為音速，m/s。

方程組(9)～(11)是對原來的3個基本方程進行組合后得到的，因而已分不清連續(xù)性方程、動量方程和能量方程了，但它在整體上仍然與原來的基本微分方程組等價。采用以逆步進法為基礎(chǔ)的有限差分法進行離散求解。離散網(wǎng)格如圖2所示。在這個方法中解點A的位置通過網(wǎng)格劃分確定。交于A點上的3條特征線向后延伸與前一時刻計算的一條線相交，以確定初始數(shù)據(jù)點G、M和H的位置。初始數(shù)據(jù)點上的流動參數(shù)用先前的解點進行插值獲得。

圖2 特征線法差分網(wǎng)格示意圖Fig.2 Difference grid diagram of method of characteristics

3 算例

3.1 條件設(shè)定

采用Kiuchi提出的輸氣系統(tǒng)模型［2］對本文的模擬結(jié)果進行驗證與分析。該系統(tǒng)如圖3所示，包括一段長度為L=5 km、內(nèi)徑D=0.5 m的管線。地溫Tamb=25℃，管內(nèi)流體與管外地層間的總傳熱系數(shù)(包括管內(nèi)流體與管壁間的對流換熱、管壁內(nèi)的導(dǎo)熱以及管壁與緊貼管壁的第一層地層單元間的換熱)K=2.84 W/(m2·K)。

圖3 輸氣管道簡圖Fig.3 Diagram of gas transmission pipeline

式中，t為時間，min。

該模型通過入口處的流量與溫度變化模擬壓縮機出口閥門瞬時開關(guān)操作對管內(nèi)參數(shù)的影響，通常用于檢驗管道快瞬變工況模擬算法的準確性。

本文中假定管道附近地層的各離散單元具有相同的熱阻與熱容。初始狀態(tài)下，相鄰地層單元之間的溫度分布與穩(wěn)態(tài)的情況一致。

3.2 模型算法對比驗證

先用簡化的穩(wěn)態(tài)傳熱模型對上述輸氣系統(tǒng)進行模擬，并與Abbaspour計算結(jié)果進行對比，入口處流量模擬結(jié)果見圖4。由圖4可以發(fā)現(xiàn)，本文模型算法的結(jié)果與Abbaspour的計算結(jié)果［14］總體上比較吻合，證明了本模型及算法的準確性。

圖4 模擬結(jié)果對比驗證Fig.4 Comparison of simulation results

3.3 不同傳熱模型的模擬結(jié)果對比

采用考慮地層蓄熱效應(yīng)的非穩(wěn)態(tài)傳熱模型，對Kiuchi輸氣系統(tǒng)進行了數(shù)值模擬，并將簡化穩(wěn)態(tài)傳熱模型與非穩(wěn)態(tài)傳熱模型對入口流量的模擬結(jié)果進行比較，如圖5、6所示。

圖5 簡化穩(wěn)態(tài)與非穩(wěn)態(tài)傳熱模型的模擬結(jié)果對比Fig.5 Comparison of simulation results using simplified steady and unsteady heat transfer model

圖6 閥門開關(guān)過程流量波動局部放大圖Fig.6 Partial enlargement of flow fluctuation when opening/closing valve

從圖6可以看出，閥門開啟和關(guān)閉的模擬過程中，非穩(wěn)態(tài)傳熱模型算得的流量比簡化穩(wěn)態(tài)傳熱模型的模擬結(jié)果更快地到達穩(wěn)定狀態(tài)，這可以通過地層的蓄熱效應(yīng)來解釋。壓縮機對氣體的加熱使得管內(nèi)氣體溫度高于外界地層溫度，對管外地層進行了加熱，由于地層蓄熱效應(yīng)的存在，在輸氣管線周圍形成了軸向溫度梯度分布，地層溫度的升高降低了管內(nèi)氣體溫度的下降幅度，與簡化穩(wěn)態(tài)傳熱模型中設(shè)定外界地層溫度不變相比，降低了計算過程中溫度跳躍幅度，使得計算過程可以更快地到達穩(wěn)定狀態(tài)。這表明，考慮地層蓄熱效應(yīng)的非穩(wěn)態(tài)模型可以顯著增加流量曲線隨時間變化的準確性。

4 結(jié)束語

將BWRS狀態(tài)方程與一維非等溫天然氣管流瞬變模型耦合，采用特征線法對數(shù)學(xué)模型進行特征變換，應(yīng)用逆步進法進行離散數(shù)值求解，驗證結(jié)果說明該模型和算法具有較高的準確性。不考慮地層蓄熱效應(yīng)的穩(wěn)態(tài)傳熱模型高估了管線中流量的波動幅度。工程應(yīng)用中應(yīng)該采用考慮地層蓄熱效應(yīng)的非穩(wěn)態(tài)傳熱模型對高壓天然氣管流進行動態(tài)模擬。

［1］CHACZYKOWSKI M.Transient flow in natural gas pipeline：the effect of pipeline thermal model［J］.Applied Mathematical Modelling，2010，34(4)：1051－1067.

［2］KIUCHI T.An implicit method for transient gas flows in pipe networks［J］.Int J Heat Fluid Flow，1994，15(5)：378－383.

［3］THORLEY A R D，TILEY C H.Unsteady and transient flow of compressible fluids in pipelines：a review of theoretical and some experimental studies［J］.Int J Heat Fluid Flow，1987，8(1)：3－15.

［4］OSIADACZ A J，CHACZYKOWSKI M.Simulation of non－isothermal transient gas flow in a pipeline［J］.Archives of Thermodynamics，2001，22(1/2)：51－70.

［5］李玉星，姚光鎮(zhèn).輸氣管道設(shè)計與管理［M］.東營：石油大學(xué)出版社，2009：49－51.

［6］陳福來，帥健，馮耀榮，等.高壓天然氣輸送管道斷裂過程中氣體減壓波速的計算［J］.中國石油大學(xué)學(xué)報：自然科學(xué)版，2009，33(4)：130－135.

CHEN Fu－lai，SHUAI Jian，F(xiàn)ENG Yao－rong，et al.Calculation of high－pressure natural gas decompression wave velocity during pipeline fracture［J］.Journal of China U－niversity of Petroleum(Edition of Natural Science)，2009，33(4)：130－135.

［7］呂巍然，王君苓，陳院生，等.一類非齊次微分方程解的增長性［J］.中國石油大學(xué)學(xué)報：自然科學(xué)版，2010，34(2)：166－174.

Lü Wei－ran，WANG Jun－ling，CHEN Yuan－sheng，et al.Growth of solutions of some non－h(huán)omogeneous differential equations［J］.Journal of China University of Petroleum(Edition of Natural Science)，2010，34(2)：166－174.

［8］ISSA R I，SPALDING D B.Unsteady one－dimensional compressible frictional flow with heat transfer［J］.J Mech Eng Sci，1972，14(6)：365－369.

［9］ POLONI M，WINTERBONE D E，NICHOLS J R.The calculation of pressure and temperature discontinuity in a pipe by method of characteristics and the two－step differential Lax－Wendroff methods［J］.ASME FED，1987，62：1－7.

［10］唐建峰，段常貴，呂文哲，等.特征線法在燃氣管道動態(tài)模擬中的應(yīng)用［J］.油氣儲運，2001，20(8)：12－17.

TANG Jian－feng，DUAN Chang－gui，Lü Wen－zhe，et al.Application of the method of characteristics in dynamic simulation of gas pipes［J］.Oil＆ Gas Storage and Transportation，2001，20(8)：12－17.

［11］曹學(xué)文，陳麗，杜永軍，等.超聲速旋流天然氣分離器的旋流特性數(shù)值模擬［J］.中國石油大學(xué)學(xué)報：自然科學(xué)版，2007，31(6)：79－86.

CAO Xue－wen，CHEN Li，DU Yong－jun，et al.Numerical simulation of swirling flow characteristics of supersonic swirling natural gas separator［J］.Journal of China University of Petroleum(Edition of Natural Science)，2007，31(6)：79－86.

［12］GREYVENSTEIN G P.An implicit method for the analysis of transient flows in pipe networks［J］.Int J Numer Meth Eng，2002，53：1127－1143.

［13］GATO L M C，HENRIQUES J C C.Dynamic behaviour of high－pressure natural－gas flow in pipelines［J］.Int J Heat Fluid Flow，2005，26：817－825.

［14］ ABBASPOUR M，CHAPMAN K S.Nonisothermal transient flow in natural gas pipeline ［J］.J Appl Mech Trans ASME，2008，75(3)：181－188.

Effect of underground heat capacity on numerical simulation accuracy of transient flow in natural gas pipelines

DING Yan－peng1，LI Yu－xing1，HAO Hong－na1，BI Yan－jun2
(1.College of Pipeline and Civil Engineering in China University of Petroleum，Qingdao 266555，China;2.China Petro ＆ Gas Pipe Research Institute，Langfang 065000，China)

One－dimensional，nonisothermal gas flow model was adopted to simulate the transient flow in high－pressure gas transmission pipelines.The effects of underground heat capacity on numerical simulation accuracy of fast transient flow in high－pressure natural gas pipelines were studied.From the engineering point of view，the basic variables of the mass，momentum and energy conservation equations were transformed equivalently.One－dimensional axial symmetry unsteady heat transfer model was established between natural gas pipe flow and the outside formation.The dynamic model based on the state equation of BWRS was solved by the method of characteristics.The results show that domestic simplified steady－state heat transfer model overestimates the fluctuation amplitude of flow rate in the pipeline，and the unsteady heat transfer model considering the effect of underground heat capacity should be used to simulate the transient flow in the practical high－pressure gas transmission pipelines.The unsteady－state model and the algorithm have high accuracy.

natural gas pipeline;fast transient flow;numerical simulation;underground heat capacity;method of characteristics

TE 832

A

10.3969/j.issn.1673－5005.2011.03.025

1673－5005(2011)03－0124－05

2010－12－06

國家自然科學(xué)基金項目(51074175);教育部新世紀人才支持計劃(NCET－07－0847)

丁延鵬(1983－)，男(漢族)，山東煙臺人，博士研究生，主要從事輸氣管道數(shù)值模擬研究。

(編輯 沈玉英)