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

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

1 引 言

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

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

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

2 溫度-滲流-應(yīng)力耦合數(shù)學模型

2.1 巖土介質(zhì)熱傳遞模型

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

式中：λs、cs、ρs分別為巖土骨架的導(dǎo)熱系數(shù)、比熱和密度；qs為單位時間內(nèi)單位體積骨架產(chǎn)生的能量；T為任意時刻的溫度；T0為初始溫度；βl為線膨脹系數(shù)；K為巖土介質(zhì)體積模量；εv為體積應(yīng)變。

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

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

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

2.2 巖土介質(zhì)流固耦合滲流模型

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

滲透系數(shù)、孔隙度與體積應(yīng)變的關(guān)系[1]為

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

2.3 巖土介質(zhì)熱流固耦合力學模型

根據(jù)Biot有效應(yīng)力原理，有效應(yīng)力可表示為

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

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

式中：De為彈性剛度張量；ε為總應(yīng)變張量；εp為塑性應(yīng)變張量；Tε為溫度變化時產(chǎn)生的應(yīng)變張量。

根據(jù)彈塑性理論，式（8）可寫為

式中：

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

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

3 耦合模型的程序?qū)崿F(xiàn)

3.1 程序結(jié)構(gòu)

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

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

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

圖1 多場耦合程序系統(tǒng)程序結(jié)構(gòu)Fig.1 Flowchart of the multifield coupling code

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

3.2 熱-固耦合模型測試

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

圖2 驗證模型有限元網(wǎng)格Fig.2 Finite element meshes of verification model

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

假設(shè)工況2和工況3中材料的熱傳導(dǎo)系數(shù)和比熱均按如下公式動態(tài)演化：

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

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

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

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

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

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

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

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

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

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

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

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

4 模型應(yīng)用

4.1 計算模型

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

圖7 測點A固結(jié)位移變化曲線比較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

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

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

表1 計算參數(shù)表Table 1 Main parameters for calculation

4.2 井壁失穩(wěn)系數(shù)

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

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

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

5.3 結(jié)果分析與討論

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

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

圖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耦合條件下井壁徑向應(yīng)力分布示意圖Fig.11 The radial stress distributions with different cases

圖12 THM耦合條件下井壁環(huán)向應(yīng)力分布示意圖Fig.12 The tangential stress distributions with different cases under THM coupling conditions

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

圖13 地層受井眼內(nèi)流體加熱作用時井壁失穩(wěn)系數(shù)比較Fig.13 Comparison of coefficients of wellbore stability under fluid heating condition

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

5 結(jié) 語

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

圖14 地層受井眼內(nèi)流體冷卻作用時井壁失穩(wěn)系數(shù)比較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] 陳衛(wèi)忠, 譚賢君, 伍國軍, 等.非飽和巖石溫度-滲流-應(yīng)力耦合力學模型研究[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] 周創(chuàng)兵, 陳益峰, 姜清輝, 等.復(fù)雜巖體多場廣義耦合分析導(dǎo)論[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] 陳波, 李寧, 糕瑞花.多孔介質(zhì)的變形場-滲流場-溫度場耦合有限元分析[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] 張玉軍.熱-水-應(yīng)力耦合彈塑性二維有限元程序的開發(fā)與應(yīng)用嘗試[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] 馮夏庭, 潘鵬志, 丁梧秀, 等.結(jié)晶巖開挖損傷區(qū)的溫度-水流-應(yīng)力-化學耦合研究[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] 陳益峰, 周創(chuàng)兵, 童富果, 等.多相流傳輸 THM 全耦合數(shù)值模型及程序驗證[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] 傅少君, 吳秋軍, 黃振科.考慮溫度效應(yīng)的固結(jié)問題有限元分析[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] 趙延林, 王衛(wèi)軍, 趙陽升, 等.雙重介質(zhì)熱-水-力三維耦合模型及應(yīng)用[J].中國礦業(yè)大學學報, 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] 白冰.巖土顆粒介質(zhì)非等溫一維熱固結(jié)特性研究[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.