垂直井筒內(nèi)超臨界蒸汽流動過程中的傳熱特性計算與分析

龐占喜，高 輝，張澤權(quán)，王陸亭

（1.油氣資源與探測國家重點實驗室，北京 102249；2.中國石油大學(xué)（北京）石油工程學(xué)院，北京 102249；3.中國中化集團有限公司中化能源物流有限公司，北京 100031）

中國稠油資源儲量豐富，如何高效開采稠油資源仍是目前面臨的主要難題［1-2］。超臨界蒸汽因其高熱量和特殊性質(zhì)，在熱力采油上具有很大的發(fā)展前景。20 世紀60 年代RAMEY 最先建立了井筒傳熱的數(shù)值計算模型［3］，但是該模型只適用于單相流，而后SAGAR 等對RAMEY 的模型進行改進［4］，使其適用于多相流。1965 年，SATTER 把從井筒到地層間傳遞熱量的過程分成從井筒到水泥層的穩(wěn)定傳熱與從水泥層到地層間的非穩(wěn)定傳熱2 部分，但是假設(shè)蒸汽壓力沒有損失且按照一維徑向流來處理傳熱問題［5］。1972 年，PACHECO 等進一步改進了SATTER 的模型，在計算壓力損失時，對地層中的非穩(wěn)態(tài)傳熱采用二維軸對稱模型［6］。中外已對超臨界蒸汽有一定的研究。1966 年，SHITSMAN 對垂直管內(nèi)超臨界蒸汽的傳熱進行了系統(tǒng)的實驗研究［7］。1970 年，NOWAK 等對超臨界蒸汽在水平圓管內(nèi)的流動換熱特性進行了理論研究［8］。2000 年，張毅等對超臨界蒸汽在注汽井筒內(nèi)的參數(shù)變化進行了模擬計算，但將超臨界蒸汽視為理想氣體［9］。超臨界蒸汽在井筒流動過程中，由于井筒熱損失，其可能變成過熱蒸汽、濕蒸汽甚至熱水。所以，研究超臨界蒸汽在井筒內(nèi)的流動特性及熱傳遞規(guī)律對該流體在稠油油藏?zé)岵芍械膽?yīng)用顯得尤為重要［10-11］。筆者首先采用IAPWS-IF97 工業(yè)標(biāo)準［12］建立了超臨界蒸汽的物理化學(xué)參數(shù)計算方法，在此基礎(chǔ)上建立了一套超臨界蒸汽在垂直井筒內(nèi)的流動過程及熱傳遞計算模型，對超臨界蒸汽在稠油油藏?zé)岵芍械膽?yīng)用與推廣具有重要的理論意義。

1 超臨界蒸汽物理化學(xué)特性

不同溫度、壓力條件下水的狀態(tài)可分為4 個區(qū)（圖1a），分別為固態(tài)區(qū)、過飽和水區(qū)、過熱蒸汽區(qū)和超臨界區(qū)；過飽和水區(qū)與過熱蒸汽區(qū)之間的界限為飽和線，飽和線起始于三相點（溫度為0 ℃，壓力為0.101 MPa），結(jié)束于臨界點（溫度為374.15 ℃，壓力為22.12 MPa）；當(dāng)溫度、壓力同時超過臨界點時，水處于超臨界狀態(tài)。在超臨界區(qū)，因高溫而熱膨脹的水相密度和因高壓而被壓縮的蒸汽密度正好相同時的蒸汽被稱為超臨界蒸汽［13-14］。

圖1 水的分區(qū)及超臨界蒸汽的物化參數(shù)Fig.1 Partition of water regions and physicochemical parameters of supercritical steam

超臨界蒸汽的物理化學(xué)參數(shù)是進行井筒熱傳遞計算的基礎(chǔ)，采用IAPWS-IF97 水和水蒸汽性質(zhì)工業(yè)標(biāo)準及其補充模型［12］，完善超臨界蒸汽的物理化學(xué)參數(shù)計算方法，以便于開展井筒內(nèi)流動規(guī)律及換熱損失的研究。利用超臨界蒸汽的無因次吉布斯自由能方程推導(dǎo)得到熱物性參數(shù)，如超臨界蒸汽的密度（比容的倒數(shù)）和比焓的表達式分別為：

