基于熱-流-固耦合模型的石油鉆井施工過程數值分析

賈善坡，鄒臣頌，王越之，譚賢君，甘新星

（1.長江大學 城市建設學院，湖北 荊州 434023；2.長江大學 油氣資源與勘探技術教育部重點實驗室，湖北 荊州 434023；3.中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071）

1 引 言

近年來針對巖石多物理場耦合特性的研究已成為國內外巖土工程界關注的熱點之一，我國實施與規劃中的深部資源開采、放射性廢料地質深埋處置、石油天然氣地下能源儲存等工程，都涉及到多孔介質中的熱-水-力（THM）耦合問題。Noorishad等基于擴展的Biot固結理論[1]，首次提出了飽和巖土介質的固-液-熱耦合基本方程。在此之后，許多學者對THM耦合理論及數值方法進行了進一步的探討和研究，特別是從20世紀90代開始，美國、英國、加拿大等國以及歐盟的 12個研究小組開始進行 DECOVALEX的國際合作研究計劃，采取不同的方法，對核廢料儲庫圍巖體THM耦合行為進行理論和試驗研究，并取得了矚目的成果[2－3]。在80年代末與 90 年代中期，陸續出現過若干 THM模型的專用程序，如THAMES、MOTIF、FRACON、FEMH、FRIP、FRACTURE、GEORACK等，這些程序大都是專用軟件，算法的適應性有限，在前后處理、單元庫和材料庫等方面不及大型商業軟件[3]；軟件FLAC、UDEC和COMSOL也可用于THM耦合問題計算，但這些程序在求解多種材料組成、復雜的施工過程和三維問題時還存在一定困難。

Rutqvist等[4]將FLAC軟件和TOUGH2軟件相結合，編制了THM耦合程序，并將該程序應用于核廢料地下儲庫的模擬。陳波等[5]推導了固-液兩相介質的溫度場、變形場、滲流場三場耦合問題的有限元格式，開發了三場耦合數值分析軟件 CDST。張玉軍[6]研發出了一個用于分析飽和巖土介質中熱-水-應力耦合彈塑性問題的二維有限元程序，并通過2個算例進行了應用嘗試。馮夏庭等[7]建立了巖體彈性、彈塑性、黏彈塑性的THMC分析模型，開發了數值分析軟件系統，并分析了開挖損傷區溫度-水流-應力-化學耦合作用行為。陳益峰等[8]通過選取位移、水壓、氣壓、蒸氣壓、溫度和孔隙率為基本未知量，建立了有限元計算格式，研發三維八自由度多相流THM全耦合有限元程序THYME3D，并采用膨潤土THM耦合Mock-up試驗對數值模型和計算程序進行驗證。傅少君等[9]根據飽和土體熱固結問題中溫度、滲流及應力之間的復雜耦合關系，系統地探討了二維有限單元方法的求解過程，并在FORTRAN 語言環境下研制了相應的計算機分析程序 APOSE-T。趙延林等[10]在建立雙重介質熱-水-力耦合微分控制方程的基礎上，開發了雙重介質熱-水-力耦合三維計算程序DTHM。

筆者在前人研究的基礎上，推導出巖土介質溫度-滲流-應力耦合的數學模型模型，然后，以MATLAB語言為平臺，對ABAQUS軟件作為求解器，編制THM耦合計算程序，并通過算例對其正確性進行了試算驗證，詳細分析了石油鉆井施工過程中井壁圍巖的熱-流-固耦合作用過程。

2 溫度-滲流-應力耦合數學模型

2.1 巖土介質熱傳遞模型

取任意微元體，設微元體的孔隙度為n，則單位體積中固相骨架所占的空間為1-n，根據能量守恒原理，結合Fourier定律，給出固相骨架的能量守恒方程[4]：

式中：λs、cs、ρs分別為巖土骨架的導熱系數、比熱和密度；qs為單位時間內單位體積骨架產生的能量；T為任意時刻的溫度；T0為初始溫度；βl為線膨脹系數；K為巖土介質體積模量；εv為體積應變。

采用的類似過程可建立流體能量守恒方程。由于流體的滲流速度很小，忽略流體的流動動能，考慮對流和傳導的能量方程可表示：

