戚菁菁，初同龍，孫伶

1.中國石油天然氣管道工程有限公司（河北 廊坊 065000）

2.中石油管道局工程有限公司 燃?xì)夥止荆ê颖崩确?065000）

3.中國石油管道分公司 科技處（河北 廊坊 065000）

0 引言

我國所產(chǎn)原油80%以上為高含蠟、高凝點、高黏度的蠟基質(zhì)原油[1]。對于易凝高黏的含蠟原油，由于管道周邊為自然溫度場，環(huán)境溫度相對較低，因此若管道在常溫下運行容易因沿線溫降過快導(dǎo)致沿程摩阻急劇升高而被迫停輸，甚至可能發(fā)生凝管事故。由于管道溫度場的變化會顯著影響熱油管道的最小輸量及安全運行[2]，因此準(zhǔn)確模擬出含蠟原油管道周圍溫度場分布至關(guān)重要。

在計算管道內(nèi)油流溫度場時，常選用一維、二維或三維溫度場計算。管線結(jié)蠟層、管壁、防腐層及土壤之間傳熱在物理上是復(fù)雜的三維問題，然而軸向溫度梯度相對徑向溫度梯度非常小，因此工程上可以忽略軸向溫降影響，這樣三維導(dǎo)熱問題就簡化為二維土壤溫度場模擬問題。文獻[3]將土壤溫度場假設(shè)為一維分布即溫度變化只發(fā)生在縱向上，由于這種方法過于簡單，計算結(jié)果同實際測量的結(jié)果相比有非常大的誤差，因此工程上不予使用。

對于一維、二維和三維溫度場，其計算方法通常包括解析法、實驗方法以及數(shù)值計算等[1-3]。解析法按Fourier 的熱傳導(dǎo)方程進行求解，能夠獲得問題的精確解，并且可以為數(shù)值計算和實驗法提供比較依據(jù)。但是該方法存在很大的局限性，即使求解簡單的導(dǎo)熱問題也相當(dāng)復(fù)雜，對于復(fù)雜幾何形狀的物體和非線性邊界條件下的導(dǎo)熱問題，應(yīng)用解析法幾乎是不可能的。實驗方法是傳熱學(xué)的基本求解方法，通過實驗可以得出相對準(zhǔn)確的結(jié)果，但是此種方法的適應(yīng)性不好，針對某個特定的問題總是要進行新的實驗，耗費的人力、物力和財力巨大，并且實驗周期長。數(shù)值計算在很大程度上彌補了解析法和實驗方法的缺點，其適應(yīng)性強，特別是對于復(fù)雜問題更顯其優(yōu)越性[2-4]。

針對非網(wǎng)格劃分對不規(guī)則區(qū)域具有良好適應(yīng)性及容易實現(xiàn)網(wǎng)格自動化生成等特點[5-6]，本文采用拉普拉斯變換將土壤的半無限大區(qū)域轉(zhuǎn)換為有限矩形區(qū)域，考慮到土壤的物性變化，采用基于非結(jié)構(gòu)化網(wǎng)格離散管線周圍的土壤溫度場。

該類型網(wǎng)格常用的基本單元是三角形和四面體，生成方法包括規(guī)則劃分法、Delaunay 三角剖分法以及陣面推進法等[5-6]。對二維而言，三角形網(wǎng)格可以充滿任意幾何形狀的區(qū)域。雖然其生成過程較為繁瑣，數(shù)據(jù)結(jié)構(gòu)相對復(fù)雜，但具有自適應(yīng)的優(yōu)點。因此選用二維非結(jié)構(gòu)化網(wǎng)格的數(shù)值計算方法模擬計算凍土區(qū)含蠟原油管道溫度場。

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

埋地管道正常運行過程中，傳熱過程包含3 部分，分別為管內(nèi)油流以熱對流方式傳熱至結(jié)蠟層，再通過結(jié)蠟層、管壁、防腐層等將熱量交換至周圍土壤，最后經(jīng)地面與大氣換熱。內(nèi)部因素和外界環(huán)境因素共同構(gòu)成了影響熱油管道散熱的影響因素。其中內(nèi)部因素有油品熱物理性質(zhì)、管道輸量、管徑以及結(jié)蠟層、管壁、防腐層等。外界因素有埋地管段土壤物性、大氣溫度、風(fēng)速等。其物理結(jié)構(gòu)模型如圖1、圖2所示。

