超臨界CO2壓裂井筒流動模型及耦合求解

李小江, 李根生, 王海柱, 田守嶒, 宋先知, 陸沛青, 劉慶嶺

(1.中國石油大學油氣資源與探測國家重點實驗室,北京 102249; 2.中國石化石油工程技術研究院,北京 100101)

CO2可以在較低的溫度(304.2 K)和壓力(7.38 MPa)下達到超臨界態(tài),是一種優(yōu)質的無水壓裂液[1-6]。考慮到超臨界CO2獨特的物理性質[2,7],在其壓裂過程中,井筒熱源項的組成及其對CO2流體流動和傳熱特性的影響尤為重要,需要開展深入研究。目前,針對井筒流動和傳熱問題,主要有解析[8-12]和數(shù)值[13-18]兩種研究方法。筆者在CO2流體物性模型的基礎上,結合流體力學基本方程組,建立考慮不同類型熱量源匯的超臨界CO2壓裂井筒流動解析模型,通過井深和徑向方向的迭代計算,實現(xiàn)模型的雙重耦合數(shù)值求解,并利用該模型研究超臨界CO2壓裂井筒流動與傳熱特性。

1 井筒流動數(shù)學模型

1.1 基本假設

假設:井筒內為穩(wěn)態(tài)傳熱,地層為非穩(wěn)態(tài)傳熱;油管外傳熱僅考慮徑向傳熱;忽略輻射傳熱和相變潛熱;油管、套管與井眼同心;井眼規(guī)則,水泥膠結質量良好,無氣竄。

1.2 流體物性模型

超臨界CO2流體的物性參數(shù)是井筒流動計算的基礎參數(shù),為保證井筒溫度壓力的精確計算,須選用精度較高的流體物性模型,將其與井筒溫度場和壓力場進行耦合計算。為此,采用適用于CO2精度較高的S-W方程[19]計算CO2的密度和比熱容等熱力學性質,采用Fenghour和Vessovic等[20-21]模型計算黏度和導熱系數(shù)等輸運性質。

1.3 壓降模型

以井筒中流體為研究對象,任意截取一段微元體,針對一維穩(wěn)態(tài)流動,質量守恒方程和動量方程[22]為

(1)

(2)

式中,z為井深,m;ν為流速,m/s;g為自由落體加速度,m/s2;θ為井斜角,(°);τw為井壁(油管壁)處剪切應力,Pa;d為油管內徑,m;Ap為油管內截面積,m2。

將質量守恒方程代入動量方程,并對壁面摩擦力項進行代換,可得CO2流體在井筒中向下流動的壓降方程為

(3)

式中,f為Darcy阻力系數(shù)。

阻力系數(shù)f的精確求解是壓力預測的關鍵,筆者采用Wang等[23]提出的專門適用于超臨界CO2摩阻系數(shù)的經(jīng)驗計算公式,據(jù)文獻報道其在常用雷諾數(shù)范圍內的計算誤差小于3.08%。

1.4 傳熱模型

針對開放系統(tǒng),考慮穩(wěn)態(tài)流動的能量守恒方程[22]為

(4)

式中,e為單位質量CO2流體的內能,m2/s2。

井筒與地層之間的熱量傳遞可以用熱量傳遞方程表征[8]為

q=πdU(Tt-Tei)dz.

(5)

式中,q為徑向方向上油管內流體與地層之間傳遞的熱量,J/s;U為總傳熱系數(shù),W/(m2·K);Tt為油管內CO2流體溫度,K;Tei為原始地層溫度,K。

結合質量守恒方程和熱量傳遞方程,可將能量方程化簡為

(6)

式中,h為CO2流體的焓,m2/s2;w為質量流量,kg/s。

考慮單一油管+套管的井身結構,總傳熱系數(shù)表達式為

(7)

式中,ht為油管內CO2流體的強迫對流換熱系數(shù),W/(m2·K);ks為油/套管導熱系數(shù),W/(m·K);dto為油管外徑,m;ha為環(huán)空流體的自然對流換熱系數(shù),W/(m2·K);dci為套管內徑,m;kc為水泥環(huán)導熱系數(shù),W/(m·K);dwb為井眼直徑,m;dco為套管外徑,m;f(t)為地層非穩(wěn)態(tài)傳熱的無因次溫度函數(shù)[10],該函數(shù)是由拉普拉斯變換求解傅里葉導熱定律得來,表示井壁處溫度隨時間的變化關系,可用于表征地層非穩(wěn)態(tài)傳熱過程,而不必對時間域進行離散計算,其具體表達式可參考文獻[10];ke為地層導熱系數(shù),W/(m·K)。