式中：λf、cf、ρf分別為流體的導熱系數、比熱和密度；qf為單位時間內單位體積流體產生的能量；v為流體的 Darcy滲流速度；(v?)T稱為對流項，表示流體質點移動時引起的溫度變化率。

根據混合物理論[3]，按物質組分比例進行疊加，得到巖土介質的總能量方程：

2.2 巖土介質流固耦合滲流模型

結合Darcy定律，給出考慮介質變形的滲流方程[1－2]：

滲透系數、孔隙度與體積應變的關系[1]為

式中：k0、n0分別為初始滲透系數和孔隙度。

2.3 巖土介質熱流固耦合力學模型

根據Biot有效應力原理，有效應力可表示為

式中：σij為總應力張量；p為孔隙流體壓力，規定以壓力為正；α為Biot 系數；δij為Kronecker符號。

為了描述巖土介質的熱膨脹及塑性變形等行為，總應力可采用增量形式：

式中：De為彈性剛度張量；ε為總應變張量；εp為塑性應變張量；Tε為溫度變化時產生的應變張量。

根據彈塑性理論，式（8）可寫為

式中：

ui為巖土骨架的位移分量。

巖土強度準則采用修正的 Mohr-Coulomb準則，屈服函數F和塑性勢函數G具體表達式[1]為

3 耦合模型的程序實現

3.1 程序結構

考慮到數值收斂和穩定性的要求，本文以解耦的方法開發軟件。熱-流-固耦合模型的求解包括溫度場、流-固耦合場分別求解以及各物理場之間的解耦關系。

雖然，流-固耦合場和溫度場的數值計算有較大差別，但本質上包含了線性化和時步離散（或載荷增量）兩個基本內容，即將流-固體場和溫度場計算按兩個獨立的系統分別設計，通過數據通訊的方式，實現在每一時步上的參數耦合；流-固體耦合是以滲透壓力變化和溫度變化為載荷的，而溫度場的計算又是以流-固體耦合變形為基礎的，在每一時步上不斷修正相關的系數，而這種相互修正又是在一系列時步或載荷上交叉進行的。

以MATLAB為平臺，將ABAQUS做為流-固耦合場和溫度場的求解器，開發巖土介質多場耦合計算軟件，各計算模塊之間數據的存儲和通訊通過編制ABQMAIN子程序來實現，程序流程見圖1。

圖1 多場耦合程序系統程序結構Fig.1 Flowchart of the multifield coupling code

另外，在多場耦合的分析中，要求溫度場的單元劃分是與固體變形計算中的單元劃分是一致的，如平面應變問題可統一取為4節點等參單元。

3.2 熱-固耦合模型測試

為了驗證文中解耦求解方法的正確性和有效性，采用厚壁圓筒熱力耦合模型進行對比驗證，計算模型如圖2所示。厚壁筒彈性模量為300 MPa，泊松比為0.25，熱傳導系數為2.1 W/(m·oC)，比熱為 800 J/(kg·oC)，密度為 2 300 kg/m3，熱膨脹系數為6.0×10-51/ oC，屈服強度為0.8 MPa。模型內半徑為0.5 cm，外半徑為5 cm，內壁溫度為100 ℃，外壁溫度為0 ℃。

圖2 驗證模型有限元網格Fig.2 Finite element meshes of verification model

分如下工況進行對比分析，工況1：熱-彈塑性力學分析，采用 ABAQUS軟件自帶耦合方法。工況2：熱-彈塑性力學分析，溫度對材料力學參數的影響不計，僅考慮塑性應變對熱傳導系數和比熱的影響，采用 ABAQUS軟件自帶耦合方法，并編制USDFLD子程序來實現熱力學參數的動態演化。工況3：熱-彈塑性力學分析，溫度對材料力學參數的影響不計，僅考慮塑性應變對熱傳導系數和比熱的影響，采用本文推薦的方法。

假設工況2和工況3中材料的熱傳導系數和比熱均按如下公式動態演化：