根據(jù)簡化，基于圖1、圖2，綜合考慮界面上各層之間的相互影響，建立熱油管道正常運行過程水力-熱力模型。

1.1 流體方程

質(zhì)量守恒方程：

式中：ρ 為密度，kg/m3；A 為管流截面面積，mm2；V 為流體平均速度，m/s；τ為時間，s；z為管道軸向位置，m。

圖1 埋地管道剖面

圖2 埋地管道示意圖

動量守恒方程：

式中：α 為管道軸向與水平方向的夾角，（°）；g 為重力加速度，m/s2；p 為油流截面平均壓力，MPa；f 為摩阻系數(shù)；D為管道內(nèi)直徑，mm。

考慮到油流的速度方向，V2應(yīng)該改為V|V|，不失一般性，本文中予以簡化。

能量方程：

式中：s 為原油比熵，J/（kg·℃）；u 為原油比內(nèi)能，J/kg；h為原油比焓，J/kg；q為單位時間內(nèi)原油在單位管壁面積上的散熱量，J/（mm2·s）。

由質(zhì)量守恒方程、動量守恒方程和能量守恒方程得油流的換熱方程：

式中：Cp為原油定壓比熱容，J/（kg·℃）；T 為原油溫度，℃。

1.2 結(jié)蠟層、管壁和防腐層的傳熱方程

式中：Cj為第j 層（結(jié)蠟層、管壁和防腐層）比熱容，J/K；λj為第j 層（結(jié)蠟層、管壁和防腐層）導(dǎo)熱系數(shù)，kcal/（m·h·℃），其中，j=1、2、3，分別表示結(jié)蠟層、管壁和防腐層，其中1 kcal=4.184 kJ；r 為徑向位置，mm；θ為環(huán)向弧度，rad。

1.3 土壤導(dǎo)熱方程

式中：下標(biāo)為s 的參數(shù)代表土壤的參數(shù)；x 為垂直于軸向的水平位置，m；y為深度，m。

1.4 關(guān)聯(lián)方程

管內(nèi)流體、蠟沉積、管壁、防腐層以及土壤的熱交換過程相互關(guān)聯(lián)，滿足如下方程：

式中：α0為流體對管內(nèi)壁的放熱系數(shù)，W/(m2·℃)；T0為管內(nèi)壁溫度，℃；Rj為各層厚度，mm。

1.5 滿足的邊界條件

埋地管道的界面在物理上是對稱的，其計算區(qū)域也具有對稱性，僅取管道截面右半部分研究，應(yīng)滿足的邊界條件如下：

式中：αa為地表向大氣的放熱系數(shù)，W/(m2·℃)；ho為管道埋深，m；H 為縱向熱力影響區(qū)，m；Ta為大氣溫度，℃；Tn為恒溫層溫度，℃；De為管道外徑，mm。

2 計算區(qū)域離散化及數(shù)值計算方法

2.1 計算區(qū)域離散化

土壤導(dǎo)熱計算區(qū)域應(yīng)為半無限大土壤介質(zhì)區(qū)域。但解析求解需對問題作簡化，因而所獲得的結(jié)果與真實值偏差較大[7]。因此可以通過控制網(wǎng)格劃分，獲得較高精確度的解。本研究的計算區(qū)域限制為原油管道的熱力影響區(qū)，采用數(shù)值計算方法。在數(shù)值計算之前要對求解區(qū)域離散（網(wǎng)格劃分），如圖3—圖5所示。

圖3 土壤區(qū)域密集網(wǎng)格劃分

圖4 土壤區(qū)域稀疏網(wǎng)格劃分

圖5 結(jié)蠟層、管壁和防腐層網(wǎng)格劃分及周圍的土壤網(wǎng)格

采用有限差分法（FDA）對含蠟原油管道進行離散化，如圖6 所示。計算節(jié)點從前一個泵站的出口（點1）到后一個泵站的入口（點n）。

圖6 管道離散圖

2.2 控制方程離散化

對結(jié)蠟層、管壁、防腐層及土壤采用控制容積法研究。

2.2.1 油流方程離散