相比常規(guī)流體,CO2對流需要考慮浮升力的影響,因而采用Liao和Zhao[24]提出的溫度壓力適用范圍較廣、精度較高的CO2圓管向下流動強迫對流換熱系數(shù)計算公式。環(huán)空自然對流換熱系數(shù)采用改進的Dropkin和Somerscales[25]關聯(lián)式。通常情況下,油管和套管壁的熱阻可近似忽略不計,而環(huán)空自然對流換熱和水泥環(huán)熱阻則相對比較重要,但是CO2對流換熱系數(shù)的計算需要考慮流體和管壁之間的溫度差異,因而模型不能忽略油管和套管壁熱阻。

考慮比焓梯度方程為

(8)

結合壓降方程(3),可得CO2流體在井筒中向下流動的傳熱方程為

(9)

式中,pfr為CO2流體摩阻損失,Pa;η為焦耳-湯姆遜系數(shù),K/Pa;cp為比定壓熱容，J/(kg·K)。

2 模型雙向耦合求解

超臨界CO2的物性參數(shù)對溫度和壓力的變化非常敏感,進行溫度或壓力某一元素的求解時,必須同時考慮另一元素的變化對物性參數(shù)的影響,為提高計算精度,模型首次采用雙向耦合求解。即在井深方向上,物性模型、壓降模型和傳熱模型三者進行耦合迭代求解;在徑向方向上為了精確求解傳熱模型的關鍵參數(shù)——總傳熱系數(shù),需要同時耦合油管-環(huán)空-地層,以計算不同管壁的內外溫度,并將其同井筒內溫度和壓力一同迭代求解。

考慮到超臨界CO2物性參數(shù)沿井深方向存在顯著差異,將井筒沿井深方向劃分成N個計算單元,計算流程如圖1所示。

圖1 計算流程Fig.1 Flowchart of algorithm

圖1中下角i和i-1代表迭代次數(shù),下標ti、to和ci分別表示油管內壁、油管外壁和套管內壁。

3 模型對比與驗證

為了驗證模型的正確性,從理論模型對比與現(xiàn)場實測數(shù)據(jù)驗證兩個方面對本文中模型的推導過程和模型精度進行驗證。

3.1 模型對比

相比壓降模型,傳熱模型涉及參數(shù)較多,模型更為復雜,因此主要進行傳熱模型的對比分析。

Ramey[8]和Kabir等[10]最早提出并發(fā)展了基于非嚴格能量守恒的井筒傳熱經(jīng)典模型,但這些模型均不考慮相變、摩擦以及氣體膨脹等效應產(chǎn)生的熱量,因而是一種無源匯的、簡化的能量守恒模型。熱量傳遞示意圖如圖2所示。

如圖2所示,模型以長度dz的井筒單元為計算單元,根據(jù)井深方向與徑向方向熱量傳遞守恒的原則,給出井筒熱量平衡方程為

qt(z)-qt(z+dz)=-qtf.

(10)

式中,qt(z)為單位時間內在井深方向流入計算單元的熱量,J/s;qt(z+dz)為單位時間內在井深方向流出計算單元的熱量,J/s;qtf為單位時間內在徑向方向井筒內流體與地層之間的熱量傳遞,J/s。

圖2 熱量傳遞示意圖Fig.2 Sketch map of heat transfer

基于井筒熱量平衡方程(10),推導可得井筒流動的傳熱模型為

(11)

模型推導過程詳見文獻[8]和[10],同時亦可參照本文中2.3進行類似推導。對比式(9)和(11)可以發(fā)現(xiàn),Ramey等的模型中沒有考慮摩擦和氣體膨脹與做功等產(chǎn)生的熱量,因而缺少摩擦生熱項、焦耳-湯姆遜效應項以及氣體膨脹做功項等熱量源匯項。為此,將摩擦生熱源項加入Ramey等的模型中,熱量平衡方程變?yōu)?/p>