其中：

（1）和（2）式中的無因次吉布斯自由能系數(shù)和指數(shù)來源于對無因次吉布斯函數(shù)的回歸，具體取值見文獻［12］。在此基礎(chǔ)上，基于ITS-90 溫標(biāo)和IAP?WS-IF97 計算得到的密度值，根據(jù)不同的溫度與壓力值可迭代計算出對應(yīng)的運動黏度值，再將運動黏度與密度相乘則可得到動力黏度值。計算得到超臨界蒸汽的比容、比焓、運動黏度等熱物性參數(shù)曲線（圖1b—1d），為下一步進行超臨界蒸汽在井筒中的流動特性及傳熱規(guī)律分析提供基礎(chǔ)數(shù)據(jù)。

2 數(shù)學(xué)模型與驗證

2.1 基本假設(shè)

超臨界蒸汽在井筒的流動過程中，從井筒到地層的傳熱過程可分成2部分（圖2），第一部分是由注汽管中心到水泥環(huán)外緣的一維穩(wěn)定傳熱，第二部分是由水泥環(huán)外緣到地層間的一維不穩(wěn)定傳熱［15］。超臨界蒸汽井筒沿程熱損失計算自井口開始，其基本假設(shè)包括：①井口處超臨界蒸汽的參數(shù)（速率、壓力、溫度）保持不變，地表溫度恒定。②垂直井筒中蒸汽的流動是等質(zhì)量流。③垂直井筒底部使用封隔器，確保蒸汽不竄入油套環(huán)空，油套環(huán)空中充滿導(dǎo)熱系數(shù)恒定的某種流體。④注汽井井筒內(nèi)的蒸汽為一維流動，且同一截面處蒸汽的壓力、溫度相等。⑤不考慮沿井身方向的縱向傳熱。

圖2 超臨界蒸汽流動過程對應(yīng)井筒結(jié)構(gòu)剖面Fig.2 Profile of wellbore structure corresponding to supercritical steam flow

2.2 數(shù)學(xué)模型建立

根據(jù)以上假設(shè)，自垂直井筒中心經(jīng)注汽管內(nèi)管、隔熱管、注汽管外管、油套環(huán)空、套管、水泥環(huán)后再經(jīng)過地層的傳熱過程需要分成2部分分別計算。

2.2.1 自注汽管中心到水泥環(huán)外緣穩(wěn)定傳熱過程

自注汽管中心到水泥環(huán)外緣單位時間內(nèi)dz微元長度上的熱流量表達式為：

2.2.2 自水泥環(huán)外緣至地層間的不穩(wěn)定傳熱過程

計算水泥環(huán)外緣到地層間的不穩(wěn)定傳熱時，引入時間變量，因傳熱量與時間成反比，則自水泥環(huán)外緣至地層間單位時間內(nèi)dz微元長度上的熱流量表達式為：

其中：

F（t）表達式可根據(jù)SATTER計算模型給出［5］：

根據(jù)能量守恒原理，自注汽管中心至水泥環(huán)外緣所傳遞的熱量等于水泥環(huán)外緣至地層間所傳遞的熱量，將（5）式、（6）式合并，可得：

2.2.3 相關(guān)方程

注汽管內(nèi)一維流動的動量方程為：

其中，重力項、摩擦損失項、加速損失項的表達式分別為：

超臨界蒸汽的特殊物理化學(xué)性質(zhì)造成其流動特征完全不同于常規(guī)濕蒸汽，因其處于高溫高壓狀態(tài)下的高速流動，必須對加速度損失項進行修正。有些研究者將超臨界蒸汽視為理想氣體進行計算，勢必造成較大的計算誤差［9，16-18］。有的學(xué)者提出了適用于超臨界蒸汽管流沿程摩擦損失梯度的計算式為［19］：

在超臨界條件下，考慮氣液兩相處于單一相態(tài)的特征，對摩擦損失項化簡變形得：

其中：

（15）—（18）式充分體現(xiàn)了超臨界蒸汽物化性質(zhì)對所受摩擦力的影響。

注汽管內(nèi)一維流動的能量方程為：

經(jīng)推導(dǎo)，得到能量守恒方程為：