采用有限差分法，對流體非穩(wěn)態(tài)水力-熱力耦合方程進行離散求解。如圖6所示，沿z方向?qū)⒐芏蝿澐譃槿舾呻x散的小管段，在ΔZ時間段內(nèi)，對方程（4）兩邊進行離散，得到：

2.2.2 對結(jié)蠟層、管壁及防腐層的控制方程離散

結(jié)蠟層、管壁、防腐層中的導(dǎo)熱方程可以寫成如下統(tǒng)一形式：

圖7 極坐標(biāo)中的網(wǎng)格系統(tǒng)

結(jié)合圖7，將式（19）整理成通用的離散化方程形式：

式中：

式中：下標(biāo)E、S、W、N分別代表東、南、西、北4個方位。

2.2.3 土壤控制方程的離散

將計算節(jié)點置于三角形的重心，如圖7所示，節(jié)點P0可視為陰影三角區(qū)域的代表，稱該三角形為P0的控制容積。離散化導(dǎo)熱方程建立起計算節(jié)點P0與相鄰點P1,、P2、P3的溫度的代數(shù)關(guān)系式。土壤導(dǎo)熱方程（6）可以針對任意的控制容積寫成如下形式：

3 算例

選取我國東部地區(qū)某含蠟原油管道，該管道管體大部分處在凍土區(qū)，存在凍脹融沉風(fēng)險，因此需要了解其溫度場變化，從而確定管道安全輸量，減小凍脹融沉風(fēng)險。該管道參數(shù)如下：管頂埋深2.2 m；通過對管道運營的土壤進行試驗分析得到，土壤分層（0～3 m為含水20%亞砂土，3～10 m為含水25%粉質(zhì)黏土，10～20 m為基巖），熱力影響范圍為20 m×15 m；無保溫層。凍土導(dǎo)熱系數(shù)λ1=5.77 kJ/（m·h·℃），融土導(dǎo)熱系數(shù)λ2=5.77 kJ/（m·h·℃），土壤比熱容為C= 1 900 J/( kg·℃) ，土壤表面對大氣的放熱系數(shù)aK= 17. 5 W/ (m2·℃) ，土壤初始溫度T0=-12.5 ℃,下邊界溫度為-5 ℃，地面空氣溫度TK=-15 ℃。假定凍土厚度6 m，根據(jù)管道冬季各站場原油進出站溫度（最低2 ℃左右，最高5 ℃左右）分別假定油溫為2 ℃和5 ℃，預(yù)測油溫下土壤溫度場及凍融情況。

根據(jù)1.1及2.1可得到計算網(wǎng)格劃分如圖8所示。

圖8 網(wǎng)格劃分

采用自主開發(fā)的網(wǎng)格計算軟件，具體參數(shù)見2.1，計算結(jié)果如圖9—圖11所示。管道周圍綠色線所圍為融化區(qū)，即不存在凍土，下方接近水平的線是多年凍土下限，同時給出油溫5 ℃和2 ℃時，0 ℃等值線，從而確定凍土區(qū)域。

圖9 2 ℃油溫條件下土壤溫度場云圖和等值線

圖10 5 ℃油溫條件下土壤溫度場云圖和等值線

圖11 0 ℃等值線

1）油溫5 ℃時，管道熱力影響較大，不存在凍融情況。

2）油溫2 ℃，管道只對周圍10 m 內(nèi)土壤產(chǎn)生熱影響，溫度場零度線在管道外側(cè)，會造成凍土融化。

本文的土壤參數(shù)、地表溫度邊界、現(xiàn)場凍土情況和油溫歷史等條件采用了經(jīng)驗數(shù)據(jù)，在進行數(shù)值模擬后得出的結(jié)果對管道安全運行具有指導(dǎo)作用。

4 結(jié)束語

基于非結(jié)構(gòu)化網(wǎng)格數(shù)值計算方法，針對含蠟原油管道外溫度場建立了二維非穩(wěn)態(tài)熱力計算的數(shù)學(xué)模型。將該模型應(yīng)用于某凍土區(qū)含蠟原油管道溫度數(shù)值模擬，通過模擬得到該管道2 ℃及5 ℃油溫條件下土壤溫度場云圖和等值線圖，了解管道凍融沉降區(qū)域，對于指導(dǎo)管道安全運行和節(jié)能降耗具有重要意義。