qt(z)-qt(z+dz)+qfr=-qtf.

(12)

式中,qfr為單位時間內計算單元內流體與管壁等摩擦而產(chǎn)生的熱量,J/s;類比電流做工生熱(P=UI),摩擦生熱項可用qfr=Δpfrwv[26]表示,其中Δpfr為計算單元的摩阻壓降,Pa;wv為體積流量,m3/s。

同樣,可得考慮摩擦生熱的井筒傳熱模型為

(13)

通過類似方法,氣體做功項和焦耳-湯姆遜效應項亦可添加到Ramey等的熱量平衡方程,最終可以得到與式(9)完全相同的方程。不同的推導方法得到相同的結果,說明本文中模型推導過程的正確性。

3.2 模型驗證

國內超臨界CO2壓裂尚處于起步階段,目前只開展了少數(shù)先導性壓裂試驗,因此以一口超臨界CO2注入井為例[27]進行模型驗證,該井的基本參數(shù)為:井深3 100 m,油管內徑62 mm,油管外徑73 mm,套管內徑124.37 mm,套管外徑137 mm,井筒直徑215.9 mm,地表溫度288.15 K,地溫梯度0.03 K/m,地層密度2 600 kg/m3,地層質量熱容837 J/(kg·K),地層導熱系數(shù)2.09 W/(m·K),油套管導熱系數(shù)44.7 W/(m·K),水泥環(huán)導熱系數(shù)0.52 W/(m·K),排量55.4 t/d,注入壓力24.5 MPa,注入溫度253.15 K,注入時間13 h,水泥返高2 074 m,完井液密度1 000 kg/m3,完井液質量熱容4 186.8 J/(kg·K),完井液導熱系數(shù)0.6 W/(m·K),完井液黏度0.6 mPa·s。

基于現(xiàn)場數(shù)據(jù),采用本文中模型計算得到的井筒溫度與壓力剖面如圖3中實線所示。

由圖3可知,在注入壓力24.5 MPa、注入排量55.4 t/d的條件下,井筒內CO2的流動摩阻很小,井筒內CO2流體壓力隨井深近似呈線性增加,并在井底達到最大值52.46 MPa,與實測值52.02 MPa的相對誤差為0.85%;而文獻[27]的計算壓力值低于本文的計算結果,計算井底壓力為50.06 MPa,與現(xiàn)場實測值的相對誤差為3.77%。

在CO2的注入溫度為253.15 K的條件下,隨著井深的增加,井筒油管內CO2溫度也逐漸增加,且從井口開始,溫度的增速逐漸變緩,直至井深1 000 m以后,溫度增加速度基本保持不變,與地溫梯度相當。這是因為CO2注入排量較小,在井口附近,流體溫度與地層溫度差異較大,因而兩者之間的熱量交換大;隨著井深的增加,溫差逐漸降低,直至達到熱量交換平衡狀態(tài)。計算井底溫度為374.53 K,與實測溫度374 K的相對誤差為0.14%(若溫度單位為攝氏度,則相對誤差為0.52%)。雖然文獻[27]計算得到的井底溫度與本文一致,但其上部井筒溫度預測值顯著高于本文中的計算結果,更接近地層溫度,井筒溫度與地溫的差值隨井深的增加先減小后增大,不符合實際的井筒傳熱規(guī)律,此外對井筒溫度的過高估計也是導致井底壓力計算值低于實測值的原因之一。

此外可知,圖3中的油管內CO2流體與油管內外壁溫度相差很小,但與套管內壁溫度存在一定的差距。上述現(xiàn)象表明油管熱阻要遠小于環(huán)空流體熱阻,在計算井筒溫壓場時應充分考慮環(huán)空流體對井筒傳熱的影響,尤其是在環(huán)空中充滿絕熱性能較好的流體的工況下。

綜上所述,模型計算得到的超臨界CO2壓裂過程中的溫度壓力場,相比于現(xiàn)場實測數(shù)據(jù),井底壓力誤差小于0.85%,井底攝氏溫度誤差小于0.52%,且井筒流動規(guī)律符合實際工況條件下的溫度壓力場變化規(guī)律。由此,可認為模型精度較高,可以開展相應的規(guī)律研究。

圖3 井筒壓力和溫度剖面Fig.3 Wellbore pressure and temperature profile

