基于PDE的凍土區(qū)輸氣管道熱影響數(shù)值研究

王翔鵬，潘 振，蘇文濤，劉 利，陳星元，黃紹暢，商麗鳳

(1.遼寧石油化工大學 石油天然氣工程學院，遼寧 撫順 113001；2.撫順市特種設(shè)備監(jiān)督檢驗所，遼寧 撫順 113000)

多年凍土是指溫度等于或低于0℃至少連續(xù)兩年的地面(土壤或巖石，包括冰)[1]。北半球四分之一的面積和17%的裸露地表被永久凍土覆蓋，并且有著大量化石燃料儲量[2]。中國政府已經(jīng)發(fā)布了2030年碳排放高峰指標(CEUP)[3]。天然氣資源比煤炭更清潔，比可再生能源更便宜，因此有很大的發(fā)展空間[4]。管道作為一種安全、可靠、環(huán)保的運輸方式，已廣泛應(yīng)用于天然氣的長途運輸[5]。近幾十年來，化石能源的開發(fā)和運輸導(dǎo)致了凍土區(qū)管道系統(tǒng)的廣泛應(yīng)用。然而，不同工況管道建設(shè)的干擾以及管道與周圍凍土之間的熱相互作用對管道系統(tǒng)的設(shè)計、施工、安全運行和維護構(gòu)成了重大挑戰(zhàn)[6-7]。根據(jù)管道運行的溫度范圍，管道在運行后可能會出現(xiàn)大量凍脹和/或融化沉降。這些問題嚴重影響凍土地區(qū)管道系統(tǒng)的安全運行，不僅增加了施工和維護成本，還增加了管道沿線油氣泄漏的可能性[8]。

本文設(shè)計了基于偏微分方程(PDE)的考慮冰水相變、水分遷移的非飽和凍土區(qū)耦合模型，以西部某凍土區(qū)正溫輸氣管道為例進行模型驗證。對新建管道或停輸再啟動190 d內(nèi)，管-土區(qū)溫度、相對飽和度、體積含水率的分布及變化進行分析。討論了輸氣管道埋深、輸送溫度變化對凍土的熱影響。為了改善這種地面熱狀況，引用了類似道路和鐵路路堤的一種緩解熱影響管道路堤(下文簡稱：管堤)的設(shè)計。管堤主要由碎石層組成。碎石層在寒冷季節(jié)通過自然空氣對流將管道中的熱量傳遞到寒冷環(huán)境中。該層已廣泛應(yīng)用于道路和鐵路路堤的施工中，以控制氣候變暖條件下多年凍土路基的溫度狀況。在冬季，高孔隙碎石層內(nèi)的自然空氣對流可能會增加永久凍土路基的熱損失，或?qū)⒏嗟睦淠軓沫h(huán)境空氣引入凍土路基[9]。為了評估設(shè)計的熱性能，對管道路堤進行模擬，以分析其傳熱過程，為凍土地區(qū)能源管道系統(tǒng)的設(shè)計和施工提供信息參考。

1 模型介紹

1.1 物理模型

圖1顯示了埋地輸氣管道和擬建的管道路堤的物理模型。管道管頂埋深設(shè)置為天然地面以下2.5 m，管道未設(shè)置保溫層，管徑和壁厚分別為660 mm，10 mm。與研究現(xiàn)場實際管線參數(shù)相同[10]。土壤分為三層，即草炭亞粘土(Ⅰ)、亞粘土(Ⅱ)和粉質(zhì)粘土(Ⅲ)。

圖1 埋地輸氣管道和管堤的物理模型

對于管堤的物理模型，路堤由0.8 m厚的碎石層(Ⅳ)，1.6 m厚的粗粒土層組成。為了方便管道維修期間車輛進入，路堤的表面寬度設(shè)計為3.5 m。輸氣管道放置在對凍脹和融沉不敏感的粗顆粒層(Ⅴ)內(nèi)。路基土壤的計算域從路堤坡腳向外延伸8 m，路基土層的地層組成、管徑、壁厚與埋地輸氣管道模型相同。

