熱固耦合下裂隙產(chǎn)狀對(duì)導(dǎo)熱系數(shù)影響的模擬分析*

渠成堃， 周 輝， 任振群， 程文武

(1.中國(guó)科學(xué)院巖土力學(xué)研究所 巖土力學(xué)與工程國(guó)家重點(diǎn)試驗(yàn)室， 武漢 430071；2.兗州煤業(yè)集團(tuán) 濟(jì)寧三號(hào)煤礦， 山東 濟(jì)寧 273500)

熱固耦合下裂隙產(chǎn)狀對(duì)導(dǎo)熱系數(shù)影響的模擬分析*

渠成堃1， 周 輝1， 任振群2， 程文武2

(1.中國(guó)科學(xué)院巖土力學(xué)研究所 巖土力學(xué)與工程國(guó)家重點(diǎn)試驗(yàn)室， 武漢 430071；2.兗州煤業(yè)集團(tuán) 濟(jì)寧三號(hào)煤礦， 山東 濟(jì)寧 273500)

為了研究裂隙產(chǎn)狀對(duì)巖石導(dǎo)熱系數(shù)的影響，通過(guò)模擬立方定律建立了裂隙產(chǎn)狀與導(dǎo)熱系數(shù)間的關(guān)系.利用多場(chǎng)耦合分析軟件Comsol Multiphysics建立四種不同的裂隙幾何模型，施加溫度邊界條件后，得到了不同裂隙產(chǎn)狀下的模型整體導(dǎo)熱系數(shù).結(jié)果表明：裂隙對(duì)熱傳導(dǎo)有明顯的阻礙效果，增大裂隙開(kāi)度，模型整體導(dǎo)熱系數(shù)呈冪函數(shù)形式減??；增大裂隙傾角，導(dǎo)熱系數(shù)呈線性增長(zhǎng)規(guī)律；減小裂隙數(shù)量，增大裂隙表面粗糙度，導(dǎo)熱系數(shù)呈現(xiàn)增大趨勢(shì).

熱固耦合；高放廢料處置；裂隙；表面粗糙度；立方定律；導(dǎo)熱系數(shù)；多場(chǎng)耦合分析軟件；數(shù)值模擬

近年來(lái)，關(guān)于高放廢料處置過(guò)程中的熱固耦合研究[1-3]成為廣大學(xué)者研究的重點(diǎn).考慮多場(chǎng)耦合作用下的巖石長(zhǎng)期穩(wěn)定性，對(duì)于高放廢料處置有重大意義，然而巖石是一種含有大量裂隙的各向異性體，考慮到裂隙對(duì)巖石整體力學(xué)特性的影響，研究巖石在熱固耦合下的長(zhǎng)期力學(xué)特性絕非易事.

本文對(duì)比滲流場(chǎng)中的立方定律[10-12]，通過(guò)多場(chǎng)耦合分析軟件Comsol Multiphysics建立不同裂隙產(chǎn)狀模型，以此獲得等效熱力學(xué)開(kāi)度與導(dǎo)熱系數(shù)間關(guān)系.

1 數(shù)學(xué)模型的構(gòu)建

1.1 整體導(dǎo)熱系數(shù)計(jì)算

q=-kT

(1)

式中：q為熱通量；k為導(dǎo)熱系數(shù)；T為不同方向的溫度梯度.

由于整體熱通量q不便計(jì)算，而在模型上加入熱邊界條件時(shí)，通過(guò)整體的熱通量與通過(guò)模型某一截面的熱通量相等.因此，在計(jì)算模型總體導(dǎo)熱系數(shù)ktot時(shí)，依據(jù)熱通量守恒公式，即

(2)