4 井筒傳熱特性分析

通過與Ramey等經(jīng)典傳熱模型的對比分析可知,Ramey等模型未考慮摩擦生熱、氣體膨脹或壓縮做功以及焦耳-湯姆遜效應等熱量源匯項。重點考察在超臨界CO2小排量壓裂(試壓,或注CO2)和大排量壓裂兩種工況下,不同類型的熱量源匯項對井筒流動和傳熱特性的影響。

圖4為小排量CO2注入壓裂過程中的井筒壓力和溫度剖面,其中注入排量為55.4 t/d,即0.035 m3/min。

由圖4可知,熱量源匯項對井筒壓力基本無影響,但對井筒溫度有一定影響。摩擦生熱對井筒傳熱貢獻很小,圖4中藍線與黑線基本重合;焦耳-湯姆遜效應的影響也較弱;相比之下,氣體做功對井筒傳熱貢獻較大,考慮氣體做功的溫度要高于不考慮氣體做功的溫度,兩者溫度差異在井底達到最大,在本算例中為2 K,相對誤差為0.53 %。一方面由于CO2注入排量較小,因而摩阻較低,摩擦生熱對井筒溫度影響很小;另一方面,由上到下,井筒壓力逐漸增加,對氣體的壓縮效應逐漸增強,因而外界對氣體做功轉化為氣體內能,溫度上升;并且隨著井深的增加,CO2流體溫度逐漸升高,接近地層溫度,導致流體密度降低,可壓縮性增強。因此考慮和不考慮氣體做功兩種情況的溫差逐漸增大,在井底達到最高值。

圖5為大排量CO2壓裂過程中的井筒壓力和溫度剖面,其中排量為3 200 t/d,即1.9 m3/min,注入壓力為60 MPa。

由圖5可知,在CO2大排量壓裂施工過程中,井筒溫壓呈現(xiàn)與CO2小排量注入完全不同的規(guī)律。隨著井深的增加,壓力先增加后降低,摩擦生熱對井筒壓力的影響最大,氣體做功次之,焦耳-湯姆遜效應最小。井筒溫度要顯著低于原始地層溫度,摩擦生熱對井筒溫度的影響也顯著強于其他因素,考慮摩擦生熱的情況同不考慮摩擦生熱相比,井底溫度相差15 K,相對誤差為5.27%;相比CO2小排量注入過程,氣體做功的影響反而減小;同樣,焦耳-湯姆遜效應影響最弱,基本可以忽略不計。

上述現(xiàn)象的原因在于:一方面由于CO2排量較大,流動摩阻損失要高于重力勢能對壓能的轉化部分,壓力出現(xiàn)反轉,而不是持續(xù)增加,因此需要提高井口注入壓力;另一方面,流動摩阻的提高也增加了摩擦生熱產(chǎn)生的熱量,因而考慮摩擦生熱的井筒溫度提高、密度降低,導致實際井筒壓力降低。氣體做功對井筒溫度的影響減弱,因大排量壓裂導致摩阻損失增加,井筒的壓力區(qū)間變化范圍縮小,隨著井深的增加,壓力增速變緩甚至出現(xiàn)反轉,氣體膨脹和壓縮效應減弱并相互抵消,對井筒溫度的貢獻減小。

圖4 小排量下井筒壓力和溫度剖面Fig.4 Wellbore pressure and temperature profile at low flow rate

圖5 大排量下井筒壓力和溫度剖面Fig.5 Wellbore pressure and temperature profile at large flow rate

井底壓力和溫度隨排量的變化如圖6所示。雖然熱量源匯項對井底壓力的影響較小,各曲線重合度較高,但仍可觀察到小排量下不考慮氣體做功曲線及大排量下的不考慮摩擦生熱曲線與其他曲線的差異。隨著排量的增加,氣體膨脹或壓縮做功對井筒溫度的貢獻先增大后減小;摩擦生熱對井筒溫度的貢獻則持續(xù)增加,并在排量為1.37 m3/min時超過氣體做功,成為主控因素;焦耳-湯姆遜效應在研究排量范圍內對井筒溫度貢獻最小,與實際溫度最大相對誤差不超過1.75%。當排量高于0.6 m3/min時,井底溫度開始低于臨界溫度,因此需要采取提高井口注入溫度等措施提高井底溫度,使CO2達到超臨界狀態(tài)。