1.2 控制方程

為了建立非飽和凍土中水運動、熱傳輸和變形的控制方程，進行了以下幾個基本假設(shè)：

(1)土壤顆粒、水和冰是均勻的、各向同性的、不可壓縮的。冰的生長/膨脹是各向同性的。水冰相變是基于平衡熱力學的。

(2)在有輸氣管道的淺層凍土中，假設(shè)始終有未凍水，水分只以液體形式遷移，遵循達西定律，不考慮冰晶的運動。

1.2.1 水分場

凍土中瞬態(tài)非飽和液體流動控制方程[11]

(1)

式中θu——土壤中未凍水的體積含量；

D——水力擴散系數(shù)；

k——土體的滲透系數(shù);

下標w,i——表示水，冰。

水力擴散系數(shù)[11]

(2)

式中c(θu)——比水容量/m-1，由滯水模型確定，本文選用Van Genuehten(VG)滯水模型；

I——阻抗因子[12]，表達式如下

I=10-10θi

(3)

引入一個無量綱含水率變量S(本文稱為“相對飽和度”)，它可由飽和含水率θs和殘余含水率θr定義成標準化含水率的形式[13]

(4)

若殘余含水率θr等于0，相對飽和度S等于土的飽和度。本文θs和θr分別取0.4和0.08。

土體的滲透系數(shù)選用Gardner滲透系數(shù)方程[11]

k(θu)=ks·Sl[1-(1-S1/m)m]2

(5)

式中ks——飽和滲透系數(shù)，9.62×10-5m/s。

比水容量[11]

(6)

式中a0、m、l——隨土質(zhì)變化而變化的本構(gòu)參數(shù)，本文分別取2,1/m；0.5；0.5。

1.2.2 溫度場

考慮冰水相變的凍土傳熱方程如下[14]

(7)

式中ρ——土壤密度/kg·m-3；

ρi——冰的密度/kg·m-3；

C——土壤的體積熱容/J·kg-1·℃-1；

λ——土壤的導(dǎo)熱系數(shù)/W·m-1·℃-1；

θ——體積含水量；

L——相變潛熱/kJ·kg-1,334.56 kJ/kg；

T——溫度/℃；

t——時間/s。

管壁中的熱傳導(dǎo)方程[14]

(8)

式中ρp——鋼管密度/kg·m-3；

Cp——鋼管的體積熱容/J·kg-1·℃-1；

λp——鋼管的導(dǎo)熱系數(shù)/W·m-1·℃-1。

1.2.3 補充方程

因為冰的體積含量θi、未凍水體積含量θu同時存在于控制方程(1)，(7)中，且含有未知變量溫度T。故還需引入一個聯(lián)系方程實現(xiàn)耦合模型求解，本文通過引入凍土孔隙中冰與未凍水之關(guān)系[13]聯(lián)立求解

(9)

式中Tf——土體凍結(jié)溫度/-1 ℃；

B——常數(shù),與土壤種類和含鹽量有關(guān)，本文取0.56。

1.2.4 應(yīng)力場

平衡方程[15]

?{σ}+[F]=0

(10)

幾何關(guān)系[15]

{ε}=?{δ}

(11)

本構(gòu)關(guān)系[15]

{σ}=[C]({ε}-{εv})

(12)

式中F——體積力；

δ——位移量；

C——彈性模量；

εv——有水分相變產(chǎn)生的應(yīng)變，表達式[13]如下

(13)

上述數(shù)值程序由有限元程序COMSOL Multiphysics實現(xiàn)，以初始條件作為初始猜測通過使用的Newton-Raphson算法來尋求解決方案。在每次Newton-Raphson迭代中，使用完全耦合方法中的并行稀疏直接解算器MUMPS求解線性化方程組，得到包含控制方程耦合產(chǎn)生的所有非對角項的系統(tǒng)雅可比矩陣。偏微分方程導(dǎo)入COMSOL的過程可參考文獻[16]。