式中：k1為完整巖石的導(dǎo)熱系數(shù)；T1為模型一斷面的溫度梯度；Ttot為模型整體溫度梯度.計(jì)算過(guò)程中模型整體的Ttot由整體溫度邊界條件得到，而T1可以由某一截面計(jì)算數(shù)據(jù)得到.由于在取截面時(shí)僅取巖石基質(zhì)，因此可以近似認(rèn)為k1為完整巖石的導(dǎo)熱系數(shù).

1.2 導(dǎo)熱機(jī)理分析

目前研究普遍認(rèn)為，巖體中熱量傳遞是以熱傳導(dǎo)形式進(jìn)行的.因此，在求解多場(chǎng)耦合中溫度場(chǎng)時(shí)，傅里葉熱傳導(dǎo)公式被廣泛應(yīng)用，其表達(dá)式為

(3)

式中：ρ為巖體密度；cp為巖體比熱容；u為熱流的對(duì)流速度；Q為熱量.等號(hào)左邊前兩項(xiàng)為熱傳導(dǎo)項(xiàng)，第三項(xiàng)為熱對(duì)流項(xiàng)，由于巖體主要為基質(zhì)，因此認(rèn)為其對(duì)流速度u為0，即可以忽略對(duì)流項(xiàng)影響.

由于巖體內(nèi)部存在裂隙，熱量通過(guò)裂隙時(shí)，除以熱傳導(dǎo)形式傳播外，還以熱輻射形式從裂隙的上表面輻射到下表面，傳遞形式如圖1所示.

圖1 裂隙巖體熱量傳遞形式Fig.1 Form of heat transfer in fractured rock

考慮到裂隙熱輻射效果，當(dāng)熱量在裂隙中進(jìn)行傳遞時(shí)，需引入熱輻射方程，即

(4)

式中：n為方向向量；ε為表面發(fā)射率；σ為相關(guān)常數(shù)；Tamb為環(huán)境溫度；T0為沿n方向的溫度差.

2 不同裂隙產(chǎn)狀模型的建立及計(jì)算

在Comsol中建立高10cm，寬5cm的長(zhǎng)方形巖石試樣.模型上邊界施加373K(100 ℃)恒定溫度，下邊界溫度為273K(0 ℃).選取完整花崗巖，其導(dǎo)熱系數(shù)為3.49W/(m·K)，花崗巖表面發(fā)射率為0.85.研究中構(gòu)建了四種不同形式的裂隙模型，分別為單裂隙平行板模型、單傾斜裂隙模型、多傾斜裂隙模型及粗糙單裂隙模型.

2.1 單裂隙平行板模型

在模型中間構(gòu)造一條水平裂隙，裂隙開(kāi)度以0.02mm梯度逐漸由0mm增大至2mm，如圖2所示.分析不同開(kāi)度下溫度變化情況，以此研究裂隙開(kāi)度對(duì)導(dǎo)熱系數(shù)的影響.

圖2 2 mm裂隙模型溫度分布圖Fig.2 Temperature distribution in model with 2 mm fracture

根據(jù)圖2及式(2)，計(jì)算不同開(kāi)度下的導(dǎo)熱系數(shù)，得到導(dǎo)熱系數(shù)與裂隙開(kāi)度關(guān)系曲線如圖3所示.

圖3 導(dǎo)熱系數(shù)與裂隙開(kāi)度關(guān)系曲線Fig.3 Relation curve between thermal conductivity and fracture aperture

由圖3可以看出，不同裂隙開(kāi)度下模型整體導(dǎo)熱系數(shù)k隨裂隙開(kāi)度b的非線性變化規(guī)律可表示為

k=1.687 5b-0.359

(5)

由此可見(jiàn)，導(dǎo)熱系數(shù)隨裂隙開(kāi)度變化呈現(xiàn)出與立方定律類(lèi)似的冪函數(shù)關(guān)系.原因在于裂隙開(kāi)度增大，相當(dāng)于在巖石中添加了一層導(dǎo)熱性差的材料，致使導(dǎo)熱系數(shù)明顯減小.當(dāng)裂隙開(kāi)度增大到一定程度后，可認(rèn)為巖石被裂隙分隔成上下兩獨(dú)立部分，此時(shí)熱傳導(dǎo)對(duì)于裂隙上表面至下表面的溫度傳遞貢獻(xiàn)較小，熱量更多是以熱輻射形式傳遞，因而導(dǎo)熱系數(shù)變化趨于平緩.