通過以上分析可知,不同排量下,氣體膨脹或壓縮做功對井筒流動和傳熱的影響較大,在工程計算中不能忽略;大排量超臨界CO2壓裂過程中,摩擦生熱效應顯著增強,不能忽略其對井筒流動和傳熱的影響,這也與文獻[18]的研究結論一致。絕大多數(shù)情況下,焦耳-湯姆遜效應對井筒流動和傳熱的影響較小,在工程上可根據(jù)計算需求選擇性忽略。在現(xiàn)場實際計算時,根據(jù)不同工況忽略影響較小的熱量源匯項,在保證計算精度的同時減少計算量。

圖6 井底壓力和溫度隨排量的變化Fig.6 Variation of bottomhole pressure and temperature with flow rate

5 結 論

(1)在考慮CO2密度等物性參數(shù)變化的基礎上,聯(lián)立質量守恒方程、動量方程和能量守恒方程,建立了考慮熱量源匯項的超臨界CO2壓裂井筒流動模型。通過溫度和壓力以及油管-環(huán)空-地層的耦合迭代計算,實現(xiàn)了井深和徑向方向的雙重耦合數(shù)值求解,得到了超臨界CO2壓裂過程中的井筒溫度和壓力分布。

(2)通過與Ramey經(jīng)典傳熱模型對比分析,驗證了模型推導的正確性;同時計算得到的井底溫度和壓力與現(xiàn)場實測值的誤差均小于1%,能夠滿足工程計算的需要。

(3)熱量源匯項在不同排量下對井筒流動和傳熱的影響程度不同。在CO2小排量注入壓裂過程中,氣體做功對井筒流動和傳熱的影響不能忽略;而在CO2大排量壓裂過程中,需要同時考慮摩擦生熱和氣體做功的影響;不同排量下,焦耳-湯姆遜效應影響較弱,工程計算上可忽略。

參考文獻:

[1] WANG H, LI G, SHEN Z. A feasibility analysis on shale gas exploitation with supercritical carbon dioxide[J]. Energy Sources: Part A, 2012,34(15):1426-1435.

[2] 李根生,王海柱,沈忠厚,等.超臨界CO2射流在石油工程中應用研究與前景展望[J].中國石油大學學報(自然科學版),2013,37(5):76-80, 87.

LI Gensheng, WANG Haizhu, SHEN Zhonghou, et al. Application investigations and prospects of supercritical carbon dioxide jet in petroleum engineering[J].Journal of China University of Petroleum(Edition of Natural Science), 2013,37(5):76-80,87.

[3] 孫寶江,孫文超.超臨界CO2增黏機制研究進展及展望[J].中國石油大學學報(自然科學版),2015,39(3):76-83.

SUN Baojiang, SUN Wenchao. Research progress and prospectives of supercritical CO2thickening technology[J]. Journal of China University of Petroleum(Edition of Natural Science), 2015,39(3):76-83.

[4] MIDDLETON R S, CAREY J W, CURRIER R P, et al. Shale gas and non-aqueous fracturing fluids: opportunities and challenges for supercritical CO2[J]. Applied Energy, 2015,147:500-509.

[5] YOST A B I, MAZZA R L, REMINGTON R E I. Analysis of production response to CO2/sand fracturing: a case study[R]. SPE 29191-MS, 1994.

[6] CHEN Y, NAGAYA Y, ISHIDA T. Observations of fractures induced by hydraulic fracturing in anisotropic granite[J]. Rock Mechanics and Rock Engineering, 2015,48(4):1455-1461.

[7] 程宇雄,李根生,王海柱,等.超臨界二氧化碳噴射壓裂孔內流場特性[J].中國石油大學學報(自然科學版),2014,38(4):81-86.

CHENG Yuxiong, LI Gensheng, WANG Haizhu, et al. Flow field character in cavity during supercritical carbon dioxide jet fracturing[J]. Journal of China University of Petroleum(Edition of Natural Science), 2014,38(4):81-86.

[8] RAMEY H J. Wellbore heat transmission[J]. Journal of Petroleum Technology, 1962,14(4):427-435.