2 邊界條件與初始參數(shù)

計算模型上表面的溫度邊界基于在研究現(xiàn)場收集的空氣和地溫測量數(shù)據(jù)、附面層理論[9,16-17]。使用表1列出的方程式確定計算域上表面的溫度。輸氣溫度設(shè)為8 ℃，初始地溫設(shè)為年平均地溫-0.9 ℃，下邊界恒溫層溫度設(shè)為-1 ℃，左右邊界設(shè)為絕熱邊界。不考慮水分補給所以水分場邊界設(shè)為零通量。左右邊界設(shè)為輥支撐，地表設(shè)為自由邊界，停輸再啟動或新建管線短期內(nèi)不會有明顯位移，所以管道和底部邊界設(shè)為固定約束。水熱力耦合模型中涉及的參數(shù)[9,17-19]見表2。路基模型中，碎石層設(shè)為多孔介質(zhì)，當外部溫度波動時，空氣在其中自然對流。考慮到冬季的厚積雪，路基的兩個邊坡假設(shè)是不透氣。

表1 計算域上表面的溫度邊界

表2 參數(shù)表

3 模型驗證

為驗證本文模型與數(shù)值計算的準確性與可靠性，按照西部某凍土區(qū)正溫輸氣管道周圍測溫數(shù)據(jù)[10]，以幾何模型中距管道豎直中心線水平距離4.5 m的鉛垂線上溫度分布為研究對象。將12月下旬、3月下旬不同深度位置的土壤自然溫度值與計算值比較。圖2即不同深度位置土壤自然溫度值與數(shù)值模擬值對比圖。從圖2中可見此地區(qū)土壤溫度場分布的計算值與實際測溫元件的實測值吻合良好，測點位置最大偏差不超過1.2℃，說明本模型構(gòu)建以及參數(shù)選取合理。

圖2 實測地溫與模擬值隨深度變化圖

4 結(jié)果與分析

4.1 埋深、輸運溫度對管周凍土溫度場的影響

將天然氣輸送溫度依次設(shè)為8 ℃，10 ℃；管道頂端埋深依次設(shè)為2.5 m，3.0 m進行190 d秋冬季停輸再啟動/新建管線的耦合場模擬分析。運行44 d、89 d以及運行187 d溫度云圖如圖3所示，管道下方土體受管道影響較大，融化范圍(融化溫度0.1℃等值線圖中用粗實線表示)逐漸擴大，融化下限在6 m深處附近，水平融化距離始終在3 m以內(nèi)。受管道影響管道上方凍土先融化未凍水向兩側(cè)凍結(jié)鋒面遷移，兩側(cè)凍土中未凍水結(jié)冰釋放潛熱導(dǎo)致管道上方溫度場不是均勻分布如圖3(c)。距管道中心線4 m以上受輸氣管道影響較小，取距管道中心1 m、1.5 m、2.5 m處鉛垂線上溫度為研究對象。當?shù)乇頊囟茸畹蜁r，2.5 m埋深，不同輸運溫度情況下該鉛垂線上土壤溫度如圖4所示。

圖3 (a)運行44天、(b)運行89天、(c)運行187天時溫度場

由圖4可知當管溫提升時，管道附近1.5 m內(nèi)受影響較大有較明顯升溫，2.5 m處升溫不明顯。1 m地溫最高處高于天然場地5 ℃，凍土原有穩(wěn)定性受到較大影響。不同埋深、相同輸運溫度時，各處溫度分布如圖5。圖5中，當管道埋深增加，熱影響區(qū)明顯下移，1 m地溫最高處從3 m下移至3.6 m。距管道中心線2.5 m處鉛垂深度超過3.6 m時，受地表溫度影響小，埋深變化導(dǎo)致溫度上升，表明埋深增加橫向影響距離變遠。

圖4 不同輸運溫度時土壤各處溫度分布