2.2 單傾斜裂隙模型

在長(zhǎng)方形模型中繪制一水平裂隙，并以15°梯度逐漸增大裂隙傾角至90°，研究一系列不同傾角下模型導(dǎo)熱系數(shù)的變化情況，如圖4所示.

圖4 不同傾角裂隙模型溫度分布圖Fig.4 Temperature distribution in models with different fracture angles

根據(jù)計(jì)算得到不同傾角裂隙溫度梯度分布，取模型某一截面，計(jì)算其溫度梯度加權(quán)平均值，最后結(jié)合式(2)得到導(dǎo)熱系數(shù)k隨裂隙傾角α的變化曲線，如圖5所示.

圖5 導(dǎo)熱系數(shù)與裂隙傾角關(guān)系曲線Fig.5 Relation curve between thermal conductivity and fracture angle

由圖5可知，導(dǎo)熱系數(shù)k與裂隙傾角α近似呈線性關(guān)系，相應(yīng)關(guān)系可表示為

k=0.02α+1.617 3

(6)

隨著裂隙傾角的增大，水平截面上裂隙占有面積逐步減小，裂隙對(duì)于熱傳遞阻礙效果明顯減弱.從模型溫度分布可以看出，以45°傾角為界，當(dāng)裂隙傾角大于45°時(shí)，被裂隙分隔的兩部分巖石呈現(xiàn)出明顯溫度變化；而裂隙傾角為90°時(shí)，可以認(rèn)為模型被裂隙分隔成左右兩獨(dú)立部分，相當(dāng)于兩塊完整巖塊，此時(shí)裂隙不再對(duì)熱傳導(dǎo)產(chǎn)生影響.

2.3 多傾斜裂隙模型

在長(zhǎng)方形模型中將多條傾斜裂隙組合以研究導(dǎo)熱系數(shù)隨裂隙數(shù)量的變化情況，如圖6所示.

圖6 不同裂隙產(chǎn)狀模型溫度分布圖Fig.6 Temperature distribution in models with different fracture occurrences

由圖6可知，模型溫度以裂隙為邊界被劃分為幾個(gè)部分，裂隙的存在對(duì)模型整體溫度分布影響較大.取一溫度變化不大的水平截面，計(jì)算此截面上的溫度梯度，并依據(jù)式(2)計(jì)算不同裂隙產(chǎn)狀下模型導(dǎo)熱系數(shù)變化，如表1所示.

表1 不同裂隙產(chǎn)狀模型的導(dǎo)熱系數(shù)Tab.1 Thermal conductivities in models with different fracture occurrences

結(jié)合圖6以及計(jì)算的導(dǎo)熱系數(shù)可知，隨著裂隙傾角的增大，裂隙對(duì)熱量傳遞阻礙作用逐漸減小.

結(jié)合圖6和表1得到的導(dǎo)熱系數(shù)可以明顯看出，45°裂隙數(shù)量增加，巖塊被裂隙分割數(shù)量增多，每塊巖體上溫度差異較大.通過(guò)計(jì)算結(jié)果分析，45°單裂隙時(shí)，巖塊最底部穩(wěn)定后溫度達(dá)到296 K，而隨著裂隙數(shù)量增加，巖塊底部最大溫度分別為292、289和287 K，表明裂隙對(duì)巖塊中的溫度傳遞有明顯阻礙作用.

2.4 粗糙單裂隙模型

在長(zhǎng)方形模型中繪制網(wǎng)格，并定義圖7中間部分網(wǎng)格為裂隙，從而得到簡(jiǎn)單粗糙裂隙模型，以此研究裂隙粗糙度對(duì)模型導(dǎo)熱系數(shù)的影響規(guī)律.