[9] HASAN A R, KABIR C S. Aspects of wellbore heat transfer during two-phase flow[J]. SPE Production & Facilities, 1994,9(3):211-216.

[10] KABIR C S, HASAN A R, KOUBA G E, et al. Determining circulating fluid temperature in drilling, work over, and well control operations[J]. SPE Drilling & Completion, 1996,11(2):74-79.

[11] HAGOORT J. Rameys wellbore heat transmission revisited[J]. SPE Journal, 2004, 9 (4):465-474.

[12] SPINDLER R. Analytical models for wellbore-temperature distribution[J]. SPE Journal, 2011,16(1):125-133.

[13] RAYMOND L R. Temperature distribution in a circulating drilling fluid[J]. Journal of Petroleum Technology, 1969,21(3):333-341.

[14] EICKMEIER J, ERSOY D, RAMEY H, et al. Wellbore temperatures and heat losses during production or injection operations[J]. Journal of Canadian Petroleum Technology, 1970,9(2):115-121.

[15] YOSHIOKA K, ZHU D, HILL A D, et al. Interpretation of temperature and pressure profiles measured in multilateral wells equipped with intelligent completions[R]. SPE 94097-MS,2005.

[16] 王海柱,沈忠厚,李根生.超臨界CO2鉆井井筒壓力溫度耦合計算[J].石油勘探與開發(fā),2011,38(1):97-102.

WANG Haizhu, SHEN Zhonghou, LI Gensheng. Wellbore temperature and pressure coupling calculation of drilling with supercritical carbon dioxide[J]. Petroleum Exploration & Development, 2011,38(1):97-102.

[17] 王瑞和,倪紅堅.二氧化碳連續(xù)管井筒流動傳熱規(guī)律研究[J].中國石油大學學報(自然科學版),2013,37(5):65-70.

WANG Ruihe, NI Hongjian. Wellbore heat transfer law of carbon dioxide coiled tubing drilling[J]. Journal of China University of Petroleum(Edition of Natural Science), 2013,37(5):65-70.

[18] 郭建春,曾冀,張然，等.井筒注二氧化碳雙重非穩(wěn)態(tài)耦合模型[J].石油學報,2015,36(8):976-982.

GUO Jianchun, ZENG Ji， ZHANG Ran, et al. A dual transient coupling model for wellbore of carbon dioxide injection well[J]. Acta Petrolei Sinica, 2015,36(8):976-982.

[19] SPAN R, WAGNER W. A new equation of state for carbon dioxide covering the fluid region from the triple-toint temperature to 1 100 K at pressures up to 800 MPa[J]. Journal of Physical and Chemical Reference Data, 1996,25(6):1509-1596.

[20] VESOVIC V, WAKEHAM W A, OLCHOWY G A, et al. The transport properties of carbon dioxide[J]. Journal of Physical and Chemical Reference Data, 1990,19(3):763-808.

[21] FENGHOUR A, WAKEHAM W A, VESOVIC V. The viscosity of carbon dioxide[J]. Journal of Physical and Chemical Reference Date, 1998,27:31-44.

[22] 汪志明,崔海清,何光渝.流體力學[M],北京:石油工業(yè)出版社,2006.

[23] WANG Z, SUN B, WANG J, et al. Experimental study on the friction coefficient of supercritical carbon dioxide in pipes[J]. International Journal of Greenhouse Gas Control, 2014,25:151-161.

[24] LIAO S M, ZHAO T S. An experimental investigation of convection heat transfer to supercritical carbon dioxide in miniature tubes[J]. International Journal of Heat and Mass Transfer, 2002,45(25):5025-5034.

[25] DROPKIN D, SOMERSCALES E. Heat transfer by natural convection in liquids confined by two parallel plates which are inclined at various angles with respect to the horizontal[J]. Journal of Heat Transfer, 1965,87(1):77-82.

[26] PETERSEN J, BJ?RKEVOLL K S, LEKVAM K, et al. Computing the danger of hydrate formation using a modified dynamic kick simulator[R]. SPE/IADC, 2001.

[27] 張勇,唐人選.CO2井筒壓力溫度的分布[J].海洋石油,2007,27(2):59-64.

ZHANG Yong, TANG Renxuan. Study of distribution of wellbore pressure and temperature of CO2injection well[J]. Offshore Oil, 2007,27(2):59-64.