計算井筒總熱阻時，除考慮油套環(huán)空內(nèi)的對流換熱與輻射換熱外，超臨界蒸汽在注汽管內(nèi)的高速流動也使得注汽管內(nèi)壁處存在對流換熱效應(yīng)，因此引入注汽管內(nèi)壁處的對流換熱系數(shù)來進一步修正總熱阻。這也是不同于其他文獻中關(guān)于井筒熱損失計算的一種修正。井筒總熱阻的表達式為：

其中的輻射傳熱可表達為：

2.3 求解思路

垂直井筒中超臨界蒸汽在流動過程中的熱損失計算流程包括：①針對井筒選擇合適的井段長度（Δz）、時間步長（Δt）和溫度變化值（ΔT），首先令壓力降（Δp）為0，計算相應(yīng)的超臨界蒸汽的比容、比焓及運動黏度。②采用（10）式計算壓力變化微分量（dp），并計算相應(yīng)的超臨界蒸汽的比容、比焓及運動黏度。③采用（21）式計算井筒總熱阻，再利用（9）式計算水泥環(huán)外緣溫度。④根據(jù)自注汽管中心到注汽管外壁導(dǎo)出的熱量與注汽管外壁到套管內(nèi)壁導(dǎo)出的熱量相等的原理，計算注汽管外壁溫度和套管內(nèi)壁溫度。⑤利用（20）式計算超臨界蒸汽的焓值變化微分量（dHcs），從而得到溫度變化微分量（dT），判斷|dT-ΔT|是否小于誤差ε；若滿足要求則進行步驟⑥，否則令ΔT=dT，返回步驟①重新迭代。⑥選擇下一個井段，按照步驟②—⑤進行迭代計算，直至計算到井底為止。⑦選擇下一個時間段，按照步驟①—⑥進行迭代計算，直至計算注汽結(jié)束為止。

2.4 模型驗證

文獻［9］將超臨界蒸汽處理為理想氣體，文獻［10］中摩擦損失計算采用了分散流型摩擦壓降系數(shù)法，基于文獻［9-10］中所采用的數(shù)學(xué)模型和計算參數(shù)，井口注入溫度為400 ℃，注入壓力為25 MPa，注汽速率選擇12 t/h，再利用（1）式、（2）式及（5）—（22）式計算得到井筒沿程的溫度和壓力。采用的熱物性參數(shù)及井筒尺寸如表2 所示。由計算結(jié)果（圖3）可見，本文計算得到的溫度和壓力變化趨勢與文獻［9-10］中的計算結(jié)果變化趨勢一致；在注汽速率為12 t/h 的情況下，本文模型與文獻［9］模型間的溫度誤差為0.367%，壓力誤差為0.678%；本文模型與文獻［10］模型間的溫度誤差僅為0.018%，壓力誤差僅為0.243%，可見本文結(jié)果與文獻［10］結(jié)果之間的誤差更小。本文與文獻［10］都是基于超臨界蒸汽的物理化學(xué)參數(shù)，均考慮了超臨界蒸汽在井筒中流動時的質(zhì)量守恒、動量守恒與能量守恒三大原理，本文與文獻［10］的差異在于本文對井筒內(nèi)超臨界蒸汽流動時的摩擦損失計算中也充分考慮了超臨界蒸汽的特殊性質(zhì)，然而文獻［9］中將超臨界蒸汽視為理想氣體，其計算結(jié)果與另外2 種模型的計算結(jié)果差異較大。

表2 相關(guān)熱物性參數(shù)及井筒尺寸數(shù)據(jù)Table2 Related thermophysical properties and wellbore dimension

圖3 本文計算結(jié)果與文獻計算結(jié)果驗證對比Fig.3 Comparison between calculation results in this paper and those in literatures

3 影響因素分析

3.1 注汽壓力