圖7 不同裂隙粗糙度模型溫度分布圖Fig.7 Temperature distribution in models with different fracture roughness

由圖7可知，粗糙裂隙兩側(cè)溫度差異明顯.同樣取一溫度變化較小的水平截面，利用式(2)計(jì)算不同粗糙度裂隙模型的導(dǎo)熱系數(shù)變化，結(jié)果如表2所示.

分析圖7中溫度變化情況可知，隨著裂隙粗糙度的增加，裂隙處顏色逐漸變淺，裂隙上表面與下表面溫差逐漸減小.由表2可知，隨著裂隙表面粗糙度的增加，巖體整體導(dǎo)熱性明顯提升.裂隙表面粗糙度增加，一定程度上增大了裂隙上表面與下表面接觸率.另外，增加表面粗糙度可理解為在原平行板裂隙中填充了基質(zhì)，顯然巖石基質(zhì)導(dǎo)熱系數(shù)遠(yuǎn)高于空氣，因此，導(dǎo)熱性較平行板模型明顯提升.

表2 不同裂隙粗糙度模型的導(dǎo)熱系數(shù)Tab.2 Thermal conductivities in models with different fracture roughness

3 模擬結(jié)果分析

根據(jù)平行板裂隙模型所得到的計(jì)算結(jié)果，可以看出導(dǎo)熱系數(shù)與裂隙開(kāi)度存在類(lèi)似于立方定律的冪函數(shù)關(guān)系.

在2 mm開(kāi)度的平行板裂隙模型上表面施加80 MPa軸向壓力，假設(shè)巖石基質(zhì)不可壓縮，受壓縮僅為裂隙部分，因此在巖石單軸壓縮過(guò)程中裂隙開(kāi)度b可以表示為

b=2-v

(7)

式中，v為模型軸向位移.代入式(5)，則導(dǎo)熱系數(shù)與裂隙開(kāi)度關(guān)系可表示為

k=1.687 5(2-v)-0.359

(8)

將此導(dǎo)熱系數(shù)k代入Comsol計(jì)算，并將結(jié)果與導(dǎo)熱系數(shù)為常數(shù)情況進(jìn)行對(duì)比，結(jié)果如圖8所示.

圖8 單溫度場(chǎng)與熱固耦合條件下溫度對(duì)比

Fig.8 Comparison in temperatures between single temperature field and thermal-solid coupling condition

由圖8可知，如果在計(jì)算溫度變化時(shí)不考慮外部荷載影響，由于裂隙阻礙了上下表面溫度傳遞，裂隙上表面至下表面溫度會(huì)出現(xiàn)明顯跳躍，兩側(cè)溫度差異較大.如果考慮熱固耦合的情況，由于施加了外部荷載，使得裂隙逐漸閉合，裂隙開(kāi)度減小，因而裂隙上下表面的溫差變小，裂隙巖體整體導(dǎo)熱性增加，這更加符合現(xiàn)場(chǎng)情況.

4 結(jié) 論

本文利用Comsol Multiphysics建立了含裂隙巖樣模型，旨在通過(guò)模擬立方定律建立裂隙開(kāi)度與導(dǎo)熱系數(shù)的關(guān)系，主要得出以下結(jié)論：

1) 導(dǎo)熱系數(shù)隨著平行板裂隙開(kāi)度的增加近似呈冪函數(shù)形式遞減.裂隙對(duì)模型自上而下的熱傳導(dǎo)有明顯阻礙效果，并且裂隙開(kāi)度越大，阻礙效果越明顯.

2) 裂隙傾角增加，模型導(dǎo)熱系數(shù)呈近似線性遞增規(guī)律，且當(dāng)裂隙傾角為90°時(shí)，可認(rèn)為模型被分割為左右兩獨(dú)立部分，其導(dǎo)熱系數(shù)通過(guò)計(jì)算與完整巖體一致.

