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

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

1 引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中：

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

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

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

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

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

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

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

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

另外，在多場(chǎng)耦合的分析中，要求溫度場(chǎng)的單元?jiǎng)澐质桥c固體變形計(jì)算中的單元?jiǎng)澐质且恢碌?，如平面?yīng)變問(wèn)題可統(tǒng)一取為4節(jié)點(diǎn)等參單元。

3.2 熱-固耦合模型測(cè)試

為了驗(yàn)證文中解耦求解方法的正確性和有效性，采用厚壁圓筒熱力耦合模型進(jìn)行對(duì)比驗(yàn)證，計(jì)算模型如圖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，屈服強(qiáng)度為0.8 MPa。模型內(nèi)半徑為0.5 cm，外半徑為5 cm，內(nèi)壁溫度為100 ℃，外壁溫度為0 ℃。

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

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

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

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

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

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

3.3 熱-流-固耦合模型測(cè)試

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

巖土體的力學(xué)、熱學(xué)和水力基本參數(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。

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

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

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

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

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

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

4 模型應(yīng)用

4.1 計(jì)算模型

圖6 測(cè)點(diǎn)B溫度變化曲線比較Fig.6 Comparison of temperatures between numerical and analytical solutions of point B

圖7 測(cè)點(diǎn)A固結(jié)位移變化曲線比較Fig.7 Comparison of displacements between numerical and analytical solutions of point A

圖8 測(cè)點(diǎn)B孔隙壓力變化曲線比較Fig.8 Comparison of pore pressures between numerical and analytical solutions of point B

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

計(jì)算分為如下幾個(gè)步驟：①初始地應(yīng)力場(chǎng)平衡計(jì)算；②鉆孔開(kāi)挖卸載，將鉆孔單元?dú)⑺溃虎蹮?流-固耦合計(jì)算。

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

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

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

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

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

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

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

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

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

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

5 結(jié) 語(yǔ)

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

圖14 地層受井眼內(nèi)流體冷卻作用時(shí)井壁失穩(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)忠, 譚賢君, 伍國(guó)軍, 等.非飽和巖石溫度-滲流-應(yīng)力耦合力學(xué)模型研究[J].巖石力學(xué)與工程學(xué)報(bào), 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ù)雜巖體多場(chǎng)廣義耦合分析導(dǎo)論[M].北京: 中國(guó)水利水電出版社, 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ì)的變形場(chǎng)-滲流場(chǎng)-溫度場(chǎng)耦合有限元分析[J].巖石力學(xué)與工程學(xué)報(bào), 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)力耦合彈塑性二維有限元程序的開(kāi)發(fā)與應(yīng)用嘗試[J].巖土力學(xué), 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é)晶巖開(kāi)挖損傷區(qū)的溫度-水流-應(yīng)力-化學(xué)耦合研究[J].巖石力學(xué)與工程學(xué)報(bào),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ù)值模型及程序驗(yàn)證[J].巖石力學(xué)與工程學(xué)報(bào), 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é)問(wèn)題有限元分析[J].土木工程學(xué)報(bào), 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)軍, 趙陽(yáng)升, 等.雙重介質(zhì)熱-水-力三維耦合模型及應(yīng)用[J].中國(guó)礦業(yè)大學(xué)學(xué)報(bào), 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].工程力學(xué), 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.