基于非結構化網格的含蠟原油管道周圍土壤溫度場模擬

2020-04-11 14:28:48戚菁菁初同龍孫伶

石油工業技術監督 2020年2期

戚菁菁，初同龍，孫伶

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

2.中石油管道局工程有限公司燃氣分公司（河北廊坊 065000）

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

0 引言

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

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

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

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

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

1 數學模型

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

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

1.1 流體方程

質量守恒方程：

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

圖1 埋地管道剖面

圖2 埋地管道示意圖

動量守恒方程：

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

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

能量方程：

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

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

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

1.2 結蠟層、管壁和防腐層的傳熱方程

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

1.3 土壤導熱方程

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

1.4 關聯方程

管內流體、蠟沉積、管壁、防腐層以及土壤的熱交換過程相互關聯，滿足如下方程：

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

1.5 滿足的邊界條件

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

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

2 計算區域離散化及數值計算方法

2.1 計算區域離散化

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

圖3 土壤區域密集網格劃分

圖4 土壤區域稀疏網格劃分

圖5 結蠟層、管壁和防腐層網格劃分及周圍的土壤網格

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

圖6 管道離散圖

2.2 控制方程離散化

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

2.2.1 油流方程離散

采用有限差分法，對流體非穩態水力-熱力耦合方程進行離散求解。如圖6所示，沿z方向將管段劃分為若干離散的小管段，在ΔZ時間段內，對方程（4）兩邊進行離散，得到：

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

結蠟層、管壁、防腐層中的導熱方程可以寫成如下統一形式：

圖7 極坐標中的網格系統

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

式中：

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

2.2.3 土壤控制方程的離散

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

3 算例

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

根據1.1及2.1可得到計算網格劃分如圖8所示。