井身長度選取1 000 m，注汽時間選取15 d，在定質(zhì)量流量和井口溫度的條件下，由不同注汽壓力時井筒沿程溫度與壓力變化（圖4）可以看出，當(dāng)注入壓力分別為23，24 和25 MPa 時，井筒內(nèi)的壓力均隨井深的增加而逐漸增加，然而溫度卻隨著井深的增加而逐漸下降，并且初始注汽壓力越高，在相同長度井段內(nèi)壓力的增值越大，而溫度降低幅度逐漸減小；當(dāng)注汽壓力為25 MPa時，井深為1 000 m 處的溫度降低幅度僅為0.43%，而當(dāng)注汽壓力為23 MPa時，井深為1 000 m 處的溫度降低幅度卻達1.71%。因此，在相同注汽溫度條件下，升高注汽壓力有利于維持井底較高的蒸汽溫度而使蒸汽保持超臨界狀態(tài)。

圖4 不同注汽壓力時井筒沿程溫度與壓力變化Fig.4 Temperature and pressure changes along the wellbore under different steam injection pressures

3.2 注汽溫度

在定質(zhì)量流量和井口壓力（25 MPa）的條件下，由不同注汽溫度時井筒沿程溫度與壓力的變化（圖5）可知，當(dāng)井口注汽溫度分別為400，425 和450 ℃時，初始注汽溫度越高，溫度下降越快，同時壓力升高越快；而當(dāng)注汽溫度為400 ℃時，井筒沿程溫度隨井深增加而出現(xiàn)略有降低趨勢，井深為1 000 m 處的溫度降低幅度僅為0.45%，而當(dāng)注汽溫度為450 ℃時，井深為1 000 m 處的溫度降低至426.7 ℃，降低幅度達到5.17%，但井底仍保持高溫狀態(tài)?？梢姡岣咦⑵麥囟扔欣诰S持較高的井底溫度而使蒸汽保持超臨界狀態(tài)。

圖5 不同注汽溫度時井筒沿程溫度與壓力變化Fig.5 Temperature and pressure changes along the wellbore at different steam injection temperatures

3.3 注汽速率

注汽溫度選擇400 ℃，注汽壓力選擇25 MPa，研究注汽速率對井筒沿程溫度與壓力的影響。由計算結(jié)果（圖6）可知，當(dāng)注汽速率低于2 t/h 時，超臨界蒸汽在井筒中流動時會發(fā)生相態(tài)轉(zhuǎn)變，相態(tài)轉(zhuǎn)變的深度約為1 000 m，呈現(xiàn)為壓力快速上升而溫度大幅下降的趨勢，即此時的水轉(zhuǎn)變?yōu)槌邷責(zé)崴?，處于過飽和水區(qū)（圖1a）。隨注汽速率的升高，轉(zhuǎn)相點逐漸消失，當(dāng)注汽速率超過4 t/h 時，完全滿足井深為1 000 m的垂直井的井底蒸汽為超臨界狀態(tài)的要求。當(dāng)注汽速率過低時，圖6 中溫度變化曲線的切線斜率隨井深不斷發(fā)生改變，這種現(xiàn)象可根據(jù)超臨界蒸汽的物理化學(xué)特性進行分析。在注汽速率一定的情況下，超臨界蒸汽在井口附近所受壓力小，使得蒸汽密度小，所以在單位時間內(nèi)通過注汽管的注汽速率就大，即蒸汽流速高，使得蒸汽內(nèi)部、蒸汽與注汽管內(nèi)壁之間的能量損失嚴重，因此自井口至200 m 井深范圍內(nèi)的井筒沿程溫度下降快，溫度曲線斜率大；而當(dāng)蒸汽進入深層注汽管后，蒸汽壓力升高，低流速時甚至變?yōu)闊崴疇顟B(tài)，其流動過程中的熱損失速率進一步增加，造成溫度下降幅度進一步加大，溫度曲線斜率增加［15-16］。

圖6 不同注汽速率時井筒沿程溫度與壓力變化Fig.6 Temperature and pressure changes along wellbore at different steam injection rates

3.4 注汽時間

注汽壓力選擇25 MPa，注汽溫度選擇400 ℃，注汽速率選擇6 與8 t/h，注汽時間分別選擇1，15 和30 d，研究注汽時間對井筒沿程溫度變化的影響。由計算結(jié)果（圖7）可見，注汽速率保持不變時，隨著天數(shù)的增加，井筒內(nèi)溫度逐漸上升；在注汽初期，井筒溫度上升速度快，而隨著注汽時間的增加，井筒溫度逐漸趨于穩(wěn)定，當(dāng)注汽時間超過15 d 后，井筒內(nèi)溫度分布基本達到穩(wěn)定狀態(tài)。隨著注汽速率的增加，井筒沿程溫度呈現(xiàn)明顯增加趨勢，即提高注汽速率有助于減少井筒熱損失，有利于維持井底蒸汽處于超臨界狀態(tài)［16］。