計算結果如圖 3、4所示。較工況 1和工況 2可以看出，變形對熱力學參數的影響對計算結果有較大影響，在熱力完全耦合模型中，應當考慮材料的變形對熱力學參數（如熱傳導系數、比熱等）的影響，但是，從目前的文獻調研發現，關于變形對熱力學參數影響的研究還不多見。比較工況2和工況3可以發現，采取本文推薦的方法時，在經過5次迭代后計算結果就達到了收斂，并且計算結果與ABAQUS自帶的計算模塊結果完全一致，說明本文方法的有效性和可靠性。

圖3 工況1和工況2計算結果比較Fig.3 Comparison of results between cases 1 and 2

圖4 工況2和工況3計算結果比較Fig.4 Comparison of results between cases 2 and 3

3.3 熱-流-固耦合模型測試

以一維熱結問題為例，進行了有限元數值解與解析解的對比[11]。

巖土體的力學、熱學和水力基本參數：彈性模量為600 kPa，泊松比為0.3，土顆粒的體積模量為20 GPa，水的體積模量為5 GPa，土顆粒的熱膨脹系數為1.5×10-5/℃，水熱膨脹系數為2.0×10-4/℃，孔隙率為0.4，土顆粒的比熱為0.8 J/(g?℃)，水的比熱為4.2 J/(g?℃)，水的密度為1.0×106g/m3，土顆粒的密度為 2 600 kg/m3，熱傳導系數為 0.5 W/(m?℃)，Biot系數為1.0，土的滲透系數為1.954 7×10-9m/s。

計算模型如圖5所示，土柱高1 m，在土柱的頂面向下施加載荷100 kPa，初始溫度為10 ℃，初始孔隙水壓力為100 kPa，頂面受60 ℃的溫度載荷作用。模型的邊界條件為：左右兩側為不透水邊界，x方向位移約束；下邊界為不透水邊界，y方向固定；上表面為滲流自由邊界。

圖5 一維熱固結計算模型Fig.5 One-dimensional thermal consolidation model

為研究的方便，定義時間因數T（無量綱）：

式中：K為巖土介質的熱傳導系數；t為時間（s）；h為土柱的高度；n為土的孔隙度；ρs、ρw分別為土顆粒和水的密度；cs、cw分別為土的比熱。

計算結果如圖 6～8所示。對于溫度而言，有限元解與解析解完全吻合；對于節點孔隙壓力和位移來說，有限元解與解析解大致吻合，而其中存在一定的差別是由于解析解是嚴格的一維固結問題，而有限元解是二維固結問題。

算例表明，本文采用解耦的數值方法處理熱-力耦合問題時比較成功的，本程序成功地和現有的大型有限元軟件模塊鏈接，利用現有的成熟的軟件資源，實現復雜的多場耦合計算。

4 模型應用

4.1 計算模型

圖6 測點B溫度變化曲線比較Fig.6 Comparison of temperatures between numerical and analytical solutions of point B

圖7 測點A固結位移變化曲線比較Fig.7 Comparison of displacements between numerical and analytical solutions of point A

圖8 測點B孔隙壓力變化曲線比較Fig.8 Comparison of pore pressures between numerical and analytical solutions of point B

以準噶爾盆地某區塊 3 000 m處泥巖地層為例，取準三維模型進行研究。井眼半徑為0.108 m，上覆巖體平均密度為2 350 kg/m3，水平最大地應力為75 MPa，水平最小地應力為54 MPa，初始溫度為120 ℃，孔隙壓力為45 MPa。該地層的平均機械鉆速為10 m/h，地層鉆開后，井壁位移無約束，液柱壓力為50 MPa，鉆井液溫度為70 ℃，泥餅質量較好，不考慮井壁的滲透性。計算參數見表1。

計算分為如下幾個步驟：①初始地應力場平衡計算；②鉆孔開挖卸載，將鉆孔單元殺死；③熱-流-固耦合計算。

表1 計算參數表Table 1 Main parameters for calculation

4.2 井壁失穩系數

基于Mohr-Coulomb準則和最大張應力理論，定義井壁坍塌系數fb和井壁破裂系數ff分別為

式中：c、φ、σt分別為巖層的黏聚力、內摩擦角和抗拉強度；σ1′、σ3′分別為最大和最小有效主應力。