圖5 不同埋深時土壤各處溫度分布

4.2 相對飽和度S與含水率的變化

非飽和凍土相對飽和度、含水率會隨水分遷移、凍結(jié)而變化。輸氣管道的散熱是否影響水分運動過程，目前實驗研究較少，本文以模擬數(shù)據(jù)分析其變化趨勢。圖6描繪了地表溫度、距管道中心線1.5 m，0.5 m深處溫度及相對飽和度S隨時間的變化。顯然圖6中(1.5,-0.5)點溫度隨氣溫(地表溫度T)變化而變化，且存在20 d左右的滯后。相對飽和度與該點溫度變化趨勢基本一致，隨地表溫度變化而波動。

圖6 地表溫度、(1.5，-0.5)點溫度相對飽和度隨時間變化圖

圖7為管道上方距地表0.5 m處溫度、相對飽和度隨時間圖。降溫初期相對飽和度與地溫變化趨勢有明顯滯后，這與土壤中的水分遷移有關(guān)。由于管道周圍土體融化區(qū)有較多未凍水持續(xù)向上方凍結(jié)區(qū)遷移減緩了該點相對飽和度的下降。同時，通過對比圖6、圖7，管道上方(0,-0.5)點相對飽和度始終大于(1.5,-0.5)。

圖7 (0，-0.5)點溫度相對飽和度隨時間變化圖

輸氣管道上方凍土是易受管道熱影響的區(qū)域，為了進一步分析管道上方水分遷移規(guī)律，取管道上方距地表1.25 m處含水率、該點溫度、地表溫度以及該點下方溫度為因變量，時間為自變量繪制圖8。

圖8 (0,-1.25)點含水率溫度隨時間變化圖

從圖8可以看出，輸運初期淺層土壤中的未凍水逐漸結(jié)冰，該點受管道影響凍土逐漸融化，土壤中的水分由未凍結(jié)區(qū)向凍結(jié)鋒面遷移并發(fā)生聚集[20]，導(dǎo)致含水量下降。管道周圍未凍水向上遷移補充該點水分，所以含水率下降緩慢。然后，該點水向上方更低溫處遷移受環(huán)境溫度影響逐漸凍結(jié)，下方未凍水向上補給，所以含水率下降緩慢。下方土體中水逐漸凍結(jié)，該點水補給降低，同時氣溫逐漸升高地溫最低處逐漸下降，該點未凍水向更低溫處遷移含水率快速下降。該點溫度成為地溫最低處后，含水率逐漸上升。當?shù)乇砩琳郎貢r冰融化，向下補給，含水率迅速上升至趨于飽和。

4.3 埋深、輸運溫度對管周凍土位移場的影響

對于凍土區(qū)管道工程而言，管道附近土體狀態(tài)對管道的穩(wěn)定性有較大影響。圖9描繪了地表溫度最低時，不同埋深、不同管溫影響下管道上方地表與距管道中心線4.5 m處地表的位移圖。圖9中，由于工況的不同，輸送溫度低、埋深大，管道對地表自然凍脹的影響??；管道溫度為8℃埋深2.5 m時，位移差可達14.5 mm。輸送溫度升高，管道上方位移減少，位移差增大，在本文模擬工況中位移差增加了1.5%。埋深增加位移差明顯減小，減小幅度達18%。研究還發(fā)現(xiàn)，當管道埋深增加，溫度場橫向影響距離變遠(見圖5)，遠處未凍水增加，地表凍結(jié)水分向上遷移，引起遠處凍脹位移增加。

圖9 不同條件下位移對比圖

4.4 緩解路堤的熱效應(yīng)