3) 隨著模型內(nèi)裂隙數(shù)量的增加，模型整體導(dǎo)熱系數(shù)呈明顯遞減趨勢(shì).通過(guò)對(duì)多裂隙模型計(jì)算，一方面對(duì)傾斜裂隙模擬結(jié)果進(jìn)行了印證，另一方面也說(shuō)明了裂隙對(duì)模型整體傳熱效果有阻礙作用.

4) 通過(guò)對(duì)幾種具有不同裂隙產(chǎn)狀模型計(jì)算可知，模型整體導(dǎo)熱系數(shù)隨裂隙粗糙度增加呈遞增趨勢(shì).

5) 對(duì)比了不耦合和耦合情況下裂隙巖體表面溫度分布，在僅考慮溫度作用時(shí)，由于裂隙存在，裂隙上表面與下表面溫度呈現(xiàn)較大差異.而施加軸向荷載后，裂隙在荷載作用下逐漸閉合，因而裂隙巖體的導(dǎo)熱性明顯提升，這比較接近真實(shí)情況.

(LIU Xue-wei，LIU Quan-sheng，LU Chao-bo，et al.Rock mechanics coupled thermo-mechanical(TM) fracture propagation numerical manifold method numerical simulation [J].Chinese Journal of Rock Mechanics and Engineering，2014，33(7)：1432-1441.)

(CAI Guo-qing，ZHAO Cheng-gang，TIAN Hui.Nu-merical simulation of coupled thermos-hydro-mechanical behavior for engineered barriers in high-level radio-active waste disposal [J].Chinese Journal of Geotechnical Engineering，2013，35(Sup1)：1-8.)

[3]Li X L.TIMODAZ：a successful international coope-ration project to investigate the thermal impact on the EDZ around a radioactive waste disposal in clay host rocks [J].Journal of Rock Mechanics and Geotechnical Engineering，2013，5(3)：231-242.

(XUE Luan-luan.A composite element model for coupled seepage-heat transfer of fractured rock mass [J].Rock and Soil Mechanics，2016，37(1)：263-268.)

[5]張永利，曹竹，肖曉春，等.溫度作用下煤體裂隙演化規(guī)律數(shù)值模擬及聲發(fā)射特性研究 [J].力學(xué)與實(shí)踐，2015，37(3)：350-354.

(ZHANG Yong-li，CAO Zhu，XIAO Xiao-chun，et al.Numerical simulation and acoustic emission characte-ristics of the evolution of coal body fissures under diffe-rent temperatures [J].Mechanics in Engineering，2015，37(3)：350-354.)

[6]賈春蘭，朱凱.復(fù)合溫度條件下石灰?guī)r多場(chǎng)耦合裂隙滲透侵蝕試驗(yàn)研究 [J].巖土工程學(xué)報(bào)，2015，37(7)：1307-1312.

(JIA Chun-lan，ZHU Kai.Evolution of rock fracture permeability in coupled processes with variable temperatures [J].Chinese Journal of Geotechnical Engineering，2015，37(7)：1307-1312.)

[7]許增光，楊榮，柴軍瑞.考慮巖體與裂隙水流熱量交換作用的溫度場(chǎng)有限元數(shù)值分析 [J].水資源與水工程學(xué)報(bào)，2014，25(4)：69-72.

(XU Zeng-guang，YANG Rong，CHAI Jun-rui.Numerical analysis of finite element on temperature field in effect of heat exchange of rock and fissure water [J].Journal of Water Resources and Water Engineering，2014，25(4)：69-72.)

[8]張巖，李寧，于海鳴，等.溫度應(yīng)力對(duì)裂隙巖體強(qiáng)度的影響研究 [J].巖石力學(xué)與工程學(xué)報(bào)，2013，32(增刊1)：2660-2668.