當fb≤1時，井壁發生坍塌破壞，且fb越小井壁破壞越嚴重；當fb＞1時，井壁處于穩定狀況。

5.3 結果分析與討論

圖9為井壁溫度分布圖。由圖可見，井壁溫度突然降低使得井壁圍巖內溫度也逐漸降低，溫度沿水平徑向的分布規律與垂直徑向的分布規律基本類似，這是由于本文給出的地層熱力學參數是常數，沒有考慮應力對熱力學參數的影響。圖 10為孔隙壓力分布圖，鉆孔卸載后，孔隙壓力在水平向（最大主應力方向）迅速下降，此后孔隙壓力梯度使得井眼孔隙壓力逐漸增大；井眼鉆開后，孔隙壓力在垂向（最小主應力方向）迅速增加，此后滲流作用使得孔隙壓力逐漸耗散。由于井壁溫度降低了 50℃，孔隙壓力產生了一個壓力降（與恒溫相比），24 h后，井壁孔隙壓力與遠場孔隙壓力的壓差約為 8 MPa。

圖 11為井壁徑向應力分布圖。對于水平向位置，由于鉆孔卸載，井壁附近徑向應力迅速降低，此后，溫度的降低引起圍巖產生收縮，加上孔隙壓力梯度引起的滲流作用，使得井壁繼續徑向應力繼續降低；對于垂向位置，卸載后井壁徑向應力迅速降低，此后，由于非均勻地應力、井壁圍巖的收縮變形以及滲流作用的綜合影響，徑向應力逐漸增加，由于井壁內巖石的收縮而拉外圈的巖石，使得井壁內壁的應力小于其內部的應力。

圖9 THM耦合條件下井壁溫度分布示意圖Fig.9 The temperature distributions with different times under THM coupling conditions

圖10 THM耦合條件下井壁孔隙壓力分布示意圖Fig.10 The pore pressure distributions with different cases under THM coupling conditions

圖11 THM耦合條件下井壁徑向應力分布示意圖Fig.11 The radial stress distributions with different cases

圖12 THM耦合條件下井壁環向應力分布示意圖Fig.12 The tangential stress distributions with different cases under THM coupling conditions

圖 12為井壁環向應力分布圖。對于水平向位置，井眼鉆開后，距離井壁0.5rw處應力最大。由于井壁內圍巖溫度的降低引起，使得圍巖產生收縮作用而使得內部應力有增大的趨勢；對于垂向位置，卸載后圍巖環向應力迅速增加，此后，井壁的收縮變形以及滲流作用使得環向應力繼續增大。通過改變鉆井液的溫度值，研究溫度波動對井壁穩定性的影響，具體計算結果見圖13、14。

圖13 地層受井眼內流體加熱作用時井壁失穩系數比較Fig.13 Comparison of coefficients of wellbore stability under fluid heating condition

由圖13可見，地層受井眼內流體加熱作用時，井壁失穩系數比較圖，溫差越大，而井壁坍塌系數和破裂系數越小，表明井壁發生剪切破壞和張性破裂的可能性越大。由圖14地層受井眼內流體流體冷卻時，井壁 失穩系數比較圖可以發現，溫差越大，井壁坍塌系數和井壁破裂系數也越大，表明井壁發生剪切破壞和張性破裂的可能性越小。

5 結 語

在有關文獻的基礎上，利用有效應力原理、能量守恒定律和達西定律推導出多孔介質溫度場、滲流場和應力場耦合的數學模型及其控制方程，以MATLAB語言為平臺，將ABAQUS程序作為一個模塊嵌入迭代算法程序中，提出了解決熱-水-力耦合問題的解耦計算方法，編制了多場耦合有限元程序，并給出了算例驗證，研究了油氣鉆井過程中的熱-流-固耦合作用過程，詳細分析了溫度波動對井壁穩定性的影響。結果表明，在商業軟件ABAQUS高效的內核求解器的輔助下，可以實現大規模多場耦合計算，本文提出的方法能解決熱-流-固耦合問題的復雜程度，完全取決于滲流-應力耦合模塊，可應用于多場耦合條件下巖體黏-彈塑性分析以及損傷分析。