由于路堤目前在實際管道工程中應(yīng)用較少，而在路堤中修建管道時地面熱狀況由環(huán)境溫度、路基、輸氣管道和凍土路基之間復(fù)雜的相互作用決定。以埋地管道模型初始條件運行1年后熱狀態(tài)為初始狀態(tài)，研究第二年的熱的相互作用。使用了第二年初冬、冬季、氣溫回升時，路基土壤的模擬溫度場。第二年初冬氣溫約為-11℃，而管溫高達8℃。由于這兩種溫度差異較大，管堤中的管堤填土層會出現(xiàn)水平方向的溫度傳遞，管道下方的碎石層中的空氣會因受熱不均勻產(chǎn)生對流。如圖10(a)所示，在初冬時，粗顆粒層與碎石層中的水平方向溫度梯度較為明顯，這說明輸氣管道向外擴散的熱量正被管道路堤傳遞到外部環(huán)境中。隨著冬季低氣溫的持續(xù)影響，空氣在碎石層孔隙中不斷運動。隨著熱量傳遞過程的進行，水平方向溫度梯度逐漸消失，粗顆粒層下邊界的溫度與天然地表溫度基本一樣且低于淺層天然凍土路基。管道路堤下部碎石層中有較為明顯的豎直方向溫度梯度如圖10(b)。這表明，管道路堤下方較高溫度優(yōu)先向上傳遞，從而較小對凍土路基的影響。當氣溫回升時，碎石層內(nèi)觀察到明顯的溫度梯度，但此時溫度梯度與冬季時圖10(b)相反。如圖10(c)所示碎石層的上部溫度比下部高，碎石層內(nèi)上層空氣的密度比下層小,分層是穩(wěn)定的，所以碎石層中的傳熱以熱傳導(dǎo)為主。由于碎石層中各級粒徑顆粒的分配不良，熱傳導(dǎo)過程將很緩慢[21]。因此，管道路堤可以減少管道中的熱量向路基土壤中傳遞，并限制管道周圍和路基上方區(qū)域熱損失。

圖10 (a)初冬時管基土壤溫度場;(b)冬季管基土壤溫度場;(c)氣溫回升時管基土壤溫度場

5 結(jié)論

本文針對埋地輸氣管道對周圍凍土的影響作用，基于PDE設(shè)計了埋地管道周圍土壤的熱水力耦合模型，并對模型進行了驗證。使用Comsol Multiphysics模擬了新建管線或停輸再啟動時，周圍土壤的溫度場、水分場、位移場的變化。為改善輸氣管道周圍土壤熱狀況，將路堤應(yīng)用在管道鋪設(shè)中，研究了緩解熱影響管堤的熱效應(yīng)。根據(jù)數(shù)值模擬結(jié)果得出以下結(jié)論：

(1)新建管線或停輸再啟動短期內(nèi)，融化范圍逐漸擴大，融化下限在6 m深處附近，水平融化距離始終在3 m以內(nèi)，當管溫提升時，管道附近1.5 m內(nèi)受影響較大有較明顯升溫，凍土原有穩(wěn)定性受到較大影響；當管道埋深增加，熱影響區(qū)明顯下移，橫向影響距離變遠。氣溫回升短時間內(nèi)，管道上方溫度場不均勻分布。

(2)隨著地表溫度變化，管道附近地表下淺層土壤溫度隨之變化，約有20 d的滯后。相對飽和度與該點溫度變化趨勢基本一致，同一深度靠近管道處的相對飽和度大于遠離管道處。管道上方含水率的變化分為四個階段：緩慢下降、下降、上升、快速上升。

(3)輸氣管道的鋪設(shè)對變飽和凍土的自然凍脹有較大影響，管道溫度為8℃埋深2.5 m時，遠離管道處地表與管道上方位移差可達14.5 mm。輸送溫度升高時，位移差增加幅度較少；埋深增加時，位移差較少幅度較大。

(4)氣溫下降時，管堤可以有效地將輸氣管道向外擴散的熱量傳遞到外部環(huán)境中，并在氣溫回升時，減少凍土路基的熱量攝入。從而達到保護凍土原有穩(wěn)定性的目的。在管道運行的短年限內(nèi)對凍土沒有較大影響。緩解熱影響管堤對凍土路基的熱狀態(tài)控制效果較好。