圖7 不同注汽時間時井筒沿程溫度的變化Fig.7 Temperature changes along wellbore over different steam injection time

3.5 井深

當(dāng)注汽速率選擇8 t/h 時，隨井深的增加，井筒沿程蒸汽壓力逐漸增加，井筒沿程溫度先降低而后在井深約為850 m 處開始增加，而當(dāng)回升一定范圍后出現(xiàn)了溫度的峰值，至約2 300 m 處又逐漸降低（圖8a）。在溫度曲線的第一個谷底處，蒸汽密度在此處會增大很快，然后趨于平緩，使得井筒中部溫度降低趨于平緩；而當(dāng)蒸汽達到深層之后，此時井筒壓力過大，蒸汽流速減小，能量損失速率進一步擴大，造成溫度曲線斜率再一次增加。分析其原因為：在井筒前段處，起主導(dǎo)作用的是熱傳導(dǎo)和對流換熱，蒸汽溫度逐漸下降；當(dāng)蒸汽進入深層井筒時，井筒壓力不斷上升，此時井筒內(nèi)壓力過大，使得超臨界蒸汽密度增大，蒸汽分子之間產(chǎn)生了相對運動，造成溫度出現(xiàn)回升現(xiàn)象［20］。由3 種不同注汽速率（4，8 和12 t/h）在3 000 m 井筒內(nèi)的溫度變化（圖8b）可知，不同注汽速率條件下井筒溫度的變化趨勢不同，注汽速率較低時（4 t/h）井筒溫度下降幅度大，其原因為注汽速率低，蒸汽體積流量小，單位時間內(nèi)經(jīng)井筒進入地層的熱損失速率大，至一定深度后蒸汽出現(xiàn)轉(zhuǎn)相點而變?yōu)榧兯啵⑵俾蕿? t/h時的相態(tài)轉(zhuǎn)變深度約為2 150 m，而當(dāng)注汽速率超過8 t/h時，在井口至3 000 m深度范圍內(nèi)未出現(xiàn)相態(tài)轉(zhuǎn)變點。因此，低流速時經(jīng)熱傳導(dǎo)和對流換熱的熱損失對超臨界蒸汽的流動特性影響程度更大。而當(dāng)蒸汽速率為12 t/h 時，井筒沿程溫度略有下降后而逐漸增加，說明高注汽速率有助于提高井筒內(nèi)超臨界蒸汽的溫度，并有助于維持其超臨界狀態(tài)。

圖8 井深為3 000 m時井筒沿程壓力和溫度的變化Fig.8 Temperature and pressure changes along wellbore at well depth of 3 000 m

4 結(jié)論

基于超臨界蒸汽的特性，運用IAPWS-IF97 水和蒸汽性質(zhì)工業(yè)標(biāo)準及其補充模型，實現(xiàn)超臨界蒸汽的物理化學(xué)性質(zhì)的數(shù)學(xué)描述。在此基礎(chǔ)上，建立超臨界蒸汽在垂直井筒中流動時的熱損失計算數(shù)學(xué)模型，該模型綜合考慮超臨界蒸汽的物理化學(xué)性質(zhì)、超臨界蒸汽的質(zhì)量守恒、動量守恒與能量守恒以及超臨界蒸汽在井筒中流動時的導(dǎo)熱規(guī)律。

分析了超臨界蒸汽井筒流動與熱損失的影響因素，包括注汽壓力、注汽溫度、注汽速率、注汽時間、井深等。井筒內(nèi)超臨界蒸汽的壓力隨井深的增加而增加，而溫度卻隨著井深的增加而逐漸下降；當(dāng)注汽溫度為400 ℃時在1 000 m 處的溫度降低幅度為0.45%，注汽溫度為450 ℃時溫度降低幅度達到5.17%；升高注汽壓力或增加注汽溫度均有利于維持井底蒸汽的超臨界狀態(tài)。