圖14 地層受井眼內流體冷卻作用時井壁失穩系數比較Fig.14 Comparison of coefficients of wellbore stability under fluid cooling condition

[1] NOORISHAD J.TSANG C F.WITHERSPOON P A.Coupled thermal-hydraulic-mechanical phenomena in saturated fractured porous rocks: Numerical approach[J].Journal of Geophysical Research, 1984, 89: 10365－10373.

[2] 陳衛忠, 譚賢君, 伍國軍, 等.非飽和巖石溫度-滲流-應力耦合力學模型研究[J].巖石力學與工程學報, 2007,26(12): 2395－2403.CHEN Wei-zhong, TAN Xian-jun, WU Guo-jun, et al.A thermo-hydro-mechanical coupling model for unsaturated porous rock[J].Chinese Journal of Rock Mechanics and Engineering, 2007, 26(12): 2395－2403.

[3] 周創兵, 陳益峰, 姜清輝, 等.復雜巖體多場廣義耦合分析導論[M].北京: 中國水利水電出版社, 2008.

[4] RUTQVIST J, WU Y S, TSANG C F, et al.A modelingapproach for analysis of coupled multiphase fluid flow,heat transfer, and deformation in fractured porous rock[J].International Journal of Rock Mechanics & Mining Sciences, 2002, 39: 429－442.

[5] 陳波, 李寧, 糕瑞花.多孔介質的變形場-滲流場-溫度場耦合有限元分析[J].巖石力學與工程學報, 2001,20(4): 467－472.CHEN Bo, LI Ning, ZHUO Rui-hua.Finite element analysis of fully coupled thermo-hydro-mechanic behavior of porous media[J].Chinese Journal of Rock Mechanics and Engineering, 2001, 20(4): 467－472.

[6] 張玉軍.熱-水-應力耦合彈塑性二維有限元程序的開發與應用嘗試[J].巖土力學, 2005, 26(2): 169－174.ZHANG Yu-jun.Development and use try of 2D FEM code for coupled thermo-hydro-mechanical elastoplastic analysis[J].Rock and Soil Mechanics, 2005, 26(2): 169－174.

[7] 馮夏庭, 潘鵬志, 丁梧秀, 等.結晶巖開挖損傷區的溫度-水流-應力-化學耦合研究[J].巖石力學與工程學報,2008, 27(4): 656－663.FENG Xia-ting, PAN Peng-zhi, DING Wu-xiu, et al.Thermo-hydro-mechano-chemical coupling studies on excavation damage zone in crystalline rocks[J].Chinese Journal of Rock Mechanics and Engineering, 2008,27(4): 656－663.

[8] 陳益峰, 周創兵, 童富果, 等.多相流傳輸 THM 全耦合數值模型及程序驗證[J].巖石力學與工程學報, 2009,28(4): 649－665.CHEN Yi-feng, ZHOU Chuang-bing, TONG Fu-guo, et al.A numerical model for fully coupled THM processes with multiphase flow and code validation[J].Chinese Journal of Rock Mechanics and Engineering, 2009, 28(4): 649－665.

[9] 傅少君, 吳秋軍, 黃振科.考慮溫度效應的固結問題有限元分析[J].土木工程學報, 2009, 42(1): 95－100.FU Shao-jun, WU Qiu-jun, HUANG Zhen-ke.Finite element analysis of thermal consolidation[J].China Civil Engineering Journal, 2009, 42(1): 95－100.

[10] 趙延林, 王衛軍, 趙陽升, 等.雙重介質熱-水-力三維耦合模型及應用[J].中國礦業大學學報, 2010, 39(5):709－715.ZHAO Yan-lin, WANG Wei-jun, ZHAO Yang-sheng, et al.3D dual medium model of thermal-hydro-mechanical coupling and its application[J].Journal of China University of Mining & Technology, 2010, 39(5): 709－715.

[11] 白冰.巖土顆粒介質非等溫一維熱固結特性研究[J].工程力學, 2005, 22(5): 186－191.BAI Bing.One-dimensional thermal consolidation scharacteristics of geotechnical media under nonisothermal condition[J].Engineering Mechanics, 2005,22(5): 186－191.