(ZHANG Yan，LI Ning，YU Hai-ming，et al.Research on influence of thermal stress on fractured rock mass strength [J].Chinese Journal of Rock Mechanics and Engineering，2013，32(Sup1)：2660-2668.)

[9]于永江，張華，張春會(huì)，等.溫度及應(yīng)力對(duì)成型煤樣滲透性的影響 [J].煤炭學(xué)報(bào)，2013，38(6)：936-941.

(YU Yong-jiang，ZHANG Hua，ZHANG Chun-hui，et al.Effects of temperature and stress on permeability of standard coal briquette specimen [J].Journal of China Coal Society，2013，38(6)：936-941.)

[10]王道遠(yuǎn)，劉剛.監(jiān)測(cè)技術(shù)在深基坑開(kāi)挖工程中的應(yīng)用 [J].地質(zhì)學(xué)刊，2010，34(1)：67-72.

(WANG Dao-yuan，LIU Gang.Application of monitoring technology in excavation engineering of deep foundation [J].Journal of Geology，2010，34(1)：67-72.)

[11]盧占國(guó)，姚軍，王殿生，等.平行裂縫中立方定律修正及臨界速度計(jì)算 [J].實(shí)驗(yàn)室研究與探索，2010，29(4)：14-16.

(LU Zhan-guo，YAO Jun，WANG Dian-sheng，et al.Correction of cubic law and calculation of critical velocity in parallel fractures [J].Laboratory Research and Exploration，2010，29(4)：14-16.)

[12]徐維生，陳興周，李娟娟.巖體粗糙裂隙滲流研究 [J].西安理工大學(xué)學(xué)報(bào)，2010，26(2)：218-222.

(XU Wei-sheng，CHEN Xing-zhou，LI Juan-juan.Research on seepage in coarse fractures of rock mass [J].Journal of Xi’an University of Technology，2010，26(2)：218-222.)

(責(zé)任編輯：鐘 媛 英文審校：尹淑英)

Simulation analysis for influence of fracture occurrence on thermal conductivity under thermal-solid coupling

QU Cheng-kun1,ZHOU Hui1,REN Zhen-qun2,CHENG Wen-wu2

(1.State Key Laboratory of Geomechanics and Geotechnical Engineering,Institute of Rock and Soil Mechanics,Chinese Academy of Sciences,Wuhan 430071,China;2.Jining III Coal Mine,Yanzhou Coal Company,Jining 273500,China)

In order to study the influence of fracture occurrence on the thermal conductivity,the relationship between the fracture occurrence and thermal conductivity was established through stimulating the cubic law.With the multi-field coupling analysis software Comsol Multiphysics,four different fracture geometry models were established.After the thermal boundary conditions were added,the thermal conductivities of whole model under different fracture occurrences were obtained.The results show that the fracture has obvious blocking effect on the thermal conduction.When the fracture aperture increases,the thermal conductivities of whole model obviously decrease in the form of power function.With increasing the fracture angle,the thermal conductivities will linearly increase.With decreasing the fracture number and increasing the fracture surface roughness,the thermal conductivities show an increasing trend.

thermal-solid coupling;disposal of high-level nuclear waste;fracture;surface roughness;cubic law;thermal conductivity;multi-field coupling analysis software;numerical simulation

2016-03-09.

國(guó)家自然科學(xué)基金資助項(xiàng)目(51209085，51579093，51479193).

渠成堃(1989-)，男，山東濟(jì)寧人，博士，主要從事圍巖開(kāi)挖擾動(dòng)區(qū)和巖體多場(chǎng)耦合等方面的研究.

22 17∶40在中國(guó)知網(wǎng)優(yōu)先數(shù)字出版.

http：∥www.cnki.net/kcms/detail/21.1189.T.20161222.1740.032.html

10.7688/j.issn.1000-1646.2017.02.18

TU 45

A

1000-1646(2017)02-0219-06