超臨界蒸汽流動過程中，井筒內(nèi)蒸汽壓力快速上升而溫度大幅度下降，隨井深的增加壓力逐漸增加，而井筒內(nèi)溫度呈現(xiàn)波動特征甚至相態(tài)轉(zhuǎn)變特征。如果注汽速率過低或井筒深度過大，造成超臨界蒸汽轉(zhuǎn)變?yōu)榧兯疇顟B(tài)，處于過飽和水區(qū)；當(dāng)注汽速率為2 t/h時的相態(tài)轉(zhuǎn)變深度約為1 000 m，而注汽速率為4 t/h時的相態(tài)轉(zhuǎn)變深度約為2 150 m。

符號解釋

a——溫度梯度，℃/m；

A——流動面積，m2；

dp——壓力變化微分量，MPa；

dHcs——超臨界蒸汽的焓值變化微分量，kJ/kg；

dQ——單位時間內(nèi)dz微元長度上的熱流微分量，J/s；

dT——溫度變化微分量，℃；

dz——垂直井筒的微元長度，m；

f——阻力系數(shù)；

f0——不考慮浮力時的強迫對流摩擦阻力系數(shù)；

fv——考慮重力加速效應(yīng)的流動摩擦阻力系數(shù)；

g——重力加速度，m/s2，取值為9.8；

Hcs——超臨界蒸汽的比焓，kJ/kg；

i——無因次吉布斯自由能系數(shù)或指數(shù)的取值標(biāo)號；

ics——超臨界蒸汽的質(zhì)量流速，kg/s；

Ii——無因次吉布斯自由能壓力指數(shù)；

Ji——無因次吉布斯自由能溫度指數(shù)；

ni——無因次吉布斯自由能系數(shù)；

p——壓力，MPa；

p*——參考壓力，MPa，取值為1；

pr——無因次壓力；

Δp——壓力降，MPa；

Q——單位時間內(nèi)的熱流量，J/s；

r1——注汽管內(nèi)管半徑，m；

r2——隔熱管內(nèi)半徑，m；

r3——隔熱管外半徑，m；

r4——注汽管外管半徑，m；

rci——套管內(nèi)半徑，m；

rco——套管外半徑，m；

rh——水泥環(huán)外緣半徑，m；

R——氣體常數(shù)，J/（K?mol），其值為8.314；

Re——雷諾數(shù)；

t——注汽時間，s；

T——絕對溫度，K；

T1，T2——油套環(huán)空內(nèi)、外界面溫度，℃；

T*——參考溫度，K，取值為540；

Te——地層溫度，℃；

Th——水泥環(huán)外緣溫度，℃；

Tm——地表溫度，℃；

Tr——無因次溫度；

Ts——蒸汽溫度，℃；

Δt——時間步長，s；

ΔT——溫度變化值，℃；

vcs——超臨界蒸汽的流動速度，m/s；

W——單位時間內(nèi)摩擦力所做功，J/s；

Δz——井段長度，m；

z——井深，m；

α——熱擴散系數(shù)，m2/s；

γ1——注汽管內(nèi)壁處的對流換熱系數(shù)，W/（m2·℃）；

γc——油套環(huán)空內(nèi)對流換熱系數(shù)，W/（m2·℃）；

γr——油套環(huán)空內(nèi)輻射換熱系數(shù)，W/（m2·℃）；

ε——誤差，℃；

?——黑度，一般取0.7；

θw——井筒總熱阻，（m?s?℃）/J；

λcas——套管導(dǎo)熱系數(shù)，W/（m·℃）；

λcem——水泥環(huán)導(dǎo)熱系數(shù)，W/（m·℃）；

λe——地層導(dǎo)熱系數(shù)，W/（m·℃）；

λins——絕熱層材料導(dǎo)熱系數(shù)，W/（m·℃）；

λtub——注汽管導(dǎo)熱系數(shù)，W/（m·℃）；

μcs——超臨界蒸汽的黏度，mPa?s；

μw——水的黏度，mPa?s；

ρcs——超臨界蒸汽的密度，kg/m3；

ρw——水的密度，kg/m3。

σ——輻射系數(shù)，W/（m2·℃4），可取值為0.000 88；

τf——超臨界蒸汽管流沿程摩擦損失梯度，